analysis_congreso_forestal.Rmd

---
title: "Prepare analysis and figures for VIII Spanish Forest Meeting"
author: "ajpelu"
date: "2022-01-13"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
---

## Prepara Datos

```{r setup, echo=FALSE}
knitr::opts_chunk$set(echo = FALSE, 
                      warning = FALSE, 
                      message = FALSE,
                      fig.width=10, fig.height=7)
```

```{r pkg}
library(tidyverse)
library(here)
library(readxl)
library(ggstatsplot)
library(Metrics)
library(kableExtra)
library(ggtext)
library(ggpubr)
library(modelr)
library(rpart)
library(party)
library(ggparty)
library(partykit)
```

- Usamos datos de cobertura vegetal de las parcelas de campo (***cob.campo***) y datos derivados de dron (***cob.dron***). 

- De los datos de dron, utilizamos el método de estimación denominado COBERTURA (ver [análisis preliminar](compara_drone_campo.html))

- Los datos de campo corresponden al muestreo realizado el 19/05/2021.

- El vuelo del dron se realizó el día 21/05/2021. 

- Los ***rangos*** de cobertura se han reclasifiaco de acuerdo a:

| RANGO_INFOCA | Nombre           | Cobertura |
| :----------- | :--------------- | :-------- |
| 1            | "Matorral claro" | "<25%"    |
| 2            | "Matorral medio" | "25-50%"  |
| 3            | "Espartal denso" | ">75%"    |
| 4            | "Aulagar denso"  | ">75%"    |


```{r}
cob.raw <- read_excel(path=here::here("data/test_drone.xlsx"),
                 sheet = "COBERTURA") %>% 
  mutate(cob.campo = COB_TOTAL_M2*100, 
         cob.dron = COBERTURA*100)

diversidad <- read_excel(path=here::here("data/test_drone.xlsx"),
                 sheet = "SHANNON") %>% mutate(shannon = abs(I_SHANNON))

richness <- read_excel(path=here::here("data/riqueza_19_05_21.xlsx")) %>% 
  rename(QUADRAT = GEO_QUADRAT.NOMBRE) %>% 
  dplyr::select(QUADRAT, rich = RIQUEZA, rich_cor = RIQUEZA_COR)

slope <- read_csv(here::here("data/slopes_quadrat.csv")) %>% 
  rename(QUADRAT = NOMBRE, slope = Slope)

df <- cob.raw %>% inner_join(diversidad) %>% 
    mutate(coverclass = 
           case_when(
             RANGO_INFOCA == 1 ~ "Matorral claro (<25%)",
             RANGO_INFOCA == 2 ~ "Matorral medio (25-50%)",
             RANGO_INFOCA == 3 ~ "Espartal denso (>75%)",
             RANGO_INFOCA == 4 ~ "Matorral denso (>75%)"
           )) %>% 
  dplyr::select(QUADRAT, RANGO_INFOCA, coverclass, 
                cob.campo, cob.dron, shannon) %>% 
  mutate(coverclass = fct_relevel(coverclass, c("Matorral claro (<25%)","Matorral medio (25-50%)","Matorral denso (>75%)","Espartal denso (>75%)"))) %>% 
  inner_join(richness) %>% 
  inner_join(slope)
```

## Correlación General

```{r general-correlation-ggscat, dev=c("png","jpeg")}
ggscatterstats(df, 
               x= "cob.campo", 
               xlab = expression('Cobertura'['campo']*' (%)'), 
               y = "cob.dron", 
               ylab = expression('Cobertura'['dron']*' (%)'), 
               results.subtitle = TRUE,
               point.args = list(
                size = 3, 
                alpha = 0.4,
                colour = "blue"
               ),
               smooth.line.args = 
                 list(size = 1, color = "black"),
               xfill = "gray", yfill = "gray",
               marginal = TRUE, 
               ggplot.component = 
                 list(geom_abline(slope = 1, colour="blue"), 
                      xlim(0,100)))
```

```{r rmse-global}
df.rmse_global <- df %>% 
  summarise(rmse = round(
    Metrics::rmse(cob.dron, cob.campo),4),
    min = min(cob.campo), 
    max = max(cob.campo),
    rmsen.minmax = rmse / (max(cob.campo) - min(cob.campo))*100)
```

```{r}
# https://stackoverflow.com/questions/17022553/adding-r2-on-graph-with-facets
lm_eqn = function(df){
  m = lm(cob.dron ~ cob.campo, df);
  eq <- substitute(r2,
                   list(r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));                 
}
```

```{r general-correlation, dev=c("png","jpeg")}
g <- df %>% 
  ggplot(aes(x=cob.campo, y = cob.dron)) + 
  geom_point(size=3, alpha=.6, colour="blue") + 
  geom_abline(slope=1) +
  xlab(expression('Cobertura'['campo']*' (%)')) + 
  ylab(expression('Cobertura'['dron']*' (%)')) +
  xlim(0,100) + ylim(0,100) + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  ggtitle("Cobertura vegetal (%): Dron vs. Campo") + 
  theme(plot.title = element_text(size=12)) + 
  annotate("text", x= 20, y = 90, 
           label = paste0("R^2 == ", lm_eqn(df)), 
           parse = TRUE) +
  annotate("text", x= 20, y = 80, 
           label = paste0("RMSE = ", round(df.rmse_global$rmse, 2))) 
```

```{r}
ggsave(g,
       filename = here::here("output/general_correlacion.jpg"),
       device = "jpg", height = 9, width = 9, unit="cm", dpi = "print")
```


## Correlación por Rangos

- Explorar como varía la correlación en los diferentes rangos de  cobertura

- Computar el RMSE, y el RMSE normalizado. El RMSE es dependiente de la escala, por lo que necesitaríamos normalizar para poder comparar entre las clases de cobertura.

```{r}
df.rmse_groups <- df %>% group_by(coverclass) %>% 
  summarise(rmse = round(
    Metrics::rmse(cob.dron, cob.campo),4),
    min = min(cob.campo), 
    max = max(cob.campo),
    rmsen.minmax = rmse / (max(cob.campo) - min(cob.campo))*100)

# see also hydroGOF pkg for RMSE et al. 
```

```{r}
df.rmse_groups %>% 
  kbl(col.names = c("Rango de cobertura",
                               "RMSE",
                               "min",
                               "max",
                               "norm. RMSE %"),
                 digits = c(0,2,0,0,2)) %>% 
  kable_material()
```

- Generamos las ecuaciones para la gráfica

```{r}
eqns <- by(df, df$coverclass, lm_eqn)
df.label <- data.frame(eq = unclass(eqns), coverclass = names(eqns))
df.label$lab = paste(df.label$coverclass, "R^2 =", df.label$eq, sep=" ")


r2_labeller <- function(variable,value){
  return(df.label$lab)
}
```


```{r correlation-coverclass, dev=c("png","jpeg")}
g.rangos <- df %>% 
  ggplot(aes(x=cob.campo, y = cob.dron, color=as.factor(coverclass)))  + 
  geom_abline(slope=1) + 
  geom_point(size=3, alpha = .5) + 
  #  facet_wrap(~coverclass, labeller = r2_labeller) + 
  facet_wrap(~coverclass, labeller = label_value) + 
  theme_bw() + 
  ylab("Dron") + xlab("Campo") + 
  # xlab(expression('Cobertura'['campo']*' (%)')) + 
  # ylab(expression('Cobertura'['dron']*' (%)')) +
  xlim(0,100) + ylim(0,100) + 
  theme(
    legend.position = "none", 
    panel.grid = element_blank(), 
    strip.background = element_rect(fill="white"), 
    strip.text = element_text(face = "bold"), 
    axis.title = element_text(face = "bold")
  ) + 
  ggtitle("Cobertura vegetal (%)") +
  geom_richtext(data = df.rmse_groups, size=3,
            aes(x = 30, y = 90, 
                label = paste0(
                  "RMSE<sub>norm.</sub> = ",
                  round(rmsen.minmax,2), " %")),
            fill = NA, label.color = NA)

g.rangos
```

```{r}
ggsave(g.rangos,
       filename = here::here("output/correlacion_rangos.jpg"),
       device = "jpg", height = 13, width = 13, unit="cm", dpi = "print")
```


## Influencia de otras variables en la Variación de la correlación

¿Existe alguna relación entre la correlación y otras variables? Podría interesarnos explorar cómo otras variables podrían influir en la correlación dron-campo, *por ejemplo* la riqueza o la pendiente. Se pueden utilizar varios enfoques (análisis exploratorio, residuos, etc.). En nuestro caso utilizamos la correlación entre los residuos de la correlación y las diferentes variables.  

- Calculamos los residuos y los residuos absolutos 

```{r, echo = TRUE}
m <- lm(cob.dron ~ cob.campo, data=df)
df <- df %>% modelr::add_residuals(m) %>% 
  mutate(resid.abs = abs(resid))


dfres <- df %>% dplyr::select(coverclass, Diversidad = shannon, Riqueza  = rich, Pendiente = slope, resid, resid.abs) %>% 
  pivot_longer(cols = c("Diversidad", "Riqueza", "Pendiente")) %>% 
  mutate(variable = fct_relevel(name, c("Diversidad", "Riqueza", "Pendiente")))
  
```

- Hacemos gráfico de las tres variables 

```{r residuos-variables, dev=c("png","jpeg")}
p <- ggpubr::ggscatter(dfres,
  x = "value", y = "resid.abs",
  color = "black",
  alpha = 0.7,
  xlab = "", 
  ylab = expression(paste("|","Residuos","|")),
  add = "reg.line",
  add.params = list(color = "blue", fill = "lightgray"),
  conf.int = TRUE, 
  facet.by = "variable"
) + 
  stat_cor(
    label.y.npc="top", label.x.npc = "left",
    aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))
    )

p.resid <- ggpubr::facet(p,
  facet.by = "variable", scales = "free_x",
  panel.labs.background = list(fill = "white")
)
p.resid
```

```{r}
ggsave(p.resid,
       filename = here::here("output/residus_variables.jpg"),
       device = "jpg", height = 8, width = 18, unit="cm", dpi = "print")
```



```{r, echo=FALSE, eval=FALSE}
# alternativa
dfres %>% 
  ggplot(aes(x=value, y=resid.abs)) + 
  geom_point() + 
  facet_wrap(~variable, scales = "free_x") +
  geom_smooth(method = lm, fill = "lightgray") + 
  ggpubr::stat_cor(
    label.y.npc="top", label.x.npc = "left",
    aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))
    ) +
  ggpubr::

```



### Correlación general vs. Riqueza 

```{r}
p.riqueza <- df %>% 
  rename(Riqueza = rich) %>% 
  ggplot(aes(x=cob.campo, y = cob.dron)) + 
  geom_point(aes(
    size=Riqueza), alpha = .4) +  
  geom_abline(slope=1) +
  theme_bw() + 
  xlab(expression('Cobertura'['campo']*' (%)')) + 
  ylab(expression('Cobertura'['dron']*' (%)')) +
  theme(legend.position = "bottom") + 
  # scale_colour_gradient(low="#fee8c8", high = "red") + 
  ggtitle("Dron vs. Campo | Riqueza")
```

```{r}
p.riqueza
```



```{r}
ggsave(p.riqueza,
       filename = here::here("output/correlacion_size_riqueza.jpg"),
       device = "jpg", height = 10, width = 10, unit="cm", dpi = "print")
```


## Comparison dronQ vs Quadrat 

- ¿Es la relación entre dronQ-Quadrat uniforme a lo largo de todos los valores de cobertura o existen grupos dentro de esta correlación? Dicho de otro modo, ¿la correlación entre las dos variables es homogénea a lo largo de todos los valores de cobertura?

- Para ello aplicamos una técnica de partición de la varianza, y hacemos árboles de clasificación. 

### Primer modelo 

```{r}
ct <- partykit::ctree(cob.dron ~ cob.campo, data=df)
ct
sctest(ct) 
```

Nos indica que existen al menos 5 grupos. 

```{r}
ggparty(ct) + 
  geom_edge() +
  geom_edge_label(colour = "grey", size = 4) +
  geom_node_plot(gglist = list(
    geom_point(aes(x = cob.campo,y = cob.dron)),
    geom_smooth(aes(x=cob.campo, y=cob.dron), method = lm, se = FALSE),
    geom_abline(slope = 1, colour = "gray"),
    theme_bw(base_size = 10),
    xlab(expression('Cobertura'['campo']*' (%)')),
    ylab(expression('Cobertura'['dron']*' (%)')) 
    ),
                scales = "fixed", 
                shared_axis_labels = TRUE,
                shared_legend = TRUE,
                legend_separator = TRUE, 
                id = "terminal") +
   geom_node_label(aes(col = splitvar),
                  line_list = list(aes(label = paste("Node", id)),
                                   aes(label = splitvar),
                                   aes(label = scales::pvalue(p.value, accuracy = 0.001, 
                                                              decimal.mark = ".", 
                                                              add_p = TRUE))), 
                                                     
          line_gpar = list(list(size = 8, col = "black", fontface = "bold"),
                                   list(size = 8),
                                   list(size = 8)),
                  ids = "inner") + 
  geom_node_label(aes(label = paste0("Node ", id, ", (n= ", nodesize, ")")),
                  fontface = "bold",
                  ids = "terminal",
                  size = 2, 
                  nudge_y = 0.01) +
  theme(legend.position = "none")

```

pero, ¿no estaremos sobreparametrizando? 


```{r}
ggplot(df, aes(x=cob.campo, y = cob.dron)) +
  geom_point() + 
  geom_abline(yintercept = 1) + 
  geom_smooth(
    data = (df %>% filter(cob.campo <=  13)),
    aes(x=cob.campo, y = cob.dron),
    method = "lm") + 
  geom_smooth(
    data = (df %>% filter(cob.campo > 13 & cob.campo <= 23)),
    aes(x=cob.campo, y = cob.dron),
    method = "lm") + 
  geom_smooth(
    data = (df %>% filter(cob.campo > 23 & cob.campo <= 36)),
    aes(x=cob.campo, y = cob.dron),
    method = "lm") + 
  geom_smooth(
    data = (df %>% filter(cob.campo > 36 & cob.campo <= 50)),
    aes(x=cob.campo, y = cob.dron),
    method = "lm") +
  geom_smooth(
    data = (df %>% filter(cob.campo >50)),
    aes(x=cob.campo, y = cob.dron),
    method = "lm")

```


#### Exploprar el Overfitting (sobreparametrización)

- Hacer crecer un árbol sobre cierto nivel de complejidad, puede llevarnos a sobreparametrizar el modelo. ¿En que nodo cortamos el árbol? 

- Es preferible tener un modelo más simple que sobreparametrizar el ajuste. Para ello podemos utilizar el Complexity parameter (CP), que controla el tamaño del árbol de decisión. Si el coste de añadir otro nodo (u otra variable) a partir del nodo actual está por encima del valor del cp, entonces el árbol no continúa creciendo. Mas información [aquí](https://statinfer.com/203-3-10-pruning-a-decision-tree-in-r/)


```{r}
set.seed(123)
ctrpart <- rpart(cob.dron ~ cob.campo, data=df)
ctrpart
printcp(ctrpart)
```

Al explorar este modelo observamos que el error (xerror) no se disminuye sustancialmente a partir del 3 item (nsplit = 2). 

```{r}
plotcp(ctrpart)
```

Volvemos a modelar, usando un cp = 0.07. Obtenemos, que no reducimos mucho el error. 

```{r}
set.seed(123)
ctrpart2 <- rpart(cob.dron ~ cob.campo, data=df, control = 
               rpart.control(minsplit = 2, cp = .07))
ctrpart2
printcp(ctrpart2)
plotcp(ctrpart2)
```

Por tanto expecificamos un complexity parameter de 0.22, y comprobamos 
```{r}
set.seed(123)
ctrpart3 <- rpart(cob.dron ~ cob.campo, data=df, control = 
               rpart.control(minsplit = 2, cp = .22))
ctrpart3
printcp(ctrpart3)
plotcp(ctrpart3)
```

En este caso, observamos que el mejor modelo es aquel que usa dos grupos (1 split) 

```{r}
set.seed(123)
ct_ok <- partykit::ctree(cob.dron ~ cob.campo, data=df, 
                          control = 
                            ctree_control(minsplit= 1,alpha = 0.05, maxdepth = 1))
plot(ct_ok)

ggparty(ct_ok) + 
  geom_edge() +
  geom_edge_label(colour = "grey", size = 4) + 
  geom_node_plot(gglist = list(
    geom_point(aes(x = cob.campo,y = cob.dron)),
    geom_smooth(aes(x=cob.campo, y=cob.dron), method = lm, se = FALSE),
    geom_abline(slope = 1, colour = "gray"),
    theme_bw(base_size = 10),
    xlab(expression('Cobertura'['campo']*' (%)')),
    ylab(expression('Cobertura'['dron']*' (%)')) 
    ),
                scales = "fixed", 
                shared_axis_labels = TRUE,
                shared_legend = TRUE,
                legend_separator = TRUE, 
                id = "terminal") +
   geom_node_label(aes(col = splitvar),
                  line_list = list(aes(label = paste("Node", id)),
                                   aes(label = splitvar),
                                   aes(label = scales::pvalue(p.value, accuracy = 0.001, 
                                                              decimal.mark = ".", 
                                                              add_p = TRUE))), 
                                                     
          line_gpar = list(list(size = 8, col = "black", fontface = "bold"),
                                   list(size = 8),
                                   list(size = 8)),
                  ids = "inner") + 
  geom_node_label(aes(label = paste0("Node ", id, ", (n= ", nodesize, ")")),
                  fontface = "bold",
                  ids = "terminal",
                  size = 2, 
                  nudge_y = 0.01) +
  theme(legend.position = "none")
```


```{r}
dff <- df %>% 
  mutate(clase = case_when(
    cob.campo <= 36 ~ as.character("cob_low"),
    TRUE ~ "cob_high"
  ))

ggscatter(dff, 
          x = "cob.campo", 
          y = "cob.dron",
          add = "reg.line",
  conf.int = TRUE, 
  color = "clase", shape = "clase", 
  xlab = expression('Cobertura'['campo']*' (%)'), 
  ylab = expression('Cobertura'['dron']*' (%)')) + 
  xlim(0,100) + ylim(0,100) + 
  stat_cor(
    label.y.npc="top", label.x.npc = "left",
    aes(color = clase, 
        label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    p.accuracy = 0.001,
    r.accuracy = 0.001
    )  + 
  geom_abline(slope = 1) 

```

Y los RMSE
```{r}
dff.rmse_groups <- dff %>% group_by(clase) %>% 
  summarise(rmse = round(
    Metrics::rmse(cob.dron, cob.campo),4),
    min = min(cob.campo), 
    max = max(cob.campo),
    rmsen.minmax = rmse / (max(cob.campo) - min(cob.campo))*100)

dff.rmse_groups

```