---
title: "Impacts of warming and fertilisation on FAPAR"
author: "Ramona Heim"
date: "2022-12-09"
output:
  html_notebook: 
    code_folding: hide
---


This is the R code to the manuscript  
  *Deciduous tundra shrubs shift toward more acquisitive light absorption strategy under climate change treatments*  
  by R. J. Heim, M. Iturrate-Garcia, M. Reji Chacko S. Karsanaev, T. C. Maximov, M. M. P. D. Heijmans, G. Schaepman-Strub

## Data preparation
```{r,results="hide", message=FALSE,warning = FALSE}
tab = read.csv("WE_traits_all.csv")
tab = subset(tab, heating == "hN" & (betula == "Y" | salix == "Y"), select = c(blk, plt, ind, betula, ftr, SLA, LDMC))
tab = na.omit(tab)
tab$species[tab$betula == "Y"] = "betula"
tab$species[tab$betula == "N"] = "salix"
tab$thickness = 1/(tab$SLA*tab$LDMC)
lops_salp = read.csv("FAPAR_salix.csv")
lops_salp$species = "salix"
lops_betn = read.csv("FAPAR_betula.csv")
lops_betn$species = "betula"
lops = as.data.frame(rbind(lops_salp, lops_betn))
agg_final = read.csv("WE_traits_all.csv")
agg_final = subset(agg_final, (betula == "Y" | salix == "Y"), select = c(blk, plt, ind, betula, cable, heating, ftr,  LNC, LPC, CN_ratio, LA, SLA, LDMC))
agg_final$thickness = 1/(agg_final$SLA*agg_final$LDMC)
agg_final$species[agg_final$betula == "Y"] = "betula"
agg_final$species[agg_final$betula == "N"] = "salix"
agg_final$Treatment = agg_final$ftr
agg_final$Plot = agg_final$plt
agg_final <- merge(agg_final, lops, by=c("species","Treatment", "Plot"))
agg_final<-as.data.frame(aggregate(cbind(FAPAR, thickness, SLA, LNC, LPC) ~ 
                               Plot + blk + species + cable + heating + ftr, data=agg_final, FUN=mean))
agg_final$heat_treatment<-as.factor(paste(agg_final$heating, agg_final$cable))
agg_final$fert_treatment<-as.factor(paste(agg_final$ftr))
agg_final$species<-as.factor(paste(agg_final$species))
```

## load packages
```{r,results="hide", message=FALSE,warning = FALSE}
library(brms)
library(rstan)
library(DHARMa)
library(ggplot2)
library(bayesplot)
library(ggExtra)
library(lme4)
library(ggforce)
```

## check data
```{r,results="hide", message=FALSE,warning = FALSE}
hist(agg_final$FAPAR)
```
skewed!

## model preparation
```{r,results="hide", message=FALSE,warning = FALSE}
mod<- lmer(FAPAR ~ species*heat_treatment*fert_treatment+ (1|blk/Plot), data=agg_final)
mod
testDispersion(mod)
simulationOutput <- simulateResiduals(fittedModel = mod, plot = F)
res<-residuals(simulationOutput)
plot(simulationOutput)
plotResiduals(res, agg_final$species)
plotResiduals(res, agg_final$fert_treatment)
plotResiduals(res, agg_final$heat_treatment)
```
Levene test for homogeneity of variance significant for species!

## model
```{r,message=FALSE,warning = FALSE}
prior1 <- c(set_prior ("normal(0,10)",class="b"),
            set_prior ("normal(0,10)",class="Intercept"))
mod_fapar<- brm(bf(FAPAR ~ species*heat_treatment*fert_treatment + (1|blk),
                   (sigma ~ species)),
                family = skew_normal(),
                data   = agg_final, 
                prior=prior1,
                iter   = 10000, 
                warmup = 5000, 
                chains = 4, 
                control=list(adapt_delta=0.99,max_treedepth=15))
mod_fapar
```

## check model
```{r,results="hide", message=FALSE,warning = FALSE}
pp_check(mod_fapar)
pp_check(mod_fapar, type = "ecdf_overlay")
posterior <- posterior <- as.array(mod_fapar)
#plot(mod_fapar,N = 12,ask = F)
```

## posterior dataframe
```{r,results="hide", message=FALSE,warning = FALSE}
newdat<-expand.grid(species=factor(c("betula", "salix"),
                                   levels=levels(agg_final$species)),
                    heat_treatment=factor(c("hN cN" ,"hN cY" ,"hY cY"),
                                          levels=levels(agg_final$heat_treatment)),
                    fert_treatment=factor(c("F" , "NF"),
                                          levels=levels(agg_final$fert_treatment)))
newdat

draws_mod_fapar <- as_draws_array(mod_fapar)
array<-posterior::summarize_draws(draws_mod_fapar)
df <- as_draws_df(draws_mod_fapar)
postfapar_fixed<-df[1:13]
postfapar_fixed<-as.matrix(postfapar_fixed)
postfapar_fixed<-postfapar_fixed[,-c(2)]
fixef<-fixef(mod_fapar)
fixef<-fixef[-c(2,14:18),]


Xmat<-model.matrix(~species*heat_treatment*fert_treatment, data=newdat)
sim<-dim(postfapar_fixed)
nsim<-sim[1]
fitmat<-matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i]<-Xmat %*% postfapar_fixed[i,]
newdat$lower<-(apply(fitmat, 1, quantile, prob=0.025))
newdat$upper<-(apply(fitmat, 1, quantile, prob=0.975))
newdat$lower89<-(apply(fitmat, 1, quantile, prob=0.055))
newdat$upper89<-(apply(fitmat, 1, quantile, prob=0.945))
newdat1 <- as.data.frame((Xmat %*% fixef))
newdat$fit<-as.numeric(newdat1$Estimate)
```
## Hypothesis testing
```{r, message=FALSE,warning = FALSE}
###check hypothesis F > NF
head(postfapar_fixed)
#Betula 
#no warming
sum((postfapar_fixed[,1]>postfapar_fixed[,1]+postfapar_fixed[,5])/nsim) #0.997
#cable
sum((postfapar_fixed[,1]+postfapar_fixed[,3]>postfapar_fixed[,1]+postfapar_fixed[,3]+postfapar_fixed[,5]+postfapar_fixed[,9])/nsim) #0.98
#warming
sum((postfapar_fixed[,1]+postfapar_fixed[,4]>postfapar_fixed[,1]+postfapar_fixed[,4]+postfapar_fixed[,5]+postfapar_fixed[,10])/nsim) #0.8345

#Salix
#no warming
sum((postfapar_fixed[,1]+postfapar_fixed[,2]>postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,5]+postfapar_fixed[,8])/nsim) #0.99525
#cable
sum((postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,3]+postfapar_fixed[,6]>postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,3]+postfapar_fixed[,5]+postfapar_fixed[,6]+postfapar_fixed[,8]+postfapar_fixed[,9]+postfapar_fixed[,11])/nsim) #0.99975
#warming
sum((postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,4]+postfapar_fixed[,7]>postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,4]+postfapar_fixed[,5]+postfapar_fixed[,7]+postfapar_fixed[,8]+postfapar_fixed[,10]+postfapar_fixed[,12])/nsim) #1

###check hypothesis heating no < heating yes
#Betula 
#F
sum((postfapar_fixed[,1]<postfapar_fixed[,1]+postfapar_fixed[,4])/nsim)
#NF
sum((postfapar_fixed[,1]+postfapar_fixed[,5]<postfapar_fixed[,1]+postfapar_fixed[,4]+postfapar_fixed[,5]+postfapar_fixed[,10])/nsim)
#Salix
#F
sum((postfapar_fixed[,1]+postfapar_fixed[,2]<postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,4]+postfapar_fixed[,7])/nsim)
#NF
sum((postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,5]+postfapar_fixed[,8]<postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,4]+postfapar_fixed[,5]+postfapar_fixed[,7]+postfapar_fixed[,8]+postfapar_fixed[,10]+postfapar_fixed[,12])/nsim)

###check hypothesis no treatment < all treatments
#betula
sum((postfapar_fixed[,1]+postfapar_fixed[,5]<postfapar_fixed[,1]+postfapar_fixed[,4])/nsim)
#salix
sum((postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,5]+postfapar_fixed[,8]<postfapar_fixed[,1]+postfapar_fixed[,2]+postfapar_fixed[,4]+postfapar_fixed[,7])/nsim)
```


```{r, message=FALSE,warning = FALSE}
newdat
```

## plot preparation
```{r,results="hide", message=FALSE,warning = FALSE}
legend_title<-"Warming treatment:"
labs <- c("Betula nana", "Salix pulchra")
names(labs) <- c("betula", "salix")

text_nh <- data.frame(
  label = c("0.99", "0.98"),
  species   = c("betula", "salix")
)
text_c <- data.frame(
  label = c("0.95", "1.00"),
  species   = c("betula", "salix")
)
text_yh <- data.frame(
  label = c("0.85", "1.00"),
  species   = c("betula", "salix")
)

plot1<-ggplot(data=agg_final, aes(x = fert_treatment, y = FAPAR, color=heat_treatment))+
  labs(y="FAPAR")+
  labs(x="")+
  facet_grid(~species, labeller = labeller(species=labs))+
  #geom_point(shape=16, alpha=0.8, size=2,
   #          position = position_jitterdodge(jitter.width=0.2,
   #                                          dodge.width = 0.9, jitter.height=0))+
  geom_violin(data=agg_final, aes(x = fert_treatment, y = FAPAR, color=heat_treatment, fill=heat_treatment),alpha=0.2) +
  scale_x_discrete(limits = rev(levels(agg_final$fert_treatment)), labels = c("Not fertilised","Fertilised"))+
  #scale_y_continuous(limits = c(0.832,0.923))+
  #geom_boxplot()+
  geom_linerange(data=newdat, aes(x = fert_treatment, y = fit, ymin=lower, 
                                  ymax=upper, group=heat_treatment), 
                 position = position_dodge(width=0.9), col="black")+
  #geom_linerange(data=newdat, aes(x = fert_treatment, y = fit, ymin=lower89, 
   #                               ymax=upper89, group=heat_treatment), 
   #              position = position_dodge(width=0.9), col="grey", size=1.2)+
  geom_point(data=newdat, aes(x =fert_treatment, y = fit, group=heat_treatment), 
             position = position_dodge(width=0.9), shape=1, size=2 ,col="black")+
  scale_color_manual(legend_title,values=c("#F69C73FF","#E83F3FFF","#4C1D4BFF"), 
                     labels = c("Control without cable", "Control with cable", "Warming"))+
  scale_fill_manual(legend_title,values=c("#F69C73FF","#E83F3FFF","#4C1D4BFF"), 
                     labels = c("Control without cable", "Control with cable", "Warming"))+
  
  #geom_segment(aes(x = 0.75 , y = 0.91, xend = 1.75, yend = 0.91), #col="#F69C73FF") +
  #geom_segment(aes(x =  1, y = 0.915, xend = 2, yend = 0.915), col="#E83F3FFF") +
  #geom_segment(aes(x = 1.25 , y = 0.92, xend = 2.25, yend = 0.92), col="#4C1D4BFF") +
 # geom_segment(aes(x = 0.75 , y = 0.91, xend = 0.75, yend = 0.909), col="#F69C73FF") +
 # geom_segment(aes(x = 1.75 , y = 0.91, xend =1.75, yend = 0.909), col="#F69C73FF") +
  #geom_segment(aes(x =  1, y = 0.915, xend = 1, yend = 0.914), col="#E83F3FFF") +
  #geom_segment(aes(x =  2, y = 0.915, xend = 2, yend = 0.914), col="#E83F3FFF") +
  #geom_segment(aes(x = 1.25 , y = 0.92, xend = 1.25, yend = 0.919), col="#4C1D4BFF") +
  #geom_segment(aes(x = 2.25 , y = 0.92, xend = 2.25, yend = 0.919), col="#4C1D4BFF") +
  #geom_text(data    = text_nh,mapping = aes(x = 1.25, y = 0.912, label = label),
   #         col="#F69C73FF",  size=3)+
  #geom_text(data    = text_c,mapping = aes(x = 1.5, y = 0.917, label = label),
    #        col="#E83F3FFF",  size=3)+
  #geom_text(data    = text_yh,mapping = aes(x = 1.75, y = 0.922, label = label),
    #        col="#4C1D4BFF",  size=3)+

  
  theme_bw()+
  # theme(legend.position="n")+
  theme(legend.background = element_rect(fill = "grey90"),
        legend.position="bottom")+
  guides(colour = guide_legend(override.aes = list(alpha=1)))+
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(face = "italic"))+
  theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm"))
```

```{r, message=FALSE,warning = FALSE}
plot1
```

```{r,results="hide", message=FALSE,warning = FALSE}
png(file="plot_FAPAR_final.png", width=15, height=10, units="cm", pointsize=12,
    res=600)
plot1
dev.off() 

posterior <- as.matrix(postfapar_fixed)
posterior

plot_title <- ggtitle("Posterior distributions",
                      "with medians and 95% intervals")
mcmc_plot<-mcmc_intervals(posterior, pars = c("b_speciessalix", "b_heat_treatmenthNcY", "b_heat_treatmenthYcY",
                                   "b_fert_treatmentNF", "b_speciessalix:heat_treatmenthNcY" ,"b_speciessalix:heat_treatmenthYcY",
                                   "b_speciessalix:fert_treatmentNF", "b_heat_treatmenthNcY:fert_treatmentNF",
                                   "b_heat_treatmenthYcY:fert_treatmentNF" ,"b_speciessalix:heat_treatmenthNcY:fert_treatmentNF",
                                   "b_speciessalix:heat_treatmenthYcY:fert_treatmentNF"),
               prob = 0.8, prob_outer = 0.95) + plot_title +
  theme_bw()+
  # theme(legend.position="n")+
  theme(legend.position=c(0.15, 0.15))+
  guides(colour = guide_legend(override.aes = list(alpha=1)))+
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(face = "italic"))+
  theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm"))

png(file="plot_MCMC.png", width=15, height=15, units="cm", pointsize=12,res=600)
mcmc_plot
dev.off() 
```
