analysis_zona_time_postFire.Rmd

---
title: "Browsings effects on Soil propierties after fire (no control)"
author: "ajpelu"
date: "2021-09-08"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
library(knitr)
library(rmdformats)

## Global options
options(max.print="100")
opts_chunk$set(echo=FALSE,
	             cache=FALSE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE,
               autodep = TRUE)
```
# Introduction
- Analysis of temporal evolution of soil parameters along time. 
- Only for Autumn treatment (*i.e.* zona == "P"; zona == "NP") 
- Interpret zona as "grazing effect": 

  - zona == "P" corresponds to *Browsing*
  - zona == "NP" corresponds to *No Browsing*


```{r pkgs}
library("tidyverse")
library("here")
library("janitor")
library("ggpubr")
library("fitdistrplus")
library("lme4")
library("lmerTest")
library("emmeans")
library("multcomp")
library("performance")
library("afex")
library("DHARMa")
library("glmmTMB")
library("glmmADMB")
library("car")
library("kableExtra")
```

# Prepare data 
```{r prepare-data, echo=TRUE}
raw_soil <- readxl::read_excel(
  here::here("data/Resultados_Suelos_2018_2021_v2.xlsx"),
  sheet = "SEGUIMIENTO_SUELOS_sin_ouliers"
) %>% 
  janitor::clean_names() %>% 
  mutate(
    treatment_name =
      case_when(
        str_detect(geo_parcela_nombre, "NP_") ~ "Autumn Burning / No Browsing",
        str_detect(geo_parcela_nombre, "PR_") ~ "Spring Burning / Browsing",
        str_detect(geo_parcela_nombre, "P_") ~ "Autumn Burning / Browsing"
      ),
    zona =
      case_when(
        str_detect(geo_parcela_nombre, "NP_") ~ "QOt_NP",
        str_detect(geo_parcela_nombre, "PR_") ~ "QPr_P",
        str_detect(geo_parcela_nombre, "P_") ~ "QOt_P"
      ),
    fecha = lubridate::ymd(fecha),
    pre_post_quema = case_when(
      pre_post_quema == "Prequema" ~ "0 preQuema",
      pre_post_quema == "Postquema" ~ "1 postQuema"
    )
  ) 

```

- Compute date as months after fire 

```{r, echo=TRUE}
autumn_fire <- lubridate::ymd("2018-12-18")

soil <- raw_soil %>% 
  filter(zona != "QPr_P") %>% 
  filter(fecha != "2018-12-11") %>% 
  mutate(zona = as.factor(zona)) %>% 
  mutate(meses = as.factor(as.character(lubridate::interval(autumn_fire, lubridate::ymd(fecha)) %/% months(1))
  )) %>% 
  mutate(pastoreo = as.factor(
    case_when(
    zona == "QOt_P" ~ "Browsing",
    zona == "QOt_NP" ~ "No Browsing"
  ))) %>% 
  relocate(pastoreo, fecha, meses) %>% 
  dplyr::select(-pre_post_quema, -tratamiento)

xtabs(~meses+pastoreo, data = soil)


# sss <- soil %>% dplyr::select(meses, pastoreo, ca_percent)
```


# Modelize 

 - For each response variable, the approach modelling is 

$Y \sim  pastoreo (Browsing|NoBrowsing)+ Fecha(0|22|29) + zona \times Fecha$ 

- using the "(1|pastoreo:geo_parcela_nombre)" as nested random effects


```{r auxiliarfunctions}
factores <- c("meses", "pastoreo", "geo_parcela_nombre")
rigth_formula <- "pastoreo * meses + (1|pastoreo:geo_parcela_nombre)"

formatAnova <- function(df, yvar){
  df$variable <- yvar
  df$factor <- rownames(df)
  df %>% 
    dplyr::select(fvalue = `F value`, p=`Pr(>F)`, variable, factor) %>%  
  assign(paste0("anova_", yvar), ., inherits = TRUE)
}

formatAnovaGamma <- function(df, yvar){
  df$variable <- yvar
  df$factor <- rownames(df)
  df %>% 
    dplyr::select(fvalue = `F`, p=`Pr(>F)`, variable, factor) %>%  
  assign(paste0("anova_", yvar), ., inherits = TRUE)
}

formatAnovaBeta <- function(df, yvar){
  df$variable <- yvar
  df$factor <- rownames(df)
  df %>% 
    dplyr::select(fvalue = `Chisq`, p=`Pr(>Chisq)`, variable, factor) %>%  
  assign(paste0("anova_", yvar), ., inherits = TRUE)
}


### Posthocs
posthoc.pastoreo <- 
  function(modelo){emmeans(modelo, list(pairwise ~ pastoreo), adjust = "tukey")}
posthoc.meses <- function(modelo){
    emmeans(modelo, list(pairwise ~ meses), adjust = "tukey")} 
posthoc.int <- function(modelo){
    emmeans(modelo, list(pairwise ~ meses | pastoreo), adjust = "tukey")} 



# custom posthocs
tabla_postHoc_meses <- function(modelo){
  posthoc <- emmeans(modelo, list(pairwise ~ meses), adjust = "tukey")
  as.data.frame(posthoc$`pairwise differences of meses`) %>% 
    rename(meses = 1) %>% 
    kbl(digits = 4) %>% 
    kable_paper("hover", full_width = F)
}

tabla_postHoc <- function(modelo){
  posthoc <- emmeans(modelo, ~ pastoreo * meses)
  contrast(posthoc, "consec", simple = "each", combine = TRUE,
           adjust = "tukey") %>% 
  kbl(digits = 4) %>% 
    kable_paper("hover", full_width = F)
}

```


# Humedad 
```{r}
yvar <- "humedad"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::lmer(form, data = df_model)

formula(m)
anova(as_lmerModLmerTest(m))
formatAnova(as.data.frame(anova(as_lmerModLmerTest(m))), yvar)
```

```{r, eval=FALSE}
plot(simulateResiduals(m))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```


# CIC 
```{r}
yvar <- "cic"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::lmer(form, data = df_model)

formula(m)
anova(as_lmerModLmerTest(m))
formatAnova(as.data.frame(anova(as_lmerModLmerTest(m))), yvar)
```

```{r, eval=FALSE}
plot(simulateResiduals(m))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```

# C 
```{r}
yvar <- "c_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::lmer(form, data = df_model)

formula(m)
anova(as_lmerModLmerTest(m))
formatAnova(as.data.frame(anova(as_lmerModLmerTest(m))), yvar)
```

```{r, eval=FALSE}
plot(simulateResiduals(m))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```



# Fe
```{r}
#Gamma
yvar <- "fe_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```





# K
```{r}
#Gamma
yvar <- "k_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```





# Mg
```{r}
#Gamma
yvar <- "mg_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```





# C/N
```{r}
#Gamma
yvar <- "c_n"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```





# MO
```{r}
#Gamma
yvar <- "mo"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```






# pH Agua
```{r}
#Gamma
yvar <- "p_h_agua_eez"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```






# pH KCl
```{r}
#Gamma
yvar <- "p_h_k_cl"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer(form, data = df_model, family= Gamma)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```


# NH4
- No data 
```{r}
soil %>% 
  dplyr::select(n_nh4, meses, pastoreo, fecha) %>% 
  drop_na() %>% group_by(meses, fecha) %>% count()
```

# NO3
- No data 
```{r}
soil %>% 
  dplyr::select(n_no3, meses, pastoreo, fecha) %>% 
  drop_na() %>% group_by(meses, fecha) %>% count()
```




# P
```{r}
#Negative Binomial
yvar <- "p"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- lme4::glmer.nb(form, data = df_model)

formula(m)
ta <- afex::mixed(form, data=df_model)
ta$anova_table
formatAnovaGamma(as.data.frame(ta$anova_table), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```


# N
```{r}
# Beta
yvar <- "n_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- glmmTMB::glmmTMB(formula = form, data = df_model,
                      family = beta_family (link = "logit"))
  
formula(m)
car::Anova(m)
formatAnovaBeta(as.data.frame(car::Anova(m)), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```


# Na
```{r}
# Beta
yvar <- "na_percent"
df_model <- soil %>% dplyr::select(one_of(c(yvar, factores)))
form <- reformulate(rigth_formula, response = yvar)
m <- glmmTMB::glmmTMB(formula = form, data = df_model,
                      family = beta_family (link = "logit"))
  
formula(m)
car::Anova(m)
formatAnovaBeta(as.data.frame(car::Anova(m)), yvar)
```

```{r, eval=FALSE}
model <- lme4::lmer(form, data = df_model)
plot(performance::check_distribution(model)) 
plot(simulateResiduals(model))
```

### *Post-hoc*
```{r}
tabla_postHoc_meses(m)
tabla_postHoc(m)
```




# General Overview
## Mean + SE table 
```{r}
library(plotrix)
library(gtsummary)
soil %>% 
  dplyr::select(pastoreo, meses, one_of("humedad","fe_percent","k_percent","mg_percent","na_percent","n_percent","c_percent","c_n","cic","p","mo","p_h_k_cl","p_h_agua_eez")) %>% 
  tbl_strata(
    strata = pastoreo, 
    .tbl_fun = 
      ~.x %>% 
      tbl_summary(
        by = meses,
        missing = "no",
        type = c(p) ~ "continuous",
       # label = list(ab ~ "Juvenile abundance"),
        statistic = all_continuous() ~ "{mean} ({std.error})",
        digits = everything() ~ 2
      ))
```


## Figures 
```{r}
pos <- position_dodge(0.15)

variables_sel <- c("humedad","fe_percent","k_percent","mg_percent","na_percent","n_percent","c_percent","c_n","cic","p","mo","p_h_k_cl","p_h_agua_eez")


medias <- soil %>% 
  dplyr::select(pastoreo, meses, one_of(variables_sel)) %>% 
  group_by(pastoreo, meses) %>% 
  summarise_at(all_of(variables_sel), ~mean(., na.rm=TRUE)) %>% 
  pivot_longer(all_of(variables_sel), names_to = "var", values_to = "mean")


errores <- soil %>% 
  dplyr::select(pastoreo, meses, one_of(variables_sel)) %>% 
  group_by(pastoreo, meses) %>% 
  summarise_at(all_of(variables_sel), ~plotrix::std.error(., na.rm=TRUE)) %>%  
  pivot_longer(all_of(variables_sel), names_to = "var", values_to = "se")


medias <- medias %>% dplyr::inner_join(errores) %>% 
  mutate(varname = case_when(
    var == "n_percent" ~ "N",
    var == "p" ~ "P",
    var == "k_percent" ~ "K",
    var == "mo" ~ "M.O.",
    var == "c_percent" ~ "C",
    var == "c_n" ~ "C/N",
    var == "humedad" ~ "Humedad",
    # var == "n_nh4" ~ "NH4",
    var == "fe_percent" ~ "Fe",
    var == "mg_percent" ~ "Mg",
    var == "na_percent" ~ "Na",
    var == "cic" ~ "C.I.C.",
    # var == "n_no3" ~ "NO3",
    var == "p_h_k_cl"  ~ "pH (KCl)",
    var == "p_h_agua_eez" ~ "pH (Agua)"
  ))

plotpastoreo <- function(df, ncol, nrow){
    ggplot(df, aes(x=meses, y = mean, 
               color = pastoreo, shape = pastoreo, group=pastoreo)) + 
    geom_point(position = pos, size = 2) + 
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), 
                  position = pos, width = .01) +
    geom_line(position = pos) +
    theme_bw() +
    theme(strip.background = element_rect(fill="white")) + 
    ylab("") + xlab("") +
  facet_wrap(~varname, scales = "free_y", ncol=ncol)
}
```

```{r}
medias %>% 
  filter(var %in% c("n_percent", "p", "k_percent",
                    "mo", "c_percent", "c_n")) %>% 
  plotpastoreo(ncol=2)
```

```{r}
medias %>% 
  filter(var %in% c("cic", "mg_percent", 
                    "na_percent", "humedad",
                    "p_h_k_cl", "p_h_agua_eez")) %>% 
  plotpastoreo(ncol=2)
```

## Anovas table
```{r}
anovas <- bind_rows(anova_c_n, 
                    anova_cic,
                    anova_c_percent,
                    anova_k_percent,
                    anova_humedad, 
                    anova_fe_percent,
                    anova_mg_percent,
                    anova_mo,
                    anova_p, 
                    anova_p_h_agua_eez, 
                    anova_p_h_k_cl, 
                    anova_n_percent, 
                    anova_na_percent) %>% 
  relocate(variable, factor)
rownames(anovas) <- NULL
```

```{r}
library(kableExtra)

anovas %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  pivot_wider(names_from = factor, values_from = c(fvalue, p)) %>% 
  relocate(variable,
           fvalue_pastoreo, p_pastoreo, 
           fvalue_meses, p_meses, 
           `fvalue_pastoreo:meses`, `p_pastoreo:meses`
           ) %>% 
  kbl(col.names = 
        c("Variables", rep(c("F", "p"),3))) %>% 
  add_header_above(c("", "zona" = 2, "fecha" = 2, "zona X fecha" = 2)) %>% 
  kable_styling() 

```

```{r}
anovas %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  pivot_wider(names_from = factor, values_from = c(fvalue, p)) %>% 
  relocate(variable,
           fvalue_pastoreo, p_pastoreo, 
           fvalue_meses, p_meses, 
           `fvalue_pastoreo:meses`, `p_pastoreo:meses`
           ) %>%  
  write_csv(here::here("output/anovas_zona_time_postFire.csv"))
```