##Supplementary file 13_Deer and non-tree cover effects on tree cover in gap plots (Figure 3)

## R packages
library(glmmTMB)
library(emmeans)
library(DHARMa)

##Loading data file from "vegetation cover data" (Supplementary file 10)

##selection of clearcut plots
clearcut<-Vegetation_cover_data[13:24,] 

##Model testing
M0<-glmmTMB(Treecover.adj~groundvegetation+Deer+groundvegetation*Deer+(1|location/plot),data = clearcut, family = beta_family(link = "logit"))
summary(M0)

M1<-glmmTMB(Treecover.adj~groundvegetation+Deer+(1|location/plot),data = clearcut, family = beta_family(link = "logit"))
M2<-glmmTMB(Treecover.adj~groundvegetation+Deer+(1|plot),data = clearcut, family = beta_family(link = "logit")) 
M3<-glmmTMB(Treecover.adj~groundvegetation+Deer+groundvegetation*Deer+(1|plot),data = clearcut, family = beta_family(link = "logit"))  

##Testing GLMM residual diagnostics using DHARMa
##Model is rejected if significant deviation is detected
simulateResiduals(M0, plot = TRUE)
simulateResiduals(M1, plot = TRUE)
simulateResiduals(M2, plot = TRUE)
simulateResiduals(M3, plot = TRUE)

##Comparison of models without significant deviation of residuals
anova(M0,M2)

##Best-fit model (M2) selected based on AIC and BIC values

summary(M2)