################################################ # Position analyses ################################################ set.seed(1206) rm(list = ls()) model_data <- read_excel("Output/Data_for_modelling/model_data.xlsx") names(model_data) trtlev <- c("shady_exclosure", "shady_control", "sunny_exclosure", "sunny_control") model_data$treatment <- factor(model_data$treatment, levels = trtlev, labels = trtlev) clev <- c("shady", "sunny") model_data$canopy <- factor(model_data$canopy, levels = clev, labels = clev) flev <- c("exclosure", "control") model_data$fence <- factor(model_data$fence, levels = flev, labels = flev) str(model_data) table(model_data$treatment, with(model_data, interaction(canopy, fence))) library(readxl) Position <- read_excel("Input/Position.xlsx") model_data <- cbind(model_data, Position) names(model_data) # species richness of herb layer Richness_pos <- glmmTMB(herb_q0 ~ pos + pos:canopy + pos:fence + pos:pH + (1|plot), family = nbinom2, data = model_data, control = glmmTMBControl( optimizer = optim, optArgs = list(method = "BFGS", maxit = 1000))) # res <- simulateResiduals(fittedModel = Richness_pos, plot = T) # testDispersion(res) # testZeroInflation(res) # plotResiduals(res) summary(Richness_pos) # Estimate Std. Error z value Pr(>|z|) # (Intercept) 0.50263 0.38472 1.306 0.191390 # possp 0.84219 0.49899 1.688 0.091455 . # pospm:canopysunny -0.34432 0.11340 -3.036 0.002395 ** # possp:canopysunny -0.03352 0.10880 -0.308 0.758041 # pospm:fencecontrol 0.16382 0.09734 1.683 0.092397 . # possp:fencecontrol 0.09924 0.09681 1.025 0.305302 # pospm:pH 0.24577 0.06678 3.680 0.000233 *** # possp:pH 0.08660 0.05932 1.460 0.144320 car::Anova(Richness_pos, type = "III") # Chisq Df Pr(>Chisq) # (Intercept) 1.7069 1 0.1913900 # pos 2.8486 1 0.0914554 . # pos:canopy 9.2205 2 0.0099496 ** # pos:fence 4.3653 2 0.1127440 # pos:pH 14.9986 2 0.0005535 *** Richness_pos.T <- glmmTMB(herb_q0 ~ pos + pos:treatment + pos:pH + (1|plot), family = nbinom2, data = model_data, control = glmmTMBControl( optimizer = optim, optArgs = list(method = "BFGS", maxit = 1000))) res <- simulateResiduals(fittedModel = Richness_pos.T, plot = T) testDispersion(res) testZeroInflation(res) plotResiduals(res) summary(Richness_pos.T) # Estimate Std. Error z value Pr(>|z|) # (Intercept) 0.52815 0.37813 1.397 0.162490 # possp 0.80214 0.49532 1.619 0.105350 # pospm:treatmentshady_control 0.09012 0.10839 0.831 0.405701 # possp:treatmentshady_control 0.21741 0.11131 1.953 0.050802 . # pospm:treatmentsunny_exclosure -0.54187 0.18621 -2.910 0.003614 ** # possp:treatmentsunny_exclosure 0.15129 0.13870 1.091 0.275400 # pospm:treatmentsunny_control -0.13708 0.14666 -0.935 0.349935 # possp:treatmentsunny_control -0.06616 0.16493 -0.401 0.688325 # pospm:pH 0.24970 0.06545 3.815 0.000136 *** # possp:pH 0.08006 0.05815 1.377 0.168555 library(multcomp) names(fixef(Richness_pos.T)$cond) K <- rbind( "sp - pm (shady_exclosure)" = c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0), "sp - pm (shady_control)" = c(0, 1, -1, 1, 0, 0, 0, 0, 0, 0), "sp - pm (sunny_exclosure)" = c(0, 1, 0, 0, -1, 1, 0, 0, 0, 0), "sp - pm (sunny_control)" = c(0, 1, 0, 0, 0, 0, -1, 1, 0, 0) ) colnames(K) <- names(fixef(Richness_pos.T)$cond) pos_g <- glht(Richness_pos.T, linfct = K) summary(pos_g, test = adjusted("holm")) # Estimate Std. Error z value Pr(>|z|) # sp - pm (shady_exclosure) == 0 0.8021 0.4953 1.619 0.2047 # sp - pm (shady_control) == 0 0.9294 0.4927 1.887 0.1777 # sp - pm (sunny_exclosure) == 0 1.4953 0.5483 2.727 0.0255 * # sp - pm (sunny_control) == 0 0.8731 0.5345 1.633 0.2047 # functional diversity ---- FD_pos <- glmmTMB(ses ~ pos + pos:canopy + pos:fence + pos:pH + (1|plot), data = model_data) # res <- simulateResiduals(fittedModel = FD_pos, plot = T) # testDispersion(res) # testZeroInflation(res) # plotResiduals(res) summary(FD_pos) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -2.8650 0.8835 -3.243 0.00118 ** # possp 2.3631 1.0725 2.203 0.02756 * # pospm:canopysunny 0.3942 0.2674 1.474 0.14036 # possp:canopysunny 0.1937 0.2673 0.725 0.46864 # pospm:fencecontrol -0.2609 0.2404 -1.085 0.27772 # possp:fencecontrol -0.3849 0.2389 -1.611 0.10712 # pospm:pH 0.2592 0.1577 1.644 0.10022 # possp:pH -0.1364 0.1406 -0.970 0.33200 FD_pos <- glmmTMB(ses ~ pos + canopy + fence + pH + (1|plot), data = model_data) # res <- simulateResiduals(fittedModel = FD_pos, plot = T) # testDispersion(res) # testZeroInflation(res) # plotResiduals(res) summary(FD_pos) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -1.56381 0.65961 -2.371 0.0177 * # possp 0.01762 0.14639 0.120 0.9042 # canopysunny 0.29510 0.21548 1.369 0.1709 # fencecontrol -0.29886 0.14634 -2.042 0.0411 * # pH 0.03680 0.11601 0.317 0.7511 # shrub layer ---- ShrubLayer_pos <- glmmTMB(CoverShrubLayer ~ pos + pos:canopy + pos:fence + pos:pH + (1|plot), family = ordbeta(link = "logit"), data = model_data) # res <- simulateResiduals(fittedModel = ShrubLayer_pos, plot = T) # testDispersion(res) # testZeroInflation(res) # plotResiduals(res) summary(ShrubLayer_pos) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -3.52514 1.17939 -2.989 0.00280 ** # possp 0.05709 1.52709 0.037 0.97018 # pospm:canopysunny 2.16710 0.31550 6.869 6.48e-12 *** # possp:canopysunny 2.32722 0.31264 7.444 9.79e-14 *** # pospm:fencecontrol -0.31572 0.27319 -1.156 0.24781 # possp:fencecontrol -0.85230 0.29878 -2.853 0.00434 ** # pospm:pH 0.28005 0.20093 1.394 0.16340 # possp:pH 0.35077 0.19840 1.768 0.07707 . car::Anova(ShrubLayer_pos, type = "III") # Chisq Df Pr(>Chisq) # (Intercept) 8.9338 1 0.002799 ** # pos 0.0014 1 0.970176 # pos:canopy 72.3401 2 < 2.2e-16 *** # pos:fence 11.7470 2 0.002813 ** # pos:pH 4.4078 2 0.110370 ShrubLayer_pos.T <- glmmTMB(CoverShrubLayer ~ pos + pos:treatment + pos:pH + (1|plot), family = ordbeta(link = "logit"), data = model_data) summary(ShrubLayer_pos.T) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -3.6366 1.1918 -3.051 0.002278 ** # possp 0.1803 1.5423 0.117 0.906945 # pospm:treatmentshady_control -0.1412 0.3649 -0.387 0.698736 # possp:treatmentshady_control -0.9586 0.4053 -2.365 0.018026 * # pospm:treatmentsunny_exclosure 2.3786 0.4320 5.507 3.66e-08 *** # possp:treatmentsunny_exclosure 2.2161 0.3831 5.785 7.27e-09 *** # pospm:treatmentsunny_control 1.8321 0.3956 4.631 3.64e-06 *** # possp:treatmentsunny_control 1.5276 0.3935 3.882 0.000104 *** # pospm:pH 0.2850 0.2013 1.416 0.156833 # possp:pH 0.3553 0.1994 1.782 0.074711 . names(fixef(ShrubLayer_pos.T)$cond) K <- rbind( "sp - pm (shady_exclosure)" = c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0), "sp - pm (shady_control)" = c(0, 1, -1, 1, 0, 0, 0, 0, 0, 0), "sp - pm (sunny_exclosure)" = c(0, 1, 0, 0, -1, 1, 0, 0, 0, 0), "sp - pm (sunny_control)" = c(0, 1, 0, 0, 0, 0, -1, 1, 0, 0) ) colnames(K) <- names(fixef(ShrubLayer_pos.T)$cond) pos_shrub <- glht(ShrubLayer_pos.T, linfct = K) summary(pos_shrub, test = adjusted("holm")) # Estimate Std. Error z value Pr(>|z|) # sp - pm (shady_exclosure) == 0 0.18029 1.54234 0.117 1 # sp - pm (shady_control) == 0 -0.63707 1.56989 -0.406 1 # sp - pm (sunny_exclosure) == 0 0.01779 1.59320 0.011 1 # sp - pm (sunny_control) == 0 -0.12420 1.56876 -0.079 1