eda.Rmd

---
title: "Exploratory data analysis"
output:
  workflowr::wflow_html:
    includes:
      in_header: header.html
editor_options:
  chunk_output_type: console
author: "Patrick Schratz"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  fig.retina = 3,
  fig.align = "center",
  fig.width = 6.93,
  fig.height = 6.13,
  out.width = "100%",
  echo = FALSE
)

library("drake")
library("magrittr")
R.utils::sourceDirectory("R")
R.utils::sourceDirectory("code")

# load drake objects
loadd(
  vi_data, nri_data,
  bands_data, coords_vi_nri_clean,

  filter_values,
  filter_info_gain_nbins,

  trees_with_bands,
  trees_with_bands_no_buffer
)
library("DataExplorer")
library("dplyr")
library("ggsci")
library("ggpubr")
library("ggplot2")
library("knitr")
library("purrr")
library("sp")
library("raster")
library("fs")
```

- VI = Vegetation Indices
- NRI = Normalized Ratio Indices
- HR = Hyperspectral reflectances (image bands)

# Datasets {.tabset .tabset-fade}

## VI

### Overview

```{r eda-1 }
intro <- introduce(vi_data)
intro_df <- data.frame(
  "Name" = c(
    "Rows", "Columns",
    "Discrete columns", "Continuous columns", "All missing columns",
    "Missing observations", "Complete Rows",
    "Total observations"
  ),
  "Value" = c(
    format(intro[["rows"]], big.mark = ","),
    format(intro[["columns"]], big.mark = ","),
    format(intro[["discrete_columns"]], big.mark = ","),
    format(intro[["continuous_columns"]], big.mark = ","),
    format(intro[["all_missing_columns"]], big.mark = ","),
    format(intro[["total_missing_values"]], big.mark = ","),
    format(intro[["complete_rows"]], big.mark = ","),
    format(intro[["total_observations"]], big.mark = ",")
  )
)
knitr::kable(intro_df)
```

### Histograms

```{r eda-2 }
plot_histogram(vi_data)
```

### PCA

```{r eda-3 }
# remove response
vi_data_no_defol_no_response <- vi_data
vi_data_no_defol_no_response$defoliation <- NULL

pca_vi <- DataExplorer::plot_prcomp(vi_data_no_defol_no_response,
  variance_cap = 0.95,
  ggtheme = ggpubr::theme_pubclean(),
  geom_label_args = list(size = 2.5, nudge_y = 0.03),
  parallel = TRUE
)
pca_vi = pca_vi[[1]] +
  labs(title = "Target POV: 95%")
pca_vi[[1]]
```

### Corr

```{r eda-4 }
plot_correlation(vi_data)
```

## NRI

### Overview

```{r eda-5 }
intro <- introduce(nri_data)
intro_df <- data.frame(
  "Name" = c(
    "Rows", "Columns",
    "Discrete columns", "Continuous columns", "All missing columns",
    "Missing observations", "Complete Rows",
    "Total observations"
  ),
  "Value" = c(
    format(intro[["rows"]], big.mark = ","),
    format(intro[["columns"]], big.mark = ","),
    format(intro[["discrete_columns"]], big.mark = ","),
    format(intro[["continuous_columns"]], big.mark = ","),
    format(intro[["all_missing_columns"]], big.mark = ","),
    format(intro[["total_missing_values"]], big.mark = ","),
    format(intro[["complete_rows"]], big.mark = ","),
    format(intro[["total_observations"]], big.mark = ",")
  )
)
kable(intro_df)
```

### Histograms

No histograms for NRI -> too many features.

### PCA

```{r eda-6 }
# remove response
nri_data_no_response <- nri_data
nri_data_no_response$defoliation <- NULL

pca_nri <- DataExplorer::plot_prcomp(nri_data_no_response,
  variance_cap = 0.95,
  ggtheme = ggpubr::theme_pubclean(),
  geom_label_args = list(size = 2.5, nudge_y = 0.03),
  parallel = TRUE
)
pca_nri = pca_nri[[1]] +
  labs(title = "Target POV: 95%")
pca_nri
```

### Corr

No correlation plot for NRI -> too many features.

## HR

### Overview

```{r eda-8 }
intro <- introduce(bands_data)
intro_df <- data.frame(
  "Name" = c(
    "Rows", "Columns",
    "Discrete columns", "Continuous columns", "All missing columns",
    "Missing observations", "Complete Rows",
    "Total observations"
  ),
  "Value" = c(
    format(intro[["rows"]], big.mark = ","),
    format(intro[["columns"]], big.mark = ","),
    format(intro[["discrete_columns"]], big.mark = ","),
    format(intro[["continuous_columns"]], big.mark = ","),
    format(intro[["all_missing_columns"]], big.mark = ","),
    format(intro[["total_missing_values"]], big.mark = ","),
    format(intro[["complete_rows"]], big.mark = ","),
    format(intro[["total_observations"]], big.mark = ",")
  )
)
kable(intro_df)
```

### Histograms

```{r eda-9 }
plot_histogram(bands_data)
```

### PCA

```{r eda-10 }
# remove response
bands_data_no_response <- bands_data
bands_data_no_response$defoliation <- NULL

pca_hr <- DataExplorer::plot_prcomp(bands_data_no_response,
  variance_cap = 0.98,
  geom_label_args = list(size = 2.5, nudge_y = 0.03),
  ggtheme = ggpubr::theme_pubclean()
)
pca_hr = pca_hr[[1]] +
  labs(title = "Target POV: 95%")
pca_hr
```

### Corr

```{r eda-11 }
plot_correlation(bands_data)
```

# Custom plots

```{r eda-12 }
vi_data_plot <- vi_data %>%
  mutate(plot = factor(rep(
    c("Laukiz 1", "Laukiz 2", "Luiando", "Oiartzun"),
    c(559, 451, 301, 497)
  )))
```

# Mean defoliation per plot

```{r eda-13 }
mean_defol <- vi_data_plot %>%
  group_by(plot) %>%
  summarise(mean(defoliation)) %>%
  pull(.)
vi_data_plot %>%
  group_by(plot) %>%
  summarise(mean(defoliation))
```

# Coefficient of variation

```{r eda-14 }
cov_defol <- vi_data_plot %>%
  group_by(plot) %>%
  summarise((sd(defoliation) / mean(defoliation)) * 100) %>%
  pull(.)
vi_data_plot %>%
  group_by(plot) %>%
  summarise((sd(defoliation) / mean(defoliation)) * 100)
```

# sd / skewness

```{r eda-15 }
sd_skewness_defol <- vi_data_plot %>%
  group_by(plot) %>%
  summarise(((sd(defoliation) / mean(defoliation)) * 100) / e1071::skewness(defoliation)) %>%
  pull(.)
vi_data_plot %>%
  group_by(plot) %>%
  summarise(((sd(defoliation) / mean(defoliation)) * 100) / e1071::skewness(defoliation))
```

```{r eda-16 }
sd_defol <- vi_data_plot %>%
  group_by(plot) %>%
  summarise(sd(defoliation)) %>%
  pull(.)
vi_data_plot %>%
  group_by(plot) %>%
  summarise(sd(defoliation))
```

```{r defoliation-distribution-plot, dev = c("png", "pdf")}
boxplot_defol <- vi_data_plot %>%
  group_by(plot) %>%
  ggboxplot(
    x = "plot", y = "defoliation", color = "plot",
    add = "jitter", add.params = list(size = "defoliation")
  ) +
  scale_size(range = c(0.5, 0.5)) +
  annotate("text",
    label = expression(bold(atop("n = 559", bar(x) ~ "= 57.23"))), x = 1,
    y = 112, size = 4, colour = "#BC3C29", fontface = 2
  ) +
  annotate("text",
    label = expression(bold(atop("n = 451", bar(x) ~ "= 13.54"))), x = 2,
    y = 112, size = 4, colour = "#0072B5", fontface = 2
  ) +
  annotate("text",
    label = expression(bold(atop("n = 301", bar(x) ~ "= 68.36"))), x = 3,
    y = 112, size = 4, colour = "#E18727", fontface = 2
  ) +
  annotate("text",
    label = expression(bold(atop("n = 497", bar(x) ~ "= 69.07"))), x = 4,
    y = 112, size = 4, colour = "#20854E", fontface = 2
  ) +
  scale_color_nejm() +
  ggpubr::theme_pubr(base_size = 14) +
  theme(legend.position = "none") +
  labs(y = "Total defoliation per tree (%)", x = "Plot")
boxplot_defol
```

# Point density

Mean average distance between all trees in a given plot.

Calculated via `raster::pointDistance(allpairs = TRUE)`.

```{r eda-18 }
plots <- list("Laukiz 1", "Laukiz 2", "Luiando", "Oiartzun")
dist <- map(plots, ~ {
  coords <- coords_vi_nri_clean %>%
    mutate(plot = factor(rep(
      c("Laukiz 1", "Laukiz 2", "Luiando", "Oiartzun"),
      c(559, 451, 301, 497)
    ))) %>%
    filter(plot == .x) %>%
    dplyr::select(-plot)
  points <- SpatialPoints(
    coords = coords,
    # EPSG: 25830
    proj4string = CRS("+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
  )
  distance <- mean(as.dist(pointDistance(points,
    allpairs = TRUE,
    lonlat = FALSE
  )))
})
set_names(dist, plots)
dist_plots <- unlist(dist)
```

```{r eda-19, echo = FALSE}
# file_move("docs/figure/eda.Rmd/defoliation-distribution-plot-1.pdf",
#           "code/98-paper/journal/")
# file_move("docs/figure/eda.Rmd/pca-bands-1.pdf",
# "code/98-paper/journal/")
```

# PCA feature set comparison

Comparing the "proportion of variance explained" (POV) for all feature sets.

```{r pca-pov-comparison, dev = c("png", "pdf")}
pca_vi / pca_nri / pca_hr *
  theme(
    plot.title = element_text(size = 10),
    axis.title.y = element_text(size = 8),
    axis.title.x = element_text(size = 10)
  ) +
  patchwork::plot_annotation(
    tag_levels = list(c("VI", "NRI", "HR")),
    title = "Proportion of variance explained by principal components",
    subtitle = "Labels indicate cumulative % of explained variance"
  )
```

# Buffer extraction scenario sketches

Standard scenario (somewhere within the pixel)

```{r buffer-scenarios-1, echo=FALSE}
knitr::include_graphics("inst/tree-buffer-vis.png")
```

Extreme scenario (at the border of one or multiple pixels)

```{r buffer-scenarios-2, echo=FALSE}
knitr::include_graphics("inst/tree-buffer-vis-max.png")
```