#Code for analyses of data for the manuscript 
#"Effect of arginine-phosphate addition on early survival and growth of Scots pine, Norway spruce and silver birch"
# Manuscript by B. Häggström, R.Lutter, T. Lundmark, F.Sjödin, A. Nordin
#Code by B.Häggström
#R version: RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) 
library(dplyr)
library(lmerTest)
library(pbkrtest)
library(lme4)
library(car)
library(emmeans)
library(ggplot2)
library(ggthemes)
library(Rmisc)
library(ggpubr)
library(ggpmisc)
library(multcomp)
library(multcompView)
library(sjPlot)
library(grid)
library(gridExtra)
require(glmmTMB)

setwd(".")

# We will need the following function for computation of Pearson chi-square/DF
# Source: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
#download and load following data files from depository,
#check that decimal is set to "comma" when importing, otherwise the measurements will not be numeric
effects_APtreatment_data_S_second_season <-
        read.csv("effects_APtreatment_data_S_second_season.csv", sep = ";", dec = ",") 
effects_APtreatment_data_NE_second_season <-
        read.csv("effects_APtreatment_data_NE_second_season.csv", sep = ";", dec = ",") 
effects_APtreatment_data_NW_second_season <- 
        read.csv("effects_APtreatment_data_NW_second_season.csv", sep = ";", dec = ",") 

#If you want to do typ-III test ANOVA, you have to have your factors coded the right way.
#The parametrisation "contr.treatment" is default,has to be changed to "contr.sum"
#prepare factors for analysis and enabling type III:
effects_APtreatment_data_S_second_season$AP_treatment<-as.factor(effects_APtreatment_data_S_second_season$AP_treatment)
effects_APtreatment_data_S_second_season$Species<-as.factor(effects_APtreatment_data_S_second_season$Species)
effects_APtreatment_data_S_second_season$Composition<-as.factor(effects_APtreatment_data_S_second_season$Composition)
contrasts(effects_APtreatment_data_S_second_season$AP_treatment)<-contr.sum
contrasts(effects_APtreatment_data_S_second_season$Species)<-contr.sum
contrasts(effects_APtreatment_data_S_second_season$Composition)<-contr.sum
effects_APtreatment_data_NE_second_season$AP_treatment<-as.factor(effects_APtreatment_data_NE_second_season$AP_treatment)
effects_APtreatment_data_NE_second_season$Species<-as.factor(effects_APtreatment_data_NE_second_season$Species)
effects_APtreatment_data_NE_second_season$Composition<-as.factor(effects_APtreatment_data_NE_second_season$Composition)
contrasts(effects_APtreatment_data_NE_second_season$AP_treatment)<-contr.sum
contrasts(effects_APtreatment_data_NE_second_season$Species)<-contr.sum
contrasts(effects_APtreatment_data_NE_second_season$Composition)<-contr.sum
effects_APtreatment_data_NW_second_season$AP_treatment<-as.factor(effects_APtreatment_data_NW_second_season$AP_treatment)
effects_APtreatment_data_NW_second_season$Species<-as.factor(effects_APtreatment_data_NW_second_season$Species)
effects_APtreatment_data_NW_second_season$Composition<-as.factor(effects_APtreatment_data_NW_second_season$Composition)
contrasts(effects_APtreatment_data_NW_second_season$AP_treatment)<-contr.sum
contrasts(effects_APtreatment_data_NW_second_season$Species)<-contr.sum
contrasts(effects_APtreatment_data_NW_second_season$Composition)<-contr.sum
#add column defining binary outcome of alive and dead seedlings for S:
effects_APtreatment_data_S_second_season$Alive[effects_APtreatment_data_S_second_season$Status=="dead"]<-"0"
effects_APtreatment_data_S_second_season$Alive[effects_APtreatment_data_S_second_season$Status=="missing"]<-"0"
effects_APtreatment_data_S_second_season$Alive[effects_APtreatment_data_S_second_season$Status=="vital"]<-"1"
effects_APtreatment_data_S_second_season$Alive[effects_APtreatment_data_S_second_season$Status=="damaged"]<-"1"
effects_APtreatment_data_S_second_season$Alive<-as.numeric(effects_APtreatment_data_S_second_season$Alive)
#Model south (S) site survival
SouthSurvivalModel<-glmer(Alive~AP_treatment*Species+(1|Species:Composition),data=effects_APtreatment_data_S_second_season,family = binomial(link = logit))
summary(SouthSurvivalModel)
#check for overdispersion (ratio should be close to 1):
overdisp_fun(SouthSurvivalModel)
#check random effect
plot_model(SouthSurvivalModel, type = c("diag"))
ggsave("SouthSurvivalModel.pdf")
#Analysis of variance of the model:
Anova(SouthSurvivalModel,type = "III")
#make data frame with letters to distinguish significance
SouthSurvivalemmeansPROBABILITY<-as.data.frame(cld(emmeans(SouthSurvivalModel, ~ AP_treatment:Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotSouthSurvival<-ggplot(SouthSurvivalemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))+
  ggtitle("South")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotSouthSurvival

#add column defining binary outcome of alive and dead seedlings for NE:
effects_APtreatment_data_NE_second_season$Alive[effects_APtreatment_data_NE_second_season$Status=="dead"]<-"0"
effects_APtreatment_data_NE_second_season$Alive[effects_APtreatment_data_NE_second_season$Status=="missing"]<-"0"
effects_APtreatment_data_NE_second_season$Alive[effects_APtreatment_data_NE_second_season$Status=="vital"]<-"1"
effects_APtreatment_data_NE_second_season$Alive[effects_APtreatment_data_NE_second_season$Status=="damaged"]<-"1"
effects_APtreatment_data_NE_second_season$Alive<-as.numeric(effects_APtreatment_data_NE_second_season$Alive)
#Model northeast (NE) site survival
NortheastSurvivalModel<-glmer(Alive~AP_treatment*Species+(1|Species:Composition),data=effects_APtreatment_data_NE_second_season,family = binomial(link = logit))
summary(NortheastSurvivalModel)
#check for overdispersion (ratio should be close to 1):
overdisp_fun(NortheastSurvivalModel)
#check random effect
plot_model(NortheastSurvivalModel, type = c("diag"))
ggsave("NortheastSurvivalModel.pdf")
#Analysis of variance of the model:
Anova(NortheastSurvivalModel,type = "III")
#make data frame with letters to distinguish significance
NESurvivalemmeansPROBABILITY<-as.data.frame(cld(emmeans(NortheastSurvivalModel, ~ AP_treatment:Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotNESurvival<-ggplot(NESurvivalemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northeast (NE)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotNESurvival

#add column defining binary outcome of alive and dead seedlings for NW:
effects_APtreatment_data_NW_second_season$Alive[effects_APtreatment_data_NW_second_season$Status=="dead"]<-"0"
effects_APtreatment_data_NW_second_season$Alive[effects_APtreatment_data_NW_second_season$Status=="missing"]<-"0"
effects_APtreatment_data_NW_second_season$Alive[effects_APtreatment_data_NW_second_season$Status=="vital"]<-"1"
effects_APtreatment_data_NW_second_season$Alive[effects_APtreatment_data_NW_second_season$Status=="damaged"]<-"1"
effects_APtreatment_data_NW_second_season$Alive<-as.numeric(effects_APtreatment_data_NW_second_season$Alive)
#Model northwest (NW) site survival
NorthwestSurvivalModel<-glmer(Alive~AP_treatment+Species+(1|Species:Composition),data=effects_APtreatment_data_NW_second_season,family = binomial(link = logit))
#boundary (singular) fit: see ?isSingular
#→ No big problem. This usually means that some
#of the variances were estimated to zero.
summary(NorthwestSurvivalModel)
#check for overdispersion (ratio should be close to 1):
overdisp_fun(NorthwestSurvivalModel)
#check random effect
plot_model(NorthwestSurvivalModel, type = c("diag")) #the random effect add no variance
ggsave("NorthwestSurvivalModel.pdf")
#checking result with random effect removed gives very similar outcome:
NorthwestSurvivalModel.2<-glm(Alive~AP_treatment+Species,data=effects_APtreatment_data_NW_second_season,family = binomial(link = logit))
#Analysis of variance of the model:
Anova(NorthwestSurvivalModel.2)
#Decision made to keep random effect in model to make the method description in manuscript coherent

#make data frame with letters to distinguish significance
NWSurvivalemmeansPROBABILITY<-as.data.frame(cld(emmeans(NorthwestSurvivalModel, ~ AP_treatment+Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotNWSurvival<-ggplot(NWSurvivalemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northwest (NW)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotNWSurvival

#Present all three sites together
grid.arrange(PlotSouthSurvival,PlotNESurvival,PlotNWSurvival,ncol=3)
SURVIVALPLOTS <- gridExtra::arrangeGrob(PlotSouthSurvival,PlotNESurvival,PlotNWSurvival,ncol=3,
                                        left=grid::textGrob("Survival (%)", gp=gpar(fontsize=16),rot=90))
grid::grid.newpage() 
grid::grid.draw(SURVIVALPLOTS)
ggsave("survivalplots.pdf", SURVIVALPLOTS, width=9)
#(to show as in manuscript, select copy plot under export and change width to 1000)

#Analysis of pine weevil damage on living seedlings at the south site
# Create subset with only living seedlings on south site:
SouthAlive<-effects_APtreatment_data_S_second_season%>%filter(effects_APtreatment_data_S_second_season$Alive=="1")
#add column defining binary outcome of pineweevil-attacked vs non attacked seedlings:
SouthAlive$Weevil<-"0"
SouthAlive$Weevil[SouthAlive$Damage_cause=="Pineweevil"]<-"1"
#set column "Weevil" as numeric to enable analysis
SouthAlive$Weevil<-as.numeric(SouthAlive$Weevil)
#Model weevil damaged 
SouthWeevilModel<-glmer(Weevil~AP_treatment*Species+(1|Species:Composition),data=SouthAlive,family = binomial(link = logit))
#check for overdispersion (ratio should be close to 1):
overdisp_fun(SouthWeevilModel)
#check random effect
pdf("SouthWeevilModel.pdf")
plot_model(SouthWeevilModel, type = c("diag"))
dev.off()
#Analysis of variance of the model:
Anova(SouthWeevilModel,type="III")
#make data frame with letters to distinguish significance
SouthWeevilemmeansPROBABILITY<-as.data.frame(cld(emmeans(SouthWeevilModel, ~ AP_treatment:Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotSouthWeevil<-ggplot(SouthWeevilemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  ggtitle("South, pine weevil")+guides(fill = guide_legend(title = "AP treatment"))+
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9))+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))
PlotSouthWeevil

#Analysis of significant damage (only south site pine weevil + browsed and NW browsed)

#Analysis of browsing damage on the south site
#add column defining binary outcome of browsed vs non non-browsed seedlings:
SouthAlive$browsed<-"0"
SouthAlive$browsed[SouthAlive$Damage_cause=="Browsed"]<-"1"
#set column "Weevil" as numeric to enable analysis
SouthAlive$browsed<-as.numeric(SouthAlive$browsed)
#Model browsed on south site
SouthBrowsedModel<-glmer(browsed~AP_treatment*Species+(1|Species:Composition),data=SouthAlive,family = binomial(link = logit))
summary(SouthBrowsedModel)
#check for overdispersion (ratio should be close to 1):
overdisp_fun(SouthBrowsedModel) 
#check random effect
pdf("SouthBrowsedModel.pdf")
plot_model(SouthBrowsedModel, type = c("diag"))
dev.off()
#Analysis of variance of the model:
Anova(SouthBrowsedModel,type="III")
#make data frame with letters to distinguish significance
SouthBrowsedemmeansPROBABILITY<-as.data.frame(cld(emmeans(SouthBrowsedModel, ~ AP_treatment:Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotSouthBrowsed<-ggplot(SouthBrowsedemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  ggtitle("South, browsed")+guides(fill = guide_legend(title = "AP treatment"))+
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9))+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))
PlotSouthBrowsed

#Analysis of browsing damage on the northwest site
# Create subset with only living seedlings on northwest site:
NWAlive<-effects_APtreatment_data_NW_second_season%>%filter(effects_APtreatment_data_NW_second_season$Alive=="1")
#add column defining binary outcome of browsed vs non non-browsed seedlings:
NWAlive$browsed<-"0"
NWAlive$browsed[NWAlive$Damage_cause=="Browsed"]<-"1"
#set column "Weevil" as numeric to enable analysis
NWAlive$browsed<-as.numeric(NWAlive$browsed)
#Model browsed on NW site
NWBrowsed<-glmer(browsed~AP_treatment+Species+(1|Species:Composition),data=NWAlive,family = binomial(link = logit))
summary(NWBrowsed)
#check for overdispersion (ratio should be close to 1):
overdisp_fun(NWBrowsed)
#check random effect
pdf("NWBrowsedModel.pdf")
plot_model(NWBrowsed, type = c("diag"))
dev.off()
#Analysis of variance of the model:
Anova(NWBrowsed)
#make data frame with letters to distinguish significance
NWBrowsedemmeansPROBABILITY<-as.data.frame(cld(emmeans(NWBrowsed, ~ AP_treatment+Species,type="response"),Letters = letters,sort=FALSE))
#make plot of the data from the model
PlotNWBrowsed<-ggplot(NWBrowsedemmeansPROBABILITY,aes(x=Species,y=prob*100,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  ggtitle("Northwest (NW), browsed")+guides(fill = guide_legend(title = "AP treatment"))+
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=25,ymax=100)+
  theme_pubr(base_size = 14)+scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme(axis.title.x = element_blank(),axis.title.y = element_blank())+
  geom_text(aes(label = .group,y=prob*100+SE*100), vjust = -0.5,position=position_dodge(0.9))+
  geom_errorbar(aes(ymin=prob*100-SE*100,ymax=prob*100+SE*100),width = 0.12,position=position_dodge(0.9))
PlotNWBrowsed

#Present all significant damage together:
grid.arrange(PlotSouthWeevil,PlotSouthBrowsed,PlotNWBrowsed,ncol=3)
DAMAGEDPLOTS <- gridExtra::arrangeGrob(PlotSouthWeevil,PlotSouthBrowsed,PlotNWBrowsed,ncol=3,
                                       left=grid::textGrob("Damaged out of living (%)", gp=gpar(fontsize=16),rot=90))
grid::grid.newpage() 
grid::grid.draw(DAMAGEDPLOTS)
ggsave("damagedplots.pdf", DAMAGEDPLOTS, width=10)
#(to show as in manuscript, select copy plot under export and change width to 1000)

#For analysis of growth variables, subset vital seedlings and remove any unmeasured seedlings:
#remove missing data from height
#also check if height and diameter are set as numeric variables
SouthVital<-SouthAlive%>%filter(Status=="vital")
SouthVital<-SouthVital[!(is.na(SouthVital$Height_cm)),]
SouthVital<-SouthVital[!(is.na(SouthVital$Diameter_mm)),]

NEVital<-effects_APtreatment_data_NE_second_season%>%filter(Status=="vital")
NEVital<-NEVital[!(is.na(NEVital$Height_cm)),]
NEVital<-NEVital[!(is.na(NEVital$Diameter_mm)),]

NWVital<-NWAlive%>%filter(Status=="vital")
NWVital<-NWVital[!(is.na(NWVital$Height_cm)),]
NWVital<-NWVital[!(is.na(NWVital$Diameter_mm)),]

#Model fitting, south site, height with species composition as random factor, lmer
SouthHeight<-lmer(Height_cm~AP_treatment+Species+(1|Species:Composition),data=SouthVital)
summary(SouthHeight)
#check residuals
pdf("SouthHeight.pdf")
plot(SouthHeight)
dev.off()
pdf("SouthHeightQQ.pdf")
qqnorm(residuals(SouthHeight))
qqline(residuals(SouthHeight))
dev.off()
#Analysis of variance of the model:
Anova(SouthHeight)
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(SouthHeight,ddf="Kenward-Roger",type="II")
#since the AP-treatment was not significant and species highly significant in both
#models, reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and damage analyses

#make data frame with letters to distinguish significance
SouthHeightemmeans <- as.data.frame(cld(emmeans(SouthHeight,~ AP_treatment+Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotHeightSouth<-ggplot(SouthHeightemmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("South")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=130)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotHeightSouth

#Model fitting, northeast (NE) site, height growth with species composition as random factor, lmer
NEHeight<-lmer(Height_cm~AP_treatment*Species+(1|Species:Composition), data=NEVital)
#boundary (singular) fit: see ?isSingular
#→ No big problem. This usually means that some
#of the variances were estimated to zero.
summary(NEHeight)
#check residuals
pdf("NEHeight.pdf")
plot(NEHeight)
dev.off()
pdf("NEHeightQQ.pdf")
qqnorm(residuals(NEHeight))
qqline(residuals(NEHeight))
dev.off()
#Analysis of variance of the model:
Anova(NEHeight,type="III")
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(NEHeight,ddf="Kenward-Roger",type="III")
#since the AP-treatment was not significant and species highly significant in both
#models, reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and danage analyses

#make data frame with letters to distinguish significance
NEHeightemmeans <- as.data.frame(cld(emmeans(NEHeight,~ AP_treatment:Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotHeightNE<-ggplot(NEHeightemmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northeast (NE)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=130)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotHeightNE

#Model fitting, northeast (NW) site, height growth with species composition as random factor, lmer
NWHeight<-lmer(Height_cm~AP_treatment*Species+(1|Species:Composition), data=NWVital)
#boundary (singular) fit: see ?isSingular
#→ No big problem. This usually means that some
#of the variances were estimated to zero.
summary(NWHeight)
#check residuals
pdf("NWHeight.pdf")
plot(NWHeight)
dev.off()
pdf("NWHeightQQ.pdf")
qqnorm(residuals(NWHeight))
qqline(residuals(NWHeight))
dev.off()
#Analysis of variance of the model:
Anova(NWHeight,type="III")
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(NWHeight,ddf="Kenward-Roger",type="III")
#since the AP-treatment was not significant and species highly significant in both
#models, reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and danage analyses

#make data frame with letters to distinguish significance
NWHeightemmeans <- as.data.frame(cld(emmeans(NWHeight,~ AP_treatment:Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotHeightNW<-ggplot(NWHeightemmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northwest (NW)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=130)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotHeightNW

#Present all significant damage together:
grid.arrange(PlotHeightSouth,PlotHeightNE,PlotHeightNW,ncol=3)
HEIGHTPLOTS <- gridExtra::arrangeGrob(PlotHeightSouth,PlotHeightNE,PlotHeightNW,ncol=3,
                                      left=grid::textGrob("Height (cm)", gp=gpar(fontsize=16),rot=90))
grid::grid.newpage() 
grid::grid.draw(HEIGHTPLOTS)
ggsave("heightplots.pdf", HEIGHTPLOTS, width=9)
#(to show as in manuscript, select copy plot under export and change width to 1000)

#Model fitting, south site, diameter at stem base  with species composition as random factor, lmer
SouthDiameter<-lmer(Diameter_mm~AP_treatment+Species+(1|Species:Composition),data=SouthVital)
summary(SouthDiameter)
#check residuals
pdf("SouthDiameter.pdf")
plot(SouthDiameter)
dev.off()
pdf("SouthDiameterQQ.pdf")
qqnorm(residuals(SouthDiameter))
qqline(residuals(SouthDiameter))
dev.off()
#Analysis of variance of the model:
Anova(SouthDiameter)
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(SouthDiameter,ddf="Kenward-Roger",type="II")
#since the AP-treatment was not significant and species significant in both
#models, albeit strongly in the Anova and moderately in the anova, 
#reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and damage analyses

#make data frame with letters to distinguish significance
SouthDiameteremmeans <- as.data.frame(cld(emmeans(SouthDiameter,~ AP_treatment+Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotDiameterSouth<-ggplot(SouthDiameteremmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("South")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=15)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotDiameterSouth

#Model fitting, northeast (NE) site, height growth with species composition as random factor, lmer
NEDiameter<-lmer(Diameter_mm~AP_treatment*Species+(1|Species:Composition), data=NEVital)
summary(NEDiameter)
#check residuals
pdf("NEDiameter.pdf")
plot(NEDiameter)
dev.off()
pdf("NEDiameterQQ.pdf")
qqnorm(residuals(NEDiameter))
qqline(residuals(NEDiameter))
dev.off()
#Analysis of variance of the model:
Anova(NEDiameter,type="III")
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(NEDiameter,ddf="Kenward-Roger",type="III")
#since the AP-treatment was not significant and species highly significant in both
#models, albeit strongly for species in the Anova and moderately in the anova,
#reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and danage analyses

#make data frame with letters to distinguish significance
NEDiameteremmeans <- as.data.frame(cld(emmeans(NEDiameter,~ AP_treatment:Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotDiameterNE<-ggplot(NEDiameteremmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northeast (NE)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=15)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotDiameterNE

#Model fitting, northeast (NW) site, height growth with species composition as random factor, lmer
NWDiameter<-lmer(Diameter_mm~AP_treatment*Species+(1|Species:Composition), data=NWVital)
summary(NWDiameter)
#check residuals
pdf("NWDiameter.pdf")
plot(NWDiameter)
dev.off()
pdf("NWDiameterQQ.pdf")
qqnorm(residuals(NWDiameter))
qqline(residuals(NWDiameter))
dev.off()
#Analysis of variance of the model:
Anova(NWDiameter,type="III")
#checking with basic anova to see if there is any big difference with sum square instead of chisquare:
anova(NWDiameter,ddf="Kenward-Roger",type="III")
#since the AP-treatment was not significant and species highly significant in both
#models, albeit strongly for species in the Anova and moderately in the anova,
#reporting the result from Anova with chisquare was chosen to make tables uniform
#with results from survival and danage analyses

#make data frame with letters to distinguish significance
NWDiameteremmeans <- as.data.frame(cld(emmeans(NWDiameter,~ AP_treatment:Species),Letters = letters, sort = FALSE))
#make plot of the data from the model
PlotDiameterNW<-ggplot(NWDiameteremmeans,aes(x=Species,y=emmean,fill=AP_treatment))+
  geom_bar(stat="identity",color="black",position=position_dodge(0.9),width=0.8)+
  geom_errorbar(aes(ymin=emmean-SE,ymax=emmean+SE),width = 0.12,position=position_dodge(0.9))+
  ggtitle("Northwest (NW)")+guides(fill = guide_legend(title = "AP treatment"))+
  geom_text(aes(label = .group,y=emmean+SE), vjust = -0.5,position=position_dodge(0.9)) +
  scale_fill_manual(values = c("lightgrey","white"))+expand_limits(ymin=0,ymax=15)+
  scale_x_discrete(labels=c("Pine","Spruce","Birch"))+
  theme_pubr(base_size = 14)+theme(axis.title.x = element_blank(),axis.title.y = element_blank())
PlotDiameterNW

#Present all significant damage together:
grid.arrange(PlotDiameterSouth,PlotDiameterNE,PlotDiameterNW,ncol=3)
DIAMETERPLOTS <- gridExtra::arrangeGrob(PlotDiameterSouth,PlotDiameterNE,PlotDiameterNW,ncol=3,
                                        left=grid::textGrob("Diameter (mm)", gp=gpar(fontsize=16),rot=90))
grid::grid.newpage() 
grid::grid.draw(DIAMETERPLOTS)
ggsave("diameterplots.pdf", DIAMETERPLOTS, width=9)
#(to show as in manuscript, select copy plot under export and change width to 1000)

cat("\n")
sessionInfo()