4_figures_Monitoring-Paper.Rmd

---
title: "R Notebook"
output:
  html_document:
    toc: true
    toc_float: true
    theme: "yeti"
    df_print: paged
  pdf_document: default
---
This document presents the subset of the figures used for paper about monitoring.

## Load packages

```{r message=FALSE, warning=FALSE}
library(MASS)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)
library(rgdal)
library(multcompView)
library(car)
library(ggpmisc)
library(chisq.posthoc.test)
library(vcd)
library(grid)
library(ggpattern)
library(car)
library(gridExtra)
library(cowplot)
```

Sampling data correspond to the data collected with kobo and previously cleaned with the script `1_preprocesamiento_datos_kobo.Rmd`.

## Load data

```{r, message=FALSE}
# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)

#Add column block to muestreo_tidy
muestreo_tidy<-muestreo_tidy%>% mutate(block = case_when(plot %in% c(1,2,3) ~1,
                 plot %in% c(1,2,3) ~1,
                 plot %in% c(4,5,6) ~2,
                 plot %in% c(7,8,9) ~3,
                 plot %in% c(10,11,12) ~4,
                 plot %in% c(13,14,15) ~5,
                 plot %in% c(16,17,18) ~6,
                 plot %in% c(19,20,21) ~7,
                 plot %in% c(22,23,24) ~8,
                 plot %in% c(25,26,27) ~9,
                 plot %in% c(28,29,30) ~10,
                 plot %in% c(31,32,33) ~11,
                 plot %in% c(34,35,36) ~12,
                 plot %in% c(37,38,39) ~13,
                 plot %in% c(40,41,42) ~14,
                 plot %in% c(43,44,45) ~15,
                 plot %in% c(46,47,48) ~16)) %>%
                mutate(block=as.factor(block)) # store as factor, instead of numeric (numbers are IDs)

parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)

# pivot long plots data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy, 
                            cols = healthy:worm, 
                            names_to = "tree_health_simplified",
                            values_to = "n_trees")
```

Data analyzed here correspond only to the trees that were approved during the validation by manually reviewing the photographs in kobotoolbox. Total of `r nrow (muestreo_tidy)` trees sampled, `r sum (muestreo_tidy$X_validation_status=="validation_status_approved")` were approved in the validation.

```{r}
muestreo_tidy<- filter(muestreo_tidy, X_validation_status=="validation_status_approved")

# write.csv(muestreo_tidy,  file="../../../../Desktop/muestreo_tidy.csv")
```

Color palettes:

```{r, message=FALSE}
# Make a nice color pallete and legend order for all plots

my_cols=c("darkgreen", 
              "darkred", 
              "orangered1", 
              "cadetblue", 
              "tan", 
              "beige", 
            #  "burlywood4", 
              "coral", 
              "aquamarine3", 
              "gray70", 
              "black")

desired_order=c("healthy", 
                "ozone", 
                "ozone_and_other", 
                "others_combined", 
                "drougth", 
                "fungi", 
             #   "insect", 
                "worm", 
                "acid_rain", 
                "other", 
                "dead")

desired_names=c("healthy", 
                "ozone", 
                "ozone and other", 
                "others combined", 
                "drougth", 
                "fungi", 
             #   "insect", 
                "worm", 
                "acid rain", 
                "other", 
                "dead")

spanish_labels=c("Sano", 
                  "Ozono", 
                  "Ozono y otros", 
                  "Otros combinados no-ozono", 
                  "Sequía", 
                  "Hongos", 
               #   "Insectos", 
                  "Gusano de seda", 
                  "Lluvia acida", 
                  "Otro", 
                  "Muerto")

# For ozone damage percentage 
 my_cols2<-c("darkgreen", "gold2", "chocolate1", "orangered", "red4", "darkorchid4")
 
desired_order_percentage<-c("0%","less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")


```

Multiplot fun:

```{r}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```

Configure google api for maps:

``` {r}
# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-g0ogle-map.html

## configure google api

# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).

# if you get the following error when running get_map():

#"Error in aperm.default(map, c(2, 1, 3)) : 
#  invalid first argument, must be an array " 

# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

##  load and register api
api <- readLines("api_key.api")
register_google(key = api)
```

Map and monitoring figures presented in the paper:

# Figure 2

## Plot 2a: PNDL location on CDMX map

```{r, fig.height=6, fig.width=6, message=FALSE}
# get cdmx shape
CDMX<-readOGR(dsn="../data/spatial", layer="CDMX")
CDMX<-fortify(CDMX)

# get PNDL shape
PNDL<-readOGR(dsn="../data/spatial", layer="Desierto_Leones_Geo_ITRF08")
PNDL<-fortify(PNDL)

# get background map
sat_map = get_map(location = c(lon = -99.133549, lat = 19.3), zoom = 10, maptype = 'terrain-background', source = "google")

## plot
p_a<-ggmap(sat_map) + 
            geom_polygon(data = CDMX,
                         aes(x = long, y = lat, group = group),
                         color="black", fill=NA, size=1.5) +
            geom_polygon(data = PNDL,
                         aes(x = long, y = lat, group = group),
                         color="red", fill=NA, size=1.5) +
            geom_point(aes(x=-98.95, y=19.6), 
                       shape=0, stroke=2, size=5, color="black") +
            geom_point(aes(x=-98.95, y=19.55), 
                       shape=0, stroke=2, size=5, color="red") +
            geom_text(aes(label="CDMX", x=-98.87, y=19.6), 
                      color="Black", fontface="bold", size=5) +
            geom_text(aes(label="PNDL", x=-98.87, y=19.55), 
                      color="Black", fontface="bold", size=5) +
            theme(text = element_text(size = 20))+
  ggtitle("a)")
```

## Plot 2b: Satellite image and surroundings of the PNDL

```{r, fig.height=6, fig.width=6, message=FALSE}
# get background map
sat_map = get_map(location = c(lon = -99.30, lat = 19.31), zoom = 13, maptype = 'satellite', source = "google")

## add towns names
towns<-data.frame(nombre=c("San Bartolo Ameyalco", 
                           "Santa Rosa Xochiac", 
                           "San Mateo Tlaltenango"),
                  long=c(-99.270, -99.29, -99.276),
                  lat=c(19.333, 19.325, 19.346))



## plot
p_b<-ggmap(sat_map) + 
            geom_polygon(data = PNDL,
                         aes(x = long, y = lat, group = group),
                         color="red", fill=NA, size=1.5) +
            geom_point(data=towns, aes(x=long, y=lat), colour="red", size=1.5) +
            geom_text(data=towns, aes(label=nombre, x=long, y=lat), 
                      color="white", fontface="bold",
                      size=5, nudge_y=0.003) +
  # add Cruz de Coloxtitla (CX), and Convento (Cn) landmarks
            geom_text(aes(label="X", x=-99.3014, y=19.286068), 
                      color="white", fontface="bold", size=4) +
            geom_text(aes(label="C", x=-99.31, y=19.3133), 
                      color="white", fontface="bold", size=4) +
            theme(text = element_text(size = 20))+
  ggtitle("b)")
```

## Plot 2c: This is the distribution of the 48 plots:

```{r, message=FALSE, fig.height=15, fig.width=15}

## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")

# plot sampled plots
p_c <-  ggmap(sat_map)
p_c <- p_c + geom_point(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude),
                      color="red") +
          geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                     check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5) +
    theme(text = element_text(size = 20))+
  ggtitle("c)")
p_c 

```


## Plot 2d: Distribution of tree health status by plot

The following figure shows the total number of trees sampled in each 10x10 m plot, and how many of these are under some category of damage:

```{r, message=FALSE, fig.height=9, fig.width=11}

parcelas_HOO<-parcelas_long%>%
  filter(tree_health_simplified ==  "ozone" )
  sum(parcelas_HOO$n_trees)

parcelas_HOO<-parcelas_long%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" | tree_health_simplified == "ozone_and_other" )


p_d <- ggplot(parcelas_HOO, aes(x=plot, y=n_trees,     fill=tree_health_simplified)) +
  geom_bar(stat="identity") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= desired_names,
                    name= "Health status") 
  

p_d <- p_d + theme_bw() +
  labs(x="Plots", y= "Number of trees") +
  theme(text = element_text(size = 20)) +
  ggtitle("d)")+
  coord_flip()

```

# Multiplot Fig 2

```{r, message=FALSE, fig.height=12, fig.width=16}

fig_2<-multiplot(p_a, p_c, p_b, p_d, cols=2)
fig_2
ggsave("../outputs/Fig_2.png",
       dpi = 600)
```

# ANOVAs related to  Figure 3

One-way ANOVAS and 2-way ANOVA
```{r}

#### subset of data only with reforestation, health and exposition
refo_expo<-muestreo_tidy %>% 
              dplyr::select(plot, block, tree_health_simplified, reforested, tree_exposition) %>% 
           # discard other typos of damage
             filter(tree_health_simplified == "healthy" |
                    tree_health_simplified == "ozone" |
                    tree_health_simplified == "ozone_and_other") %>% 
          # make "ozone" and "ozone and other" a single category = damaged
             mutate(tree_health_simplified = replace(tree_health_simplified, tree_health_simplified == "ozone", "damaged")) %>% 
             mutate(tree_health_simplified = replace(tree_health_simplified, tree_health_simplified == "ozone_and_other", "damaged"))
      


### ONE-way ANOVA
library(car)

#### By reforestation
summary_by_ref<- refo_expo %>%  
  group_by(plot, block, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

#now we have n for each category
head(summary_by_ref)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_ref<-transform(summary_by_ref, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_ref) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_ref<-summary_by_ref %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore

# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK





## By exposition
summary_by_exp<- refo_expo %>%  
  group_by(plot, block, tree_exposition,  tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

#now we have n for each category
head(summary_by_exp)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_exp<-transform(summary_by_exp, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_exp) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_exp<-summary_by_exp %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore



#### One way ANOVAS

# for reforestation
Anova(lm(summary_by_ref$perc.damage.log ~ summary_by_ref$reforested + summary_by_ref$block))


# for exposition
Anova(lm(summary_by_exp$perc.damage.log ~ summary_by_exp$tree_exposition + summary_by_exp$block))



#### For TWO-way ANOVA

### Summarise and estimate n and % of DAMAGED trees (normalized) by category and plot

summary_by_both<- refo_expo %>%  
  group_by(plot, block, tree_exposition, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

#now we have n for each category
head(summary_by_both)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_both<-transform(summary_by_both, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_both) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_both<-summary_by_both %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore


# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK

head(summary_by_both) 

##### TWO WAY ANOVA By BOTH exposition and tree origin 

### get data in correct format
summary_by_both<- refo_expo %>%  
  group_by(plot, block, tree_exposition, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

#now we have n for each category
head(summary_by_both)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_both<-transform(summary_by_both, 
                           perc.damage = ave(n,
                                             plot,
                                             FUN = prop.table))  %>% 
  
  # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_both) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_both<-summary_by_both %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore


# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK

head(summary_by_both) 

```

# Figure S10

```{r}
### TWO-way ANOVA
Anova(lm(summary_by_both$perc.damage.log ~ summary_by_both$tree_exposition * summary_by_both$reforested + summary_by_both$block))

```


## Figure 3 Boxplot of data used for ANOVA

```{r, fig.height=4, fig.width=4, message=FALSE}

# Boxplot tree expostion and reforestation origin on percentage of damage
p3<- ggplot(summary_by_both, aes(x = factor(tree_exposition), 
                                 y = perc.damage.log, color = factor(reforested))) +
  geom_boxplot() +
  ylab("Log-transformed percentage of damage by plot") +
  scale_color_manual(name= "Tree origin",
                     breaks = c("no", "yes"),
                     labels = c("Natural", "Reforested"),
                     values=c("#9f964b", "#9a5ea1")) + # nice colour blind colors
  scale_x_discrete(name="Tree exposition",
                   breaks = c("cover", "exposed"),labels = c("Nursed", "Exposed")) + 
  theme_bw()+ 
  theme(text = element_text(size = 20))

p3

```

# Figure 4

## Plot 4a

```{r, fig.height=6, fig.width=6, message=FALSE}


condition_HOO<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |
           tree_health_simplified == "ozone_and_other" )


p <- filter(condition_HOO, tree_heigth<15, tree_nodes>0) %>% 
     ggplot(.) +
     scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= desired_names,
                    name= "Health status") +
theme_bw()

p4_a <- p + geom_histogram(aes(x=tree_nodes, fill=tree_health_simplified))  +
    labs(x="Tree age (years)", y= "Number of trees") +
    theme(text = element_text(size = 20)) +
    theme(plot.title = element_text(lineheight=1.1))+
  ggtitle("a)")
p4_a
```

## Plot 4c

```{r, fig.height=6, fig.width=6, message=FALSE}
## base data
# Definir plantas sanas y dañadas por otra cosa que no fuera ozono
# cond_PO<- se  refiere a condition Percentage damage by Ozone 
cond_PO<-as_data_frame(muestreo_tidy)
# Asignar 0% de daño por ozono a los árboles healthy
cond_PO$ozone_damage_percentage = ifelse(cond_PO$tree_health == "healthy", "0%", cond_PO$ozone_damage_percentage)

# Filtrar por porcentaje de daño
condition_PO<-cond_PO%>%
  filter(ozone_damage_percentage == "0%" | ozone_damage_percentage == "less than 10%" | ozone_damage_percentage == "10 to 40%" | ozone_damage_percentage == "40 to 50%"| ozone_damage_percentage == "50 to 70%" | ozone_damage_percentage == "more than 70%")

condition_PO$ozone_damage_percentage <- as.factor(condition_PO$ozone_damage_percentage)

library(sf)
library(ggpattern)
library(magrittr)
library(ggplot2)
# Plot
p_od<- condition_PO %>% filter(!is.na(ozone_damage_percentage)) %>%
  ggplot(., aes(x = tree_nodes, fill = ozone_damage_percentage)) +
  geom_bar_pattern(aes(pattern = ozone_damage_percentage),
                   color = "black",
                   pattern_fill = "black",
                   pattern_density = 0.1,
                   pattern_spacing = 0.025,
                   pattern_key_scale_factor = 0.6) +
  scale_fill_manual(values= my_cols2, 
                    breaks = desired_order_percentage,
                    labels = c("0%","less 10%", "10 to 40%", "40 to 50%",
                               "50 to 70%", "more 70%"),
                    name= "Ozone damage\n per tree") +
  scale_pattern_manual(name= "Ozone damage\n per tree",
                       values =c("none","none", "stripe", "none",
                                 "none", "none"),
                       labels = c("0%","less 10%", "10 to 40%", "40 to 50%",
                                  "50 to 70%", "more 70%"))+
  theme_bw()+ theme(text = element_text(size = 20)) 


p_od
p4_c <- p_od +
  labs(x="Tree age (years)", y= "Number of trees") +
  theme(legend.title.align = 0.5)+
  theme(text = element_text(size = 20)) +
  theme(plot.title = element_text(lineheight=1.1))+
  ggtitle("c)")

p4_c

```

## Plot 4b

```{r, fig.height=6, fig.width=6, message=FALSE}

# Filtrar por categoría de daño
condition_HOO<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" | tree_health_simplified == "ozone_and_other" )

condition_HOO$tree_health_simplified <- as.factor(condition_HOO$tree_health_simplified)



# Data distribution

# Los datos tienen a graficar es el número de nodos para cada categoria de salud.
# Los datos son continuos discretos, por lo tanto el analisis a seguir para buscar diferencias entre los grupos son:

shapiro.test(condition_HOO$tree_nodes) #NO HAY NORMALIDAD 

# Procedo a hacer una prueba no paramétrica (kruskal)

# Debe tener Homocedasticidad: la varianza debe de ser constante entre todos los grupos.
# If the p-value is less than our chosen significance level, we can reject the null hypothesis and conclude that we have enough evidence to state that the variance among the groups is not equal.
leveneTest(sqrt(tree_nodes) ~ tree_health_simplified, data = condition_HOO) #Opera con medias

# As the p value obtained from the Shapiro-Wilk test and Levene’s test is significant (p < 0.05), we conclude that the data is not normally distributed and does not have equal variance. Kruskal-Wallis test is more appropriate for analyzing differences.

kruskal.test(sqrt(tree_nodes) ~ tree_health_simplified, data = condition_HOO) #opera con mediana

#As the p value obtained from the Kruskal-Wallis test test is significant , we conclude that there are significant differences.

# For the Kruskal-Wallis test, epsilon-squared is a method of choice for effect size measurement. The epsilon-squared is 0.74 and suggests a very strong effect of plant varieties on yield

# calculate effect size. Hay efecto relativamente fuerte
library(rcompanion)
epsilonSquared(x = sqrt(condition_HOO$tree_nodes), g = condition_HOO$tree_health_simplified)

# https://peterstatistics.com/CrashCourse/3-TwoVarUnpair/NomOrd/NomOrd3c.html
# To know which plant varieties are significantly different from each other, we will perform the Dunn’s test as post-hoc test for significant Kruskal-Wallis test. As there are multiple comparisons, we will correct the p values using Benjamini-Hochberg FDR method for multiple hypothesis testing at a 5% cut-off. The other tests that can be used for post-hoc test includes Conover and Nemenyi tests. Dunn’s test is more appropriate post-hoc than the Mann-Whitney U test for significant Kruskal-Wallis test as it retains the rank sums of the Kruskal-Wallis.

#poshoc para saber que grupos difieren
pairwise.wilcox.test(x = sqrt(condition_HOO$tree_nodes), g = condition_HOO$tree_health_simplified, p.adjust.method = "holm" )


# Perorm pairwise comparisons
library(ggpubr)
my_comparisons <- list( c("healthy", "ozone"), c("healthy", "ozone_and_other"), c("ozone", "ozone_and_other") )


p4_b<-condition_HOO%>%
   ggplot(aes(y= tree_health_simplified, x= tree_nodes))+
        scale_color_manual(values=  my_cols, labels= desired_names,
                    name= "")+
        geom_point(position="jitter",aes(color = tree_health_simplified), alpha=0.5, size= 2.5)+
  geom_boxplot(color="black", notch = F, alpha = 0.1)+
        xlab("Tree age (years)")+ ylab("Health status")+
  theme_bw()+
  ggtitle("b)")+
  theme(text = element_text(size = 20), axis.text.y=element_blank())+
  theme(plot.title = element_text(lineheight=1.1))+
  scale_x_continuous(breaks=c(5,10,15, 20, 25))+
  stat_compare_means(comparisons = my_comparisons, label = "p.signif")
p4_b


```

## Plot 4d

```{r, fig.height=6, fig.width=6, message=FALSE}
# Perorm pairwise comparisons
# Run Shapiro-Wilk test
shapiro.test(condition_PO$tree_nodes) #NO HAY NORMALIDAD

# Procedo a hacer un kruskal
# Debe tener homogeneidad. Si es mayor a 0.05 No hay evidencias en contra de la homogeneidad de varianzas. 
leveneTest(sqrt(tree_nodes) ~ ozone_damage_percentage, data = condition_PO, center = "median")

kruskal.test(sqrt(tree_nodes) ~ ozone_damage_percentage, data = condition_PO)

epsilonSquared(x = sqrt(condition_PO$tree_nodes), g = condition_PO$ozone_damage_percentage)


#poshoc que grupos difieren

#order levels 
condition_PO$ozone_damage_percentage<- ordered(condition_PO$ozone_damage_percentage, levels=c("0%","less than 10%", "10 to 40%", "40 to 50%", "50 to 70%","more than 70%"))

#Prueba
pairwise.wilcox.test(x = sqrt(condition_PO$tree_nodes), g = condition_PO$ozone_damage_percentage, p.adjust.method = "bonferroni" )

my_comparisons <- list(c("0%", "less than 10%"), c("0%", "10 to 40%"), c("0%", "40 to 50%"),
                        c("0%", "50 to 70%"), c("0%", "more than 70%"),
                        c("less than 10%", "50 to 70%"), 
                       c("10 to 40%", "50 to 70%"))

# Plot
p4_d<-condition_PO%>% filter(!is.na(ozone_damage_percentage)) %>% 
   ggplot(aes(y=ozone_damage_percentage , x= tree_nodes))+
        scale_color_manual(values=  my_cols2,labels = c("0%","less 10%", "10 to 40%", "40 to 50%",
                                         "50 to 70%", "more 70%"))+
        geom_point(position="jitter",aes(color = ozone_damage_percentage), alpha=0.5, size= 2.5)+
        xlab("Tree age (years)")+ ylab("Ozone damage per tree")+
  labs(color = "")+
          geom_boxplot(color="black", notch = F, alpha = 0.1)+
  theme_bw()+
  ggtitle("d)")+
  theme(legend.title.align = 0.5)+
  theme(text = element_text(size = 20),axis.text.y=element_blank())+
  theme(plot.title = element_text(lineheight=1.1))+
  scale_x_continuous(breaks=c(5,10,15, 20, 25))+
  ggpubr::stat_compare_means(comparisons = my_comparisons, label="p.signif")


p4_d
```

# Multiplot Fig 4

```{r, message=FALSE, fig.height=12, fig.width=16}
multiplot(p4_a, p4_c, p4_b, p4_d, cols=2)
```

# Figure 5

## Plot model all  data WITHOUT cohort

```{r, height=6, fig.width=10, message=FALSE}
condition_model_all<-dplyr::select(muestreo_tidy, contains(c("tree_heigth","tree_nodes","tree_health_simplified", "reforested", "tree_exposition", "plot")))%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |tree_health_simplified == "ozone_and_other" )%>%
  filter(tree_nodes < 15)

# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5all<-g0+geom_point(aes(colour=tree_health_simplified, 
                           shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified), fullrange =T)+
  scale_color_manual(values=  my_cols, labels= desired_names,
                     name= "Health status")+ 
  scale_shape_manual(name= "Health status", labels=desired_names, values = c(15,16,17)) + 
  labs(x="Tree age (years)", y= "Tree height")+
  theme_bw()+ 
  ggtitle("a)")
fig5all

```

## Select cohort data

```{r, fig.height=6, fig.width=6, message=FALSE}
cohort_data<- dplyr::select(condition_HOO, contains(c("tree_health_simplified", "tree_nodes")))%>% filter(tree_health_simplified == "ozone")
summary(cohort_data)
# We choose 7 & 11
```

## Model all data with cohort 7-11

```{r, height=6, fig.width=10, message=FALSE}

condition_model_all_C<-dplyr::select(muestreo_tidy, contains(c("tree_heigth","tree_nodes","tree_health_simplified", "reforested", "tree_exposition", "plot")))%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |tree_health_simplified == "ozone_and_other" )%>%
  filter(tree_nodes %in% (7:11))%>%
  filter(tree_nodes < 15)


########################
library(lme4)

########################

# Model all independent variables WITH health status interaction 
glm_1_coh<-glm(tree_heigth ~ tree_nodes, data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_1_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_1_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested data
glm_2_coh<-glmer(tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_2_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_2_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT tree exposition data
glm_3_coh<-glmer(tree_heigth ~ tree_nodes + tree_health_simplified + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_3_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_3_coh)


########################
# Model all independent variables WITH health status interaction WITHOUT reforested and tree exposition data
glm_4_coh<-glmer(tree_heigth ~ tree_nodes+tree_health_simplified + tree_exposition + (1|plot) , data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_4_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_4_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested and tree exposition data
glm_5_coh<-glmer(tree_heigth ~ tree_nodes + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_5_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_5_coh)


########################
glm_6_coh<-glmer(tree_heigth ~ tree_nodes*tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_6_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_6_coh)

########################
glm_7_coh<-glmer(tree_heigth ~ tree_nodes*tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_7_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_7_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested data
glm_8_coh<-glm(tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition + reforested, data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_8_coh)
#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_8_coh)


# Choose better model
AIC(glm_1_coh,glm_2_coh,glm_3_coh, glm_4_coh, glm_5_coh, glm_6_coh, glm_7_coh, glm_8_coh)


######
library(sjmisc)
library(lmerTest)
library(bruceR)
summary(glm_2_coh)
r2(glm_2_coh)
model_performance(glm_2_coh)

## Report glmer
```

# Figure S11
```{r, height=6, fig.width=10, message=FALSE}


print_table(glm_2_coh)
model_summary(glm_2_coh)

```

## Fig 5a: Plot model aLL data  with cohort


```{r, height=6, fig.width=10, message=FALSE}
#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all_C,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5a<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "inverse")), aes(color= tree_health_simplified, fullrange =T))+
  scale_color_manual(values=  my_cols, labels= desired_names,
                     name= "Health status")+ 
  scale_shape_manual(name= "Health status", labels=desired_names, values = c(15,16,17)) + 
  labs(x="Tree age (years)", y= "Tree height")+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("a)")
fig5a

```

## Fig 5b: Plot model ORIGIN data  with cohort

```{r, height=6, fig.width=10, message=FALSE}
#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

fig5b<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified, linetype = reforested), fullrange =T)+
  scale_color_manual(values=  my_cols, labels= desired_names,
                     name= "Health status")+ 
  scale_shape_manual(name= "Health status", labels=desired_names, values = c(15,16,17)) + 
  scale_linetype_manual(values= c("solid", "dashed"),labels = c("Natural \n regeneration", "Reforested"), name= "Origin")+
  labs(x="Tree age (years)", y= "Tree height")+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("b)")
fig5b


```

## Fig 5c: model EXPOSITION data  with cohort

```{r, height=6, fig.width=10, message=FALSE}
#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all_C,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5c<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
    ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified,linetype = tree_exposition), fullrange =T)+
  scale_color_manual(values=  my_cols, labels= desired_names,
                    name= "Health status")+
  scale_shape_manual(name= "Health status", labels=desired_names, values = c(15,16,17)) + 
  scale_linetype_manual(values= c("solid", "dashed"),labels = c("Nursed", "Exposed"), name= "Tree exposition",) +
  labs(x="Tree age (years)", y= "Tree height")+
  #stat_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), se = TRUE)+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("c)")
fig5c

```

# Multiplot Fig 5
```{r, fig.height=16, fig.width=10, message=FALSE}

multiplot(fig5a, fig5b, fig5c, cols=1)

```

# Figure S8
```{r, message=FALSE, fig.height=9, fig.width=11}
p_d <- ggplot(parcelas_long, aes(x=plot, y=n_trees, fill=tree_health_simplified)) +
  geom_bar_pattern(stat="identity",
                   aes(pattern = tree_health_simplified),
                   color = "black",
                   pattern_fill = "black",
                   pattern_density = 0.1,
                   pattern_spacing = 0.015,
                   pattern_key_scale_factor = 0.6) +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= desired_names,
                    name= "Health status") + 
  scale_pattern_manual(name= "Health status", breaks = desired_order,
                       values=c ("none",  "none", "none", "none", 
                                 "crosshatch", "none", "stripe", "circle", "none", "none"),
                       labels = desired_names)
  

p_d <- p_d + theme_bw() +
  labs(x="Plots", y= "Number of trees") +
  theme(text = element_text(size = 20)) +
  ggtitle("")+
  coord_flip()
p_d 

```

# Figure S12

```{r, height=20, fig.width=20, message=FALSE}
# plot pies in map
p_satmap <-  ggmap(sat_map)
p_satmap +geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = 1.5,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= desired_names,
                    name= "Health status") +
  theme(text = element_text(size = 20))

```

# Figure S13

```{r, height=6, fig.width=6, message=FALSE}
# Create new variable with percentage of ozone damage
parcelas_tidy<-parcelas_tidy %>% rowwise() %>% 
                     mutate(., 
                      total=sum(healthy,ozone,ozone_and_other,
                          drougth, acid_rain, other,
                          others_combined, dead, fungi,
                          # insect, 
                          worm)) %>%
                    mutate(perc.ozone= sum(ozone, ozone_and_other)/total)

#Statistical
leveneTest( parcelas_tidy$X_coordinates_altitude, parcelas_tidy$perc.ozone)


#plot
p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_altitude,
             y=perc.ozone))


p<- ggplot(parcelas_tidy, aes(X_coordinates_altitude, perc.ozone ))+
  geom_point(color= "grey50", size = 3, alpha = 0.6)

p<-p + 
  stat_smooth(color = "skyblue", formula = y ~ x,fill = "skyblue", method = "lm") +
  stat_poly_eq(
    aes(label = paste(..eq.label.., ..adj.rr.label.., sep = '~~~~')),
    formula = y ~ x,  parse = TRUE,
      size = 5, # Tamaño de fuente de la fórmula
             label.x = 0.1, #location, la proporción entre 0-1
      label.y = 0.95)+
  labs(x="Plot altitude", y= "Percentage of ozone damaged trees within a plot")+
  theme_bw() +
   theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 20), 
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        axis.text=element_text(size=9))

p
ggsave("../outputs/FigS3.png",
       dpi = 600)

```

# RESULTS Mean and sd of trees collected by each ranger during Participatory monitoring

```{r, height=6, fig.width=6, message=FALSE}
# No se incluye a M (37) y PP (6), ya que eran colectados por diferentes miembros del team
trees<- c(135,133,130,174,150,139,165,166,133,137,125,105)
rangers<- c("A","B","C","D","E","F","G","H", "I", "J", "K", "L")
mean(trees)
sd(trees)

```

# Percentage damage in plots

```{r, height=6, fig.width=6, message=FALSE}
#Percentage damage plots

PDP<-data.frame(parcelas_tidy$plot, parcelas_tidy$perc.ozone)

PDP<-PDP[order(PDP$parcelas_tidy.plot),]
range(PDP$parcelas_tidy.perc.ozone)
head(PDP[order(PDP$parcelas_tidy.perc.ozone),])
```

# Percentage per each health status

```{r, height=6, fig.width=6, message=FALSE}

# Create new variable with porcentage of ozone damage

status_HS<- c("healthy","ozone","ozone_and_other",
                          "drougth", "acid_rain", "other",
                          "others_combined", "dead", "fungi",
                          "worm")
count_HS<-c(sum(parcelas_tidy$healthy), sum(parcelas_tidy$ozone),
            sum(parcelas_tidy$ozone_and_other),sum(parcelas_tidy$drougth),
            sum(parcelas_tidy$acid_rain),sum(parcelas_tidy$other),
            sum(parcelas_tidy$others_combined), sum(parcelas_tidy$dead),
            sum(parcelas_tidy$fungi), sum(parcelas_tidy$worm))
sum(count_HS)

# Tree number to each health status
percentageH<-data.frame(status_HS, count_HS)
percentageH

# Percentage of ozone and ozone + other trees
(sum(c(parcelas_tidy$ozone, parcelas_tidy$ozone_and_other))*100)/sum(count_HS)

# Percentage of healthy trees
(sum(parcelas_tidy$healthy)*100)/sum(count_HS)

# Percentage of others damages in trees with dead trees
(sum(c(parcelas_tidy$drougth, parcelas_tidy$acid_rain , parcelas_tidy$other, parcelas_tidy$others_combined,parcelas_tidy$dead, parcelas_tidy$fungi, parcelas_tidy$worm))*100)/sum(count_HS)
```

```{r height=6, fig.width=6, message=FALSE}

# Filtrar por categoría de daño
condition_HOz<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone")

condition_HOz$tree_health_simplified <- as.factor(condition_HOz$tree_health_simplified)


# Damages by other NO ozone NO DEAD
levels(as.factor(muestreo_tidy$tree_health_simplified))
no_ozone<- dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "acid_rain"| tree_health_simplified == "fungi" | tree_health_simplified == "others_combined" | tree_health_simplified == "worm" | tree_health_simplified == "drougth"| tree_health_simplified == "other")

(609*100)/1765

# Percentage of DEAD trees
levels(as.factor(muestreo_tidy$tree_health_simplified))
dead<- dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "dead")
count(dead)
(43*100)/1765

sum(489,125, 499,1,9,70, 3,271,255, 43)
```
# Chi test

```{r, fig.height=6, fig.width=6, message=FALSE}

count(dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "ozone"))

count(dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "ozone_and_other"))

OBS_ozo<-c(125,499)
ESP_ozo<-c(0.5, 0.5)
chisq.test(x = OBS_ozo,
           p = ESP_ozo)

#224.1603
#p-value < 2.2e-16
# p<0.001

```