Deadwood decomposition mediated by insects influences global carbon fluxes; Seibold et al. (2021)
This document reproduces tables and figures.
This HTML file was generated by
# library("knitr")
# spin("Seibold_et_al_2021.R")
Packages
library("lme4")
library("multcomp")
library("lattice")
library("gridExtra")
library("lme4")
library("colorspace")
library("ggplot2")
library("RColorBrewer")
library("ggpubr")
library("patchwork")
library("vegan")
library("car")
library("MuMIn")
library("tripack")
library("viridis")
Results were obtained in this environment
sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=de_CH.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridis_0.5.1 viridisLite_0.3.0 tripack_1.3-9.1 MuMIn_1.43.17
## [5] car_3.0-10 carData_3.0-4 vegan_2.5-7 permute_0.9-5
## [9] patchwork_1.1.1 ggpubr_0.4.0 RColorBrewer_1.1-2 ggplot2_3.3.3
## [13] colorspace_2.0-0 gridExtra_2.3 lattice_0.20-41 multcomp_1.4-16
## [17] TH.data_1.0-10 MASS_7.3-53 survival_3.2-7 mvtnorm_1.1-1
## [21] lme4_1.1-26 Matrix_1.3-2 knitr_1.31
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.1.2 splines_4.0.4 assertthat_0.2.1 statmod_1.4.35
## [5] stats4_4.0.4 highr_0.8 cellranger_1.1.0 pillar_1.5.0
## [9] backports_1.2.1 glue_1.4.2 ggsignif_0.6.1 minqa_1.2.4
## [13] sandwich_3.0-0 pkgconfig_2.0.3 broom_0.7.5 haven_2.3.1
## [17] purrr_0.3.4 scales_1.1.1 openxlsx_4.2.3 rio_0.5.16
## [21] tibble_3.0.6 mgcv_1.8-33 generics_0.1.0 ellipsis_0.3.1
## [25] withr_2.4.1 magrittr_2.0.1 crayon_1.4.1 readxl_1.3.1
## [29] evaluate_0.14 fansi_0.4.2 nlme_3.1-152 rstatix_0.7.0
## [33] forcats_0.5.1 foreign_0.8-81 tools_4.0.4 data.table_1.14.0
## [37] hms_1.0.0 lifecycle_1.0.0 stringr_1.4.0 munsell_0.5.0
## [41] cluster_2.1.0 zip_2.1.1 compiler_4.0.4 rlang_0.4.10
## [45] grid_4.0.4 nloptr_1.2.2.2 boot_1.3-26 gtable_0.3.0
## [49] codetools_0.2-18 abind_1.4-5 DBI_1.1.1 curl_4.3
## [53] R6_2.5.0 zoo_1.8-8 dplyr_1.0.4 utf8_1.1.4
## [57] stringi_1.5.3 parallel_4.0.4 Rcpp_1.0.6 vctrs_0.3.6
## [61] tidyselect_1.1.0 xfun_0.21
Graphics options
TIFF <- interactive()
trellis.par.set(list(plot.symbol = list(col=1,pch=20, cex=0.7),
box.rectangle = list(col=1),
box.umbrella = list(lty=1, col=1),
strip.background = list(col = "lightgrey"),
layout.widths = list(right.padding = 5)))
ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme
ltheme$strip.background$col <- "lightgrey" ## change strip bg
lattice.options(default.theme = ltheme)
options(width = 100)
Initialise random seed
set.seed(29)
Data overview
# "PlotID" : unique identifier for sample site
# "logID" : unique identifier for wood sample
# "treatment" : experimental treatments: uncaged, open cage or closed cage
# "clade" : host tree division of wood sample: angiosperm or gymnosperm
# "dowel" : 1 = standardized machined dowel of Fagus sylvatica without bark;
# 0 = unprocessed wood of native tree species with bark
# "drymass_start": dry mass at start of experiment before exposure
# "expo_days" : exposure time in days
# "drymass_end" : dry mass at the end of experiment
# "marks_total" : signs of insect colonization (not available for all logs: NA)
# "biome" : classification of biome as either boreal, temperate or sub/tropical
# "WC2_T" : mean annual temperature according to worldclim
# "WC2_P" : mean annual precipitation according to worldclim
Rescale variables
### rescale response variables: dry mass at beginning or end of experiment, respectively
wm$drymass_end <- wm$drymass_end / 1000
wm$drymass_start <- wm$drymass_start / 1000
### rescale predictor variables: exposure time, temperature and precipitation
wm$expoyrs <- with(wm, expo_days / 365) ### in years
wm$Temp <- with(wm, WC2_T) ### in ° Celsius
### FIXME
wm$Prec <- with(wm, WC2_P / 100) ### dm per year
## scale and center in addition
wmplot <- subset(wm, !duplicated(PlotID))
## mean temperature in °C (rounded allowing direct textual reference)
(mT <- round(mean(wmplot$Temp)))
## [1] 15
## mean precipitation in dm / a (rounded allowing direct textual reference)
(mP <- round(mean(wmplot$Prec)))
## [1] 13
wm$TempS <- scale(wm$Temp, center = mT, scale = TRUE)
wm$PrecS <- scale(wm$Prec, center = mP, scale = TRUE)
### center only to preserve original measurement scale
wm$TempC <- scale(wm$Temp, center = mT, scale = FALSE)
wm$PrecC <- scale(wm$Prec, center = mP, scale = FALSE)
Main model presented in Extended Data Table 1 comparing decomposition rates of NATURAL WOOD between treatments uncaged and closed
### Subsetting data to focal treatments and either natural tree species or dowels
d_natural_uncaged_closed <- subset(wm, treatment %in% c("uncaged", "closed") &
dowel == 0)
d_dowel_uncaged_closed <- subset(wm, treatment %in% c("uncaged", "closed") &
dowel == 1)
### Model for natural wood: uncaged vs. closed (for angio / gymno)
m_n_uc <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS * clade) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_natural_uncaged_closed,
family = gaussian(link = "log"),
offset = log(drymass_start))
### increase iterations => model converges
ss <- getME(m_n_uc, c("theta","fixef"))
m_n_uc <- update(m_n_uc, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_n_uc)$optinfo$conv$opt
## [1] 0
Variance inflation factors
vif(m_n_uc)
## Warning in vif.merMod(m_n_uc): No intercept: vifs may not be sensible.
## GVIF Df GVIF^(1/(2*Df))
## expoyrs:TempS 1.207867 1 1.099030
## expoyrs:PrecS 1.312078 1 1.145460
## expoyrs:clade 2.653399 2 1.276294
## expoyrs:treatment 1.216151 1 1.102792
## expoyrs:TempS:PrecS 1.289184 1 1.135422
## expoyrs:TempS:clade 2.708672 1 1.645804
## expoyrs:PrecS:clade 6.615743 1 2.572109
## expoyrs:TempS:treatment 1.301636 1 1.140893
## expoyrs:PrecS:treatment 1.339201 1 1.157239
## expoyrs:TempS:PrecS:clade 7.322192 1 2.705955
## expoyrs:TempS:PrecS:treatment 1.433673 1 1.197361
Xscale <- model.matrix(m_n_uc)
# check correlation among predictors
cov2cor(var(Xscale))
## expoyrs:TempS expoyrs:PrecS expoyrs:cladeangio
## expoyrs:TempS 1.00000000 0.31263930 0.24577013
## expoyrs:PrecS 0.31263930 1.00000000 0.15209747
## expoyrs:cladeangio 0.24577013 0.15209747 1.00000000
## expoyrs:cladegymno -0.47692305 -0.20930517 -0.69528583
## expoyrs:treatmentuncaged -0.06283165 0.00332664 0.17631118
## expoyrs:TempS:PrecS -0.04161171 -0.14810179 -0.04741402
## expoyrs:TempS:cladegymno 0.58095059 0.21059745 0.55180631
## expoyrs:PrecS:cladegymno 0.23630619 0.42618509 0.22538700
## expoyrs:TempS:treatmentuncaged 0.69314498 0.21522828 0.16222674
## expoyrs:PrecS:treatmentuncaged 0.22360090 0.68246332 0.10351045
## expoyrs:TempS:PrecS:cladegymno -0.29440825 -0.39137224 -0.23840070
## expoyrs:TempS:PrecS:treatmentuncaged 0.01830656 -0.09329911 -0.02504622
## expoyrs:cladegymno expoyrs:treatmentuncaged
## expoyrs:TempS -0.47692305 -0.06283165
## expoyrs:PrecS -0.20930517 0.00332664
## expoyrs:cladeangio -0.69528583 0.17631118
## expoyrs:cladegymno 1.00000000 0.06825995
## expoyrs:treatmentuncaged 0.06825995 1.00000000
## expoyrs:TempS:PrecS 0.10648430 0.01919709
## expoyrs:TempS:cladegymno -0.80421381 -0.05037637
## expoyrs:PrecS:cladegymno -0.37503605 -0.03084315
## expoyrs:TempS:treatmentuncaged -0.32467203 -0.19498929
## expoyrs:PrecS:treatmentuncaged -0.14624837 0.02918454
## expoyrs:TempS:PrecS:cladegymno 0.37206529 0.01953970
## expoyrs:TempS:PrecS:treatmentuncaged 0.05422587 0.20823567
## expoyrs:TempS:PrecS expoyrs:TempS:cladegymno
## expoyrs:TempS -0.04161171 0.58095059
## expoyrs:PrecS -0.14810179 0.21059745
## expoyrs:cladeangio -0.04741402 0.55180631
## expoyrs:cladegymno 0.10648430 -0.80421381
## expoyrs:treatmentuncaged 0.01919709 -0.05037637
## expoyrs:TempS:PrecS 1.00000000 -0.20383802
## expoyrs:TempS:cladegymno -0.20383802 1.00000000
## expoyrs:PrecS:cladegymno -0.47516483 0.40669325
## expoyrs:TempS:treatmentuncaged 0.01518715 0.39172225
## expoyrs:PrecS:treatmentuncaged -0.10257555 0.14228505
## expoyrs:TempS:PrecS:cladegymno 0.52223245 -0.50961597
## expoyrs:TempS:PrecS:treatmentuncaged 0.69057075 -0.11543984
## expoyrs:PrecS:cladegymno expoyrs:TempS:treatmentuncaged
## expoyrs:TempS 0.23630619 0.6931449757
## expoyrs:PrecS 0.42618509 0.2152282782
## expoyrs:cladeangio 0.22538700 0.1622267419
## expoyrs:cladegymno -0.37503605 -0.3246720295
## expoyrs:treatmentuncaged -0.03084315 -0.1949892937
## expoyrs:TempS:PrecS -0.47516483 0.0151871494
## expoyrs:TempS:cladegymno 0.40669325 0.3917222535
## expoyrs:PrecS:cladegymno 1.00000000 0.1483911729
## expoyrs:TempS:treatmentuncaged 0.14839117 1.0000000000
## expoyrs:PrecS:treatmentuncaged 0.29995452 0.3121543640
## expoyrs:TempS:PrecS:cladegymno -0.91207557 -0.1845123855
## expoyrs:TempS:PrecS:treatmentuncaged -0.30713852 0.0005987126
## expoyrs:PrecS:treatmentuncaged expoyrs:TempS:PrecS:cladegymno
## expoyrs:TempS 0.22360090 -0.2944083
## expoyrs:PrecS 0.68246332 -0.3913722
## expoyrs:cladeangio 0.10351045 -0.2384007
## expoyrs:cladegymno -0.14624837 0.3720653
## expoyrs:treatmentuncaged 0.02918454 0.0195397
## expoyrs:TempS:PrecS -0.10257555 0.5222324
## expoyrs:TempS:cladegymno 0.14228505 -0.5096160
## expoyrs:PrecS:cladegymno 0.29995452 -0.9120756
## expoyrs:TempS:treatmentuncaged 0.31215436 -0.1845124
## expoyrs:PrecS:treatmentuncaged 1.00000000 -0.2711983
## expoyrs:TempS:PrecS:cladegymno -0.27119828 1.0000000
## expoyrs:TempS:PrecS:treatmentuncaged -0.13061873 0.3342147
## expoyrs:TempS:PrecS:treatmentuncaged
## expoyrs:TempS 0.0183065587
## expoyrs:PrecS -0.0932991076
## expoyrs:cladeangio -0.0250462188
## expoyrs:cladegymno 0.0542258673
## expoyrs:treatmentuncaged 0.2082356703
## expoyrs:TempS:PrecS 0.6905707533
## expoyrs:TempS:cladegymno -0.1154398423
## expoyrs:PrecS:cladegymno -0.3071385207
## expoyrs:TempS:treatmentuncaged 0.0005987126
## expoyrs:PrecS:treatmentuncaged -0.1306187326
## expoyrs:TempS:PrecS:cladegymno 0.3342147097
## expoyrs:TempS:PrecS:treatmentuncaged 1.0000000000
R2 : marginal (integrating out the random effects) and conditional (given the random effects)
r.squaredGLMM(m_n_uc)
## R2m R2c
## [1,] 0.8394324 0.9773549
Model fit: Plot observed dry mass vs fitted dry mass
d_natural_uncaged_closed$fitted <- fitted(m_n_uc)
vcol <- viridis_pal()(3)
tmp <- d_natural_uncaged_closed
levels(tmp$clade) <- c("Angiosperm", "Gymnosperm")
plt <- xyplot(fitted ~ drymass_end | clade + treatment, data = tmp,
groups = Biome, pch = 19, col = vcol,
par.settings = list(superpose.symbol = list (col=vcol, pch = 19)),
between = list(x = 1, y = 1),
xlab = "Observed dry mass", ylab = "Fitted dry mass",
auto.key = list(space = "right", title = "Biome",
lines = FALSE, points = TRUE,
col = "black"))
if (TIFF)
tiff(filename = "products/ObsFitted.tif",width = 7, height = 7.5,res=300,units="in",compression = "lzw")
plot(plt)
if (TIFF)
dev.off()
Model output with transformed variables for natural wood comparing uncaged vs. closed
summary(m_n_uc)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS * clade) +
## (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_natural_uncaged_closed
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -8427.6 -8345.9 4227.8 -8455.6 2519
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.9672 -0.3395 0.0793 0.4710 10.1158
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.003108 0.05575
## Residual 0.002095 0.04577
## Number of obs: 2533, groups: PlotID, 55
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.087044 0.023889 -3.644 0.000269 ***
## expoyrs:PrecS -0.021232 0.022495 -0.944 0.345244
## expoyrs:cladeangio -0.150504 0.022506 -6.687 2.27e-11 ***
## expoyrs:cladegymno -0.082867 0.024862 -3.333 0.000859 ***
## expoyrs:treatmentuncaged -0.029228 0.005694 -5.133 2.85e-07 ***
## expoyrs:TempS:PrecS -0.030245 0.021493 -1.407 0.159359
## expoyrs:TempS:cladegymno 0.039649 0.009882 4.012 6.02e-05 ***
## expoyrs:PrecS:cladegymno -0.002946 0.024293 -0.121 0.903473
## expoyrs:TempS:treatmentuncaged -0.032905 0.005867 -5.608 2.04e-08 ***
## expoyrs:PrecS:treatmentuncaged -0.035457 0.006249 -5.674 1.40e-08 ***
## expoyrs:TempS:PrecS:cladegymno 0.005545 0.017512 0.317 0.751508
## expoyrs:TempS:PrecS:treatmentuncaged -0.038991 0.006044 -6.451 1.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:cldn expyrs:cldg expyrs:t ex:TS:PS expyrs:TmpS:c
## expyrs:PrcS -0.287
## expyrs:cldn -0.190 -0.028
## expyrs:cldg -0.154 -0.039 0.885
## expyrs:trtm -0.036 0.021 -0.121 -0.129
## expyr:TS:PS -0.019 -0.324 -0.255 -0.218 0.030
## expyrs:TmpS:c -0.017 0.007 0.010 0.315 -0.030 -0.006
## expyrs:PrcS:c 0.032 -0.055 -0.025 0.113 -0.001 0.038 0.265
## expyrs:TmpS:t -0.117 0.039 -0.034 -0.042 0.324 -0.020 -0.034
## expyrs:PrcS:t 0.022 -0.103 0.021 0.023 -0.087 -0.040 0.009
## expyrs:TmpS:PrcS:c 0.017 -0.016 -0.008 0.106 -0.003 -0.002 0.424
## expyrs:TmpS:PrcS:t -0.022 -0.026 0.035 0.035 -0.172 -0.118 0.013
## expyrs:PrcS:c expyrs:TmpS:t expyrs:PrcS:t expyrs:TmpS:PrcS:c
## expyrs:PrcS
## expyrs:cldn
## expyrs:cldg
## expyrs:trtm
## expyr:TS:PS
## expyrs:TmpS:c
## expyrs:PrcS:c
## expyrs:TmpS:t 0.008
## expyrs:PrcS:t -0.009 -0.157
## expyrs:TmpS:PrcS:c 0.889 0.014 -0.003
## expyrs:TmpS:PrcS:t -0.004 0.168 0.428 -0.016
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00285545 (tol = 0.002, component 1)
Transform temperature and precipitation back to degree Celsius and dm per year.
TempC is temperature in °C - 15 and PrecC is precipitation in dm / a - 13.
Effects of TempC can be interpreted as effect of increasing temperature by 1 °C when precipitation is 13 dm/a.
Effects of PreC can be interpreted as effect an additional dm / a when temperature is 15 °C.
Xscale <- model.matrix(m_n_uc)
Xorig <- model.matrix(~ 0 + ### intercept zero
expoyrs:(TempC * PrecC * clade) + ### slopes varying in temp and prec
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempC * PrecC),
data = droplevels(d_natural_uncaged_closed))
K <- solve(crossprod(Xorig)) %*% crossprod(Xorig, Xscale)
g <- glht(m_n_uc, linfct = K)
## check back-transformation
a <- c(predict(m_n_uc, newdata = d_natural_uncaged_closed, re.form = NA))
stopifnot(all.equal(coef(g), coef(lm(a ~ 0 + Xorig)), check.attributes = FALSE))
stopifnot(all.equal(c(Xorig %*% coef(g)), a, check.attributes = FALSE))
Create main model output table (Extended Data Table 1) for natural wood
### Model output
(tab <- summary(g, test = univariate()))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS * clade) + (0 + expoyrs | PlotID) + expoyrs:treatment +
## expoyrs:treatment:(TempS * PrecS), data = d_natural_uncaged_closed,
## family = gaussian(link = "log"), control = glmerControl(optCtrl = list(maxfun = 2e+06)),
## start = ss, offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:TempC == 0 -0.0110069 0.0030209 -3.644 0.000269 ***
## expoyrs:PrecC == 0 -0.0031350 0.0033215 -0.944 0.345244
## expoyrs:cladeangio == 0 -0.1505040 0.0225062 -6.687 2.27e-11 ***
## expoyrs:cladegymno == 0 -0.0828672 0.0248622 -3.333 0.000859 ***
## expoyrs:treatmentuncaged == 0 -0.0292280 0.0056937 -5.133 2.85e-07 ***
## expoyrs:TempC:PrecC == 0 -0.0005647 0.0004013 -1.407 0.159359
## expoyrs:TempC:cladegymno == 0 0.0050137 0.0012496 4.012 6.02e-05 ***
## expoyrs:PrecC:cladegymno == 0 -0.0004350 0.0035868 -0.121 0.903473
## expoyrs:TempC:treatmentuncaged == 0 -0.0041608 0.0007419 -5.608 2.04e-08 ***
## expoyrs:PrecC:treatmentuncaged == 0 -0.0052353 0.0009227 -5.674 1.40e-08 ***
## expoyrs:TempC:PrecC:cladegymno == 0 0.0001035 0.0003270 0.317 0.751508
## expoyrs:TempC:PrecC:treatmentuncaged == 0 -0.0007280 0.0001128 -6.451 1.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Relative effects and confidence intervals
### confidence intervals
ci <- confint(g, calpha = univariate_calpha())$confint
### Relative effects
### Parameter interpretation
### "expoyrs:TempC": multiplicative effect of 1 °C temperature increase on
### annual dry mass _when_ precipitation is 13 dm/a
### 0.989 means that when we increase the temperature by 1 °C and keep all
### other variables fix (precipitation at 13 dm/a, clade and treatment), after one
### year we observe 98.9% of the dry mass we would have observed without
### this intervention. This is equivalent to 1.1% additional mass loss
### induced by the 1 °C higher temperature.
###
### "expoyrs:PrecC": multiplicative effect of 1 am/a precipitation increase on
### annual dry mass _when_ temperature is 13 °C
### 0.997 means that when we increase precipitation by 1 dm/a and keep all
### other variables fix (temperature at 13 °C, clade and treatment), after one
### year we observe 99.7% of the dry mass we would have observed without
### this intervention.
###
### "expoyrs:cladeangio": A multiplicative constant for angiosperm species
### "expoyrs:cladegymno": A multiplicative constant for gymnosperm species
###
### "expoyrs:treatmentuncaged": Multiplicative treatment effect comparing "uncaged" to the
### "closed" baseline
### 0.971 means that with insect access to the wood, after one year
### we observe 97.1% of the dry mass we would have observed without insects
### at a temperature of 15 °C and precipitation equal to 13 dm/a.
###
### "expoyrs:TempC:PrecC": Multiplicative interaction effect of temperature
### and precipitation.
### Only present when temperature != 15 °C and precipitation != 13 dm/a.
### If clade == "angio" and treatment == "closed", 0.999 means that a one
### unit increase in both temperature and precipitation leads to
### 0.989 x 0.997 x 0.999 x 100 % the annual dry mass we would have observed
### without such an intervention (product of main and interaction effects).
###
### "expoyrs:TempC:cladegymno": The additional temperature effect for gymnosperm species
### "expoyrs:PrecC:cladegymno": The additional precipitation effect for gymnosperm species
###
### "expoyrs:TempC:treatmentuncaged": Temperature as treatment effect
### modifier when precipitation is 13 dm/a
### 0.996 means that a temperature increase by 1 °C with insect access leads to 99.6%
### of the dry mass we would have observed (after one year) with insect
### access at a precipitation of 13 dm /a.
###
### "expoyrs:PrecC:treatmentuncaged": Precipitation as treatment effect
### modifier when temperature is 13 °C
### 0.995 means that 1 dm/a more precipitation with insect access leads to 99.6%
### of the dry mass we would have observed (after one year) with insect
### access at 13 °C
###
### "expoyrs:TempC:PrecC:cladegymno": The additional temperature x precipitation
### interaction effect for gymnosperm species
###
### "expoyrs:TempC:PrecC:treatmentuncaged": Temperature and precipitation as treatment effect
### modifiers
### 0.999 means that a one unit increase in both temperature and precipitation leads to
### 0.989 x 0.997 x 0.971 x 0.999 x 0.996 x 0.995 x 0.999 x 100% the annual dry mass we would have observed
### without such an intervention (product of main and interaction effects)
### for angiosperm species.
###
exp(ci)
## Estimate lwr upr
## expoyrs:TempC 0.9890535 0.9832148 0.9949268
## expoyrs:PrecC 0.9968699 0.9904015 1.0033807
## expoyrs:cladeangio 0.8602743 0.8231512 0.8990716
## expoyrs:cladegymno 0.9204734 0.8766948 0.9664380
## expoyrs:treatmentuncaged 0.9711950 0.9604172 0.9820937
## expoyrs:TempC:PrecC 0.9994355 0.9986497 1.0002218
## expoyrs:TempC:cladegymno 1.0050263 1.0025677 1.0074908
## expoyrs:PrecC:cladegymno 0.9995651 0.9925627 1.0066169
## expoyrs:TempC:treatmentuncaged 0.9958478 0.9944008 0.9972969
## expoyrs:PrecC:treatmentuncaged 0.9947784 0.9929810 0.9965790
## expoyrs:TempC:PrecC:cladegymno 1.0001035 0.9994628 1.0007447
## expoyrs:TempC:PrecC:treatmentuncaged 0.9992723 0.9990513 0.9994933
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
Visualisation of treatment effect in the presense of interaction effects. For a given temperature, precipitation, and host (angiosperm or gymnosperm), the expected wood mass loss after one year is “L” under closed cages. Uncaged, the wood mass loss changes to “p * L”, is p x 100% of the wood mass loss under closed cases. The relative treatment effect p is plotted here.
Interpretation: For sites with higher temperature (~ > 4 °C), increasing precipitation leads to increased dry mass loss in the uncaged case. In colder environments, increasing precipitation reduces dry mass loss and, in extremely cold places (-2 °C), precipitation of more than ~ 10 dm/a wood mass increases over the course of a year up to ~6 %.
qtm <- -.5:7 * 4
tmi <- cut(d_natural_uncaged_closed$Temp, breaks = qtm)
pr <- seq(from = min(d_natural_uncaged_closed$Prec),
to = max(d_natural_uncaged_closed$Prec), length = 20)
nd <- expand.grid(Temp = qtm, Prec = pr,
clade = unique(d_natural_uncaged_closed$clade),
treatment = unique(d_natural_uncaged_closed$treatment)[, drop = TRUE],
expoyrs = 1)
nd$Tempf <- factor(nd$Temp, levels = qtm,
labels = paste(ifelse(qtm > 0 & qtm < 10, " ", ""), qtm))
### compute convex hull of observed Temp and Prec and plot
### model only inside this convex hull
angio <- subset(subset(wm, clade == "angio"), !duplicated(PlotID))
gymno <- subset(subset(wm, clade == "gymno"), !duplicated(PlotID))
ch_a <- tri.mesh(angio$Temp, angio$Prec)
ch_g <- tri.mesh(gymno$Temp, gymno$Prec)
nd$inhull <- TRUE
nd$inhull[nd$clade == "angio"] <- in.convex.hull(ch_a,
nd[nd$clade == "angio", "Temp"], nd[nd$clade == "angio", "Prec"])
nd$inhull[nd$clade == "gymno"] <- in.convex.hull(ch_g,
nd[nd$clade == "gymno", "Temp"], nd[nd$clade == "gymno", "Prec"])
nd_1 <- subset(nd, treatment == "uncaged")
nd_2 <- subset(nd, treatment == "closed")
nd_1$TempS <- scale(nd_1$Temp, center = attr(wm$TempS, "scaled:center"),
scale = attr(wm$TempS, "scaled:scale"))
nd_2$TempS <- scale(nd_2$Temp, center = attr(wm$TempS, "scaled:center"),
scale = attr(wm$TempS, "scaled:scale"))
nd_1$PrecS <- scale(nd_1$Prec, center = attr(wm$PrecS, "scaled:center"),
scale = attr(wm$PrecS, "scaled:scale"))
nd_2$PrecS <- scale(nd_2$Prec, center = attr(wm$PrecS, "scaled:center"),
scale = attr(wm$PrecS, "scaled:scale"))
nd_1$lp <- predict(m_n_uc, newdata = nd_1, re.form = NA) -
predict(m_n_uc, newdata = nd_2, re.form = NA)
col <- rev(sequential_hcl(length(qtm), palette = "Red-Blue"))
lwd <- 2
nd_1$lp_u <- predict(m_n_uc, newdata = nd_1, re.form = NA)
levels(nd_1$clade) <- c("Angiosperm", "Gymnosperm")
# Interpretation: For sites with higher temperature (~ > 4 °C), increasing
# precipitation leads to increased dry mass loss in the uncaged case. In
# colder environments, increasing precipitation reduces dry mass loss and,
# in extremely cold places (-2 °C), precipitation of more than ~ 10 dm/a
# wood mass increases over the course of a year up to ~6 %.
p1 <- xyplot(I((1 - exp(lp_u)) * 100) ~ Prec | clade, data = nd_1, groups = Tempf, type = "l",
col = col, xlab = "Precipitation [dm / a]",
ylab = "Annual mass loss [%]", lwd = lwd,
scales = list(tck = c(1, 0), alternating = FALSE), subset = inhull,
between = list(x = 1),
panel = function(...) {
panel.abline(h = 0, col = "lightgrey", lty = 2)
panel.xyplot(...)
},
main = "a All decomposers",
par.settings = list(superpose.line = list (col=col, lty = 1, lwd = lwd),
par.main.text = list(font = 1,
just = "left",
x = grid::unit(5, "mm"))),
auto.key = list(space = "right", title = "Temp. [°C]",
columns = 1, lines = TRUE, points = FALSE,
col = "black"))
# Visualisation of treatment effect.
#
# Interpretation: When comparing uncaged vs. closed treatments,
# insect access ("uncaged") at sites with higher temperature (~ > 4 °C) leads to
# accelerated dry mass loss, the effect increases with increasing
# precipitation. In colder environments, insect access leads to a stronger
# accumulation of dry mass over the course of a year than observed without
# insect access.
p2 <- xyplot(I((1 - exp(lp)) * 100) ~ Prec | clade, data = nd_1, groups = Tempf, type = "l",
col = col, xlab = "Precipitation [dm / a]",
ylab = "Annual mass loss [%]", lwd = lwd,
scales = list(tck = c(1, 0), alternating = FALSE), subset = inhull,
between = list(x = 1),
panel = function(...) {
panel.abline(h = 0, col = "lightgrey", lty = 2)
panel.xyplot(...)
},
main = "b Net insect effect",
par.settings = list(superpose.line = list (col=col, lty = 1, lwd = lwd),
par.main.text = list(font = 1,
just = "left",
x = grid::unit(5, "mm"))),
auto.key = list(space = "right", title = "Temp. [°C]",
columns = 1, lines = TRUE, points = FALSE,
col = "black"))
if (TIFF)
tiff(filename = "products/ExtendedDataFigure4.tif",width = 7, height = 7.5,res=300,units="in",compression = "lzw")
grid.arrange(arrangeGrob(p1, p2, nrow = 2))
if (TIFF)
dev.off()
Second model presented in Table 1 comparing decomposition rates of DOWELS between treatments uncaged and closed
m_d_uc <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_dowel_uncaged_closed,
family = gaussian(link = "log"),
offset = log(drymass_start))
ss <- getME(m_d_uc, c("theta","fixef"))
m_d_uc <- update(m_d_uc, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_d_uc)$optinfo$conv$opt
## [1] 0
Model output with transformed variables for dowels comparing uncaged vs. closed
summary(m_d_uc)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS) + (0 + expoyrs |
## PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_dowel_uncaged_closed
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -1514.1 -1473.8 767.0 -1534.1 405
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2170 -0.3523 0.0901 0.4333 3.2177
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.00619 0.07868
## Residual 0.00138 0.03715
## Number of obs: 415, groups: PlotID, 54
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.14170 0.03102 -4.569 4.91e-06 ***
## expoyrs:PrecS -0.04196 0.02891 -1.451 0.14670
## expoyrs:treatmentclosed -0.11977 0.02836 -4.223 2.41e-05 ***
## expoyrs:treatmentuncaged -0.16634 0.02775 -5.995 2.04e-09 ***
## expoyrs:TempS:PrecS -0.08554 0.02688 -3.182 0.00146 **
## expoyrs:TempS:treatmentuncaged -0.03973 0.01263 -3.144 0.00166 **
## expoyrs:PrecS:treatmentuncaged -0.02435 0.01317 -1.849 0.06451 .
## expoyrs:TempS:PrecS:treatmentuncaged -0.03412 0.01471 -2.320 0.02035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:trtmntc expyrs:trtmntn ex:TS:PS ex:TS: ex:PS:
## expyrs:PrcS -0.392
## expyrs:trtmntc -0.109 -0.016
## expyrs:trtmntn -0.168 -0.013 0.917
## expyr:TS:PS 0.148 -0.219 -0.291 -0.302
## expyrs:TmS: -0.274 0.009 -0.122 0.084 -0.108
## expyrs:PrS: -0.035 -0.218 -0.013 0.023 -0.195 0.146
## expy:TS:PS: -0.093 -0.107 -0.009 0.022 -0.325 0.384 0.654
Transform temperature and precipitation back to degree Celsius and dm per year
Xscale <- model.matrix(m_d_uc)
Xorig <- model.matrix(~ 0 + ### intercept zero
expoyrs:(TempC * PrecC) + ### slopes varying in temp and prec
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempC * PrecC),
data = droplevels(d_dowel_uncaged_closed))
K <- solve(crossprod(Xorig)) %*% crossprod(Xorig, Xscale)
g <- glht(m_d_uc, linfct = K)
Create main model output table (Table 1) for dowels
### Model output
(tab <- summary(g, test = univariate()))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS) + (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS *
## PrecS), data = d_dowel_uncaged_closed, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:TempC == 0 -0.0179181 0.0039220 -4.569 4.91e-06 ***
## expoyrs:PrecC == 0 -0.0061949 0.0042685 -1.451 0.14670
## expoyrs:treatmentclosed == 0 -0.1197672 0.0283608 -4.223 2.41e-05 ***
## expoyrs:treatmentuncaged == 0 -0.1663416 0.0277477 -5.995 2.04e-09 ***
## expoyrs:TempC:PrecC == 0 -0.0015970 0.0005019 -3.182 0.00146 **
## expoyrs:TempC:treatmentuncaged == 0 -0.0050237 0.0015976 -3.144 0.00166 **
## expoyrs:PrecC:treatmentuncaged == 0 -0.0035949 0.0019446 -1.849 0.06451 .
## expoyrs:TempC:PrecC:treatmentuncaged == 0 -0.0006371 0.0002746 -2.320 0.02035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Relative effects and confidence intervals
### confidence intervals
ci <- confint(g, calpha = univariate_calpha())$confint
### relative effects
exp(ci)
## Estimate lwr upr
## expoyrs:TempC 0.9822415 0.9747199 0.9898211
## expoyrs:PrecC 0.9938243 0.9855445 1.0021736
## expoyrs:treatmentclosed 0.8871270 0.8391605 0.9378352
## expoyrs:treatmentuncaged 0.8467569 0.8019363 0.8940826
## expoyrs:TempC:PrecC 0.9984042 0.9974225 0.9993869
## expoyrs:TempC:treatmentuncaged 0.9949889 0.9918781 0.9981094
## expoyrs:PrecC:treatmentuncaged 0.9964115 0.9926211 1.0002165
## expoyrs:TempC:PrecC:treatmentuncaged 0.9993632 0.9988254 0.9999012
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
### Relative effects for comparison between treatments
### Treatment effect when temperature is 15 ° Celsius and
### when precipitation is 13 dm/a
g <- glht(m_d_uc, linfct = c("expoyrs:treatmentuncaged - expoyrs:treatmentclosed = 0"))
summary(g)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS) + (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS *
## PrecS), data = d_dowel_uncaged_closed, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:treatmentuncaged - expoyrs:treatmentclosed == 0 -0.04657 0.01145 -4.067 4.76e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(g)$confint)
## Estimate lwr upr
## expoyrs:treatmentuncaged - expoyrs:treatmentclosed 0.9544935 0.9333089 0.976159
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
Compute fitted values for annual wood decomposition rates in uncaged and closed treatment for each experimental site, seperately for natural wood (angiosperm and gymnosperm seperately where applicable) and dowels
### For natural wood, fitted values include three scenarios: main result based on model estimates used to generate Table 1
### and "Model uncertainity scenarios" based on upper and lower confidence intervals of fixed effects used for sensitivity analysis (Extended Data Table 2)
a <- subset(wm, treatment == "uncaged" & clade == "angio")
TP_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "Temp", "Prec")]
}))
g <- subset(wm, clade != "angio")
TP_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "Temp", "Prec")]
}))
TPS_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_a$expoyrs <- 1
clade <- sort(unique(wm$clade))
TPS_a$clade <- clade[1]
TPS_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_g$expoyrs <- 1
TPS_g$clade <- clade[2]
trt <- sort(unique(wm$treatment))
TPS_a$treatment <- trt[trt == "uncaged"]
TPS_a$rate_n_uncaged <- exp(predict(m_n_uc, newdata = TPS_a))
fm <- ~ 0 + expoyrs:(TempS * PrecS * clade) + ### slopes varying in temp and prec
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS)
K_a <- model.matrix(fm, data = TPS_a)[, names(fixef(m_n_uc))]
ci <- confint(glht(m_n_uc, linfct = K_a),
calpha = univariate_calpha())$confint
colnames(ci) <- paste0("rate_n_uncaged", colnames(ci))
TPS_a <- cbind(TPS_a,
exp(
### confint for FIXED effects
ci +
### add estimated RANDOM effect (w/o uncertainty!!!)
predict(m_n_uc, newdata = TPS_a, random.only = TRUE)
))
stopifnot(with(TPS_a, all.equal(rate_n_uncaged, rate_n_uncagedEstimate)))
TPS_a$treatment <- trt[trt == "closed"]
TPS_a$rate_n_closed <- exp(predict(m_n_uc, newdata = TPS_a))
K_a <- model.matrix(fm, data = TPS_a)[, names(fixef(m_n_uc))]
ci <- confint(glht(m_n_uc, linfct = K_a),
calpha = univariate_calpha())$confint
colnames(ci) <- paste0("rate_n_closed", colnames(ci))
TPS_a <- cbind(TPS_a,
exp(
### confint for FIXED effects
ci +
### add estimated RANDOM effect (w/o uncertainty!!!)
predict(m_n_uc, newdata = TPS_a, random.only = TRUE)
))
stopifnot(with(TPS_a, all.equal(rate_n_closed, rate_n_closedEstimate)))
TPS_g$treatment <- trt[trt == "uncaged"]
TPS_g$rate_n_uncaged <- exp(predict(m_n_uc, newdata = TPS_g))
K_g <- model.matrix(fm, data = TPS_g)[, names(fixef(m_n_uc))]
ci <- confint(glht(m_n_uc, linfct = K_g),
calpha = univariate_calpha())$confint
colnames(ci) <- paste0("rate_n_uncaged", colnames(ci))
TPS_g <- cbind(TPS_g,
exp(
### confint for FIXED effects
ci +
### add estimated RANDOM effect (w/o uncertainty!!!)
predict(m_n_uc, newdata = TPS_g, random.only = TRUE)
))
stopifnot(with(TPS_g, all.equal(rate_n_uncaged, rate_n_uncagedEstimate)))
TPS_g$treatment <- trt[trt == "closed"]
TPS_g$rate_n_closed <- exp(predict(m_n_uc, newdata = TPS_g))
K_g <- model.matrix(fm, data = TPS_g)[, names(fixef(m_n_uc))]
ci <- confint(glht(m_n_uc, linfct = K_g),
calpha = univariate_calpha())$confint
colnames(ci) <- paste0("rate_n_closed", colnames(ci))
TPS_g <- cbind(TPS_g,
exp(
### confint for FIXED effects
ci +
### add estimated RANDOM effect (w/o uncertainty!!!)
predict(m_n_uc, newdata = TPS_g, random.only = TRUE)
))
stopifnot(with(TPS_g, all.equal(rate_n_closed, rate_n_closedEstimate)))
biome <- do.call("rbind", tapply(1:nrow(wm), wm$PlotID, function(i)
wm[i[1], c("PlotID", "biome")]))
x_a <- merge(merge(TP_a, TPS_a, by = "PlotID"), biome, by = "PlotID")
x_g <- merge(merge(TP_g, TPS_g, by = "PlotID"), biome, by = "PlotID")
### Results for natural wood for uncaged and closed treatments
ret_n <- rbind(x_a, x_g)
no_dowel <- levels(wm$PlotID)[table(d_dowel_uncaged_closed$PlotID) == 0]
tmp <- subset(TPS_a, !(PlotID %in% no_dowel))
tmp$rate_n_uncaged <- tmp$rate_n_closed <- tmp$rate_n_uncagedEstimate <- tmp$rate_n_closedEstimate <- tmp$rate_n_uncagedlwr <- tmp$rate_n_closedlwr <- tmp$rate_n_uncagedupr <- tmp$rate_n_closedupr<- NULL
tmp$treatment <- trt[trt == "uncaged"]
tmp$rate_d_uncaged <- exp(predict(m_d_uc, newdata = tmp))
tmp$treatment <- trt[trt == "closed"]
tmp$rate_d_closed <- exp(predict(m_d_uc, newdata = tmp))
### Results for dowels for uncaged and closed treatments
ret_d <- merge(merge(TP_a, tmp, by = "PlotID"), biome, by = "PlotID")
Additional models comparing decomposition rates between open and uncaged treatments (not shown in paper)
### Subsetting data to focal treatments and either natural tree species or dowels
d_natural_uncaged_open <- subset(wm, treatment %in% c("uncaged", "open") &
dowel == 0)
d_dowel_uncaged_open <- subset(wm, treatment %in% c("uncaged", "open") &
dowel == 1)
### Model for natural wood: uncaged vs. open (for angio / gymno)
m_n_uo <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS * clade) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_natural_uncaged_open,
family = gaussian(link = "log"),
offset = log(drymass_start))
ss <- getME(m_n_uo , c("theta","fixef"))
m_n_uo <- update(m_n_uo, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_n_uo)$optinfo$conv$opt
## [1] 0
summary(m_n_uo)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS * clade) +
## (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_natural_uncaged_open
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -8023.0 -7941.7 4025.5 -8051.0 2447
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0577 -0.2816 0.0732 0.4238 9.4172
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.003573 0.05977
## Residual 0.002242 0.04735
## Number of obs: 2461, groups: PlotID, 55
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.119215 0.025587 -4.659 3.18e-06 ***
## expoyrs:PrecS -0.032905 0.024467 -1.345 0.178669
## expoyrs:cladeangio -0.166628 0.024627 -6.766 1.32e-11 ***
## expoyrs:cladegymno -0.103348 0.026690 -3.872 0.000108 ***
## expoyrs:treatmentuncaged -0.020840 0.005763 -3.616 0.000299 ***
## expoyrs:TempS:PrecS -0.044277 0.023050 -1.921 0.054740 .
## expoyrs:TempS:cladegymno 0.043437 0.009651 4.501 6.77e-06 ***
## expoyrs:PrecS:cladegymno -0.014619 0.023387 -0.625 0.531928
## expoyrs:TempS:treatmentuncaged -0.014606 0.006109 -2.391 0.016801 *
## expoyrs:PrecS:treatmentuncaged -0.009870 0.006834 -1.444 0.148664
## expoyrs:TempS:PrecS:cladegymno -0.009960 0.016892 -0.590 0.555417
## expoyrs:TempS:PrecS:treatmentuncaged -0.012208 0.006617 -1.845 0.065044 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:cldn expyrs:cldg expyrs:t ex:TS:PS expyrs:TmpS:c
## expyrs:PrcS -0.324
## expyrs:cldn -0.193 -0.012
## expyrs:cldg -0.163 -0.021 0.906
## expyrs:trtm -0.040 0.019 -0.119 -0.119
## expyr:TS:PS 0.019 -0.316 -0.289 -0.257 0.020
## expyrs:TmpS:c -0.023 0.005 0.012 0.281 -0.023 -0.003
## expyrs:PrcS:c 0.027 -0.045 -0.023 0.101 0.002 0.036 0.258
## expyrs:TmpS:t -0.127 0.025 -0.035 -0.038 0.321 -0.016 0.010
## expyrs:PrcS:t 0.012 -0.121 0.017 0.018 -0.073 -0.056 -0.001
## expyrs:TmpS:PrcS:c 0.014 -0.009 -0.008 0.094 -0.002 -0.003 0.402
## expyrs:TmpS:PrcS:t -0.016 -0.041 0.021 0.020 -0.131 -0.141 0.021
## expyrs:PrcS:c expyrs:TmpS:t expyrs:PrcS:t expyrs:TmpS:PrcS:c
## expyrs:PrcS
## expyrs:cldn
## expyrs:cldg
## expyrs:trtm
## expyr:TS:PS
## expyrs:TmpS:c
## expyrs:PrcS:c
## expyrs:TmpS:t 0.014
## expyrs:PrcS:t -0.011 -0.120
## expyrs:TmpS:PrcS:c 0.879 0.032 -0.011
## expyrs:TmpS:PrcS:t -0.018 0.139 0.440 -0.022
### Model for dowels: uncaged vs. open
m_d_uo <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_dowel_uncaged_open,
family = gaussian(link = "log"),
offset = log(drymass_start))
ss <- getME(m_d_uo, c("theta","fixef"))
m_d_uo <- update(m_d_uo, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_d_uc)$optinfo$conv$opt
## [1] 0
summary(m_d_uo)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS) + (0 + expoyrs |
## PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_dowel_uncaged_open
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -2000.2 -1957.4 1010.1 -2020.2 522
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9651 -0.2613 0.0838 0.4110 4.2295
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.020557 0.14338
## Residual 0.001366 0.03696
## Number of obs: 532, groups: PlotID, 55
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.359931 0.079924 -4.503 6.69e-06 ***
## expoyrs:PrecS -0.003131 0.077391 -0.040 0.96773
## expoyrs:treatmentopen -0.228259 0.073926 -3.088 0.00202 **
## expoyrs:treatmentuncaged -0.236642 0.074002 -3.198 0.00138 **
## expoyrs:TempS:PrecS -0.160763 0.066325 -2.424 0.01536 *
## expoyrs:TempS:treatmentuncaged -0.004508 0.010497 -0.429 0.66761
## expoyrs:PrecS:treatmentuncaged -0.007555 0.011458 -0.659 0.50967
## expoyrs:TempS:PrecS:treatmentuncaged -0.013677 0.012604 -1.085 0.27787
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:trtmntp expyrs:trtmntn ex:TS:PS ex:TS: ex:PS:
## expyrs:PrcS -0.373
## expyrs:trtmntp -0.398 0.111
## expyrs:trtmntn -0.402 0.110 0.992
## expyr:TS:PS 0.091 -0.613 -0.276 -0.277
## expyrs:TmS: -0.061 -0.010 -0.027 0.044 -0.036
## expyrs:PrS: -0.025 -0.048 -0.014 0.008 -0.057 0.231
## expy:TS:PS: -0.036 -0.028 -0.013 0.008 -0.083 0.415 0.717
### Relative effects for comparison between treatments
### Treatment effect when temperature is 15 ° Celsius and
### when precipitation is 13 dm/a
g <- glht(m_d_uo, linfct = c("expoyrs:treatmentuncaged - expoyrs:treatmentopen = 0"))
summary(g)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS) + (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS *
## PrecS), data = d_dowel_uncaged_open, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:treatmentuncaged - expoyrs:treatmentopen == 0 -0.008383 0.009202 -0.911 0.362
## (Adjusted p values reported -- single-step method)
exp(confint(g)$confint)
## Estimate lwr upr
## expoyrs:treatmentuncaged - expoyrs:treatmentopen 0.9916517 0.9739272 1.009699
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
### Compute fitted values for annual wood decomposition rates in uncaged and
### open treatment for each experimental site, seperately for natural wood
### (angiosperm and gymnosperm seperately were applicable) and dowels
a <- subset(wm, treatment == "uncaged" & clade == "angio")
TP_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "Temp", "Prec")]
}))
g <- subset(wm, clade != "angio")
TP_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "Temp", "Prec")]
}))
TPS_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_a$expoyrs <- 1
clade <- sort(unique(wm$clade))
TPS_a$clade <- clade[1]
TPS_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_g$expoyrs <- 1
TPS_g$clade <- clade[2]
trt <- sort(unique(wm$treatment))
TPS_a$treatment <- trt[trt == "uncaged"]
TPS_a$rate_n_uncaged_uo <- exp(predict(m_n_uo, newdata = TPS_a))
TPS_a$treatment <- trt[trt == "open"]
TPS_a$rate_n_open_uo <- exp(predict(m_n_uo, newdata = TPS_a))
TPS_g$treatment <- trt[trt == "uncaged"]
TPS_g$rate_n_uncaged_uo <- exp(predict(m_n_uo, newdata = TPS_g))
TPS_g$treatment <- trt[trt == "open"]
TPS_g$rate_n_open_uo <- exp(predict(m_n_uo, newdata = TPS_g))
biome <- do.call("rbind", tapply(1:nrow(wm), wm$PlotID, function(i)
wm[i[1], c("PlotID", "biome")]))
x_a <- merge(merge(TP_a, TPS_a, by = "PlotID"), biome, by = "PlotID")
x_g <- merge(merge(TP_g, TPS_g, by = "PlotID"), biome, by = "PlotID")
### Results for natural wood for uncaged and open treatments
ret_n_uo <- rbind(x_a, x_g)
no_dowel <- levels(wm$PlotID)[table(d_dowel_uncaged_open$PlotID) == 0]
tmp <- subset(TPS_a, !(PlotID %in% no_dowel))
tmp$rate_n_uncaged_uo <- tmp$rate_n_open_uo <- NULL
tmp$treatment <- trt[trt == "uncaged"]
tmp$rate_d_uncaged_uo <- exp(predict(m_d_uo, newdata = tmp))
tmp$treatment <- trt[trt == "open"]
tmp$rate_d_open_uo <- exp(predict(m_d_uo, newdata = tmp))
### Results for dowels for uncaged and open treatments
ret_d_uo <- merge(merge(TP_a, tmp, by = "PlotID"), biome, by = "PlotID")
Additional models presented in Table S1-1 comparing decomposition rates between open and closed treatments
### Subsetting data to focal treatments and either natural tree species or dowels
d_natural_closed_open <- subset(wm, treatment %in% c("closed", "open") &
dowel == 0)
d_dowel_closed_open <- subset(wm, treatment %in% c("closed", "open") &
dowel == 1)
### Model for natural wood: open vs. closed (for angio / gymno)
m_n_co <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS * clade) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_natural_closed_open,
family = gaussian(link = "log"),
offset = log(drymass_start))
ss <- getME(m_n_co , c("theta","fixef"))
m_n_co <- update(m_n_co, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_n_co)$optinfo$conv$opt
## [1] 0
summary(m_n_co)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS * clade) +
## (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_natural_closed_open
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -8863.5 -8781.9 4445.8 -8891.5 2508
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7120 -0.3539 0.0625 0.4657 6.2847
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.002807 0.05298
## Residual 0.001741 0.04172
## Number of obs: 2522, groups: PlotID, 55
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.078338 0.023221 -3.374 0.000742 ***
## expoyrs:PrecS -0.029418 0.022134 -1.329 0.183817
## expoyrs:cladeangio -0.147431 0.022041 -6.689 2.25e-11 ***
## expoyrs:cladegymno -0.066557 0.024463 -2.721 0.006514 **
## expoyrs:treatmentopen -0.006362 0.005138 -1.238 0.215625
## expoyrs:TempS:PrecS -0.032866 0.020791 -1.581 0.113920
## expoyrs:TempS:cladegymno 0.044702 0.009777 4.572 4.83e-06 ***
## expoyrs:PrecS:cladegymno -0.001696 0.021595 -0.079 0.937413
## expoyrs:TempS:treatmentopen -0.012574 0.005333 -2.358 0.018379 *
## expoyrs:PrecS:treatmentopen -0.020498 0.005477 -3.742 0.000182 ***
## expoyrs:TempS:PrecS:cladegymno 0.004832 0.015669 0.308 0.757773
## expoyrs:TempS:PrecS:treatmentopen -0.026013 0.005465 -4.760 1.94e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:cldn expyrs:cldg expyrs:t ex:TS:PS expyrs:TmpS:c
## expyrs:PrcS -0.312
## expyrs:cldn -0.182 -0.017
## expyrs:cldg -0.149 -0.026 0.879
## expyrs:trtm -0.028 0.010 -0.108 -0.109
## expyr:TS:PS 0.002 -0.328 -0.275 -0.235 0.023
## expyrs:TmpS:c -0.017 0.008 0.003 0.339 -0.013 -0.001
## expyrs:PrcS:c 0.033 -0.051 -0.023 0.088 -0.002 0.033 0.215
## expyrs:TmpS:t -0.104 0.025 -0.029 -0.035 0.307 -0.017 -0.042
## expyrs:PrcS:t 0.015 -0.096 0.010 0.012 -0.081 -0.037 0.011
## expyrs:TmpS:PrcS:c 0.019 -0.013 -0.006 0.087 0.000 -0.007 0.372
## expyrs:TmpS:PrcS:t -0.018 -0.026 0.023 0.024 -0.158 -0.111 -0.006
## expyrs:PrcS:c expyrs:TmpS:t expyrs:PrcS:t expyrs:TmpS:PrcS:c
## expyrs:PrcS
## expyrs:cldn
## expyrs:cldg
## expyrs:trtm
## expyr:TS:PS
## expyrs:TmpS:c
## expyrs:PrcS:c
## expyrs:TmpS:t -0.004
## expyrs:PrcS:t 0.001 -0.143
## expyrs:TmpS:PrcS:c 0.884 -0.016 0.007
## expyrs:TmpS:PrcS:t 0.012 0.117 0.415 0.004
### Transform temperature and precipitation back to degree Celsius and dm per year
Xscale <- model.matrix(m_n_co)
Xorig <- model.matrix(~ 0 + ### intercept zero
expoyrs:(TempC * PrecC * clade) + ### slopes varying in temp and prec
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempC * PrecC),
data = droplevels(d_natural_closed_open))
K <- solve(crossprod(Xorig)) %*% crossprod(Xorig, Xscale)
g <- glht(m_n_co, linfct = K)
### Create main model output table (Table S1-1) for natural wood
### Model output
(tab <- summary(g, test = univariate()))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS * clade) + (0 + expoyrs | PlotID) + expoyrs:treatment +
## expoyrs:treatment:(TempS * PrecS), data = d_natural_closed_open,
## family = gaussian(link = "log"), control = glmerControl(optCtrl = list(maxfun = 2e+06)),
## start = ss, offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:TempC == 0 -9.906e-03 2.936e-03 -3.374 0.000742 ***
## expoyrs:PrecC == 0 -4.344e-03 3.268e-03 -1.329 0.183817
## expoyrs:cladeangio == 0 -1.474e-01 2.204e-02 -6.689 2.25e-11 ***
## expoyrs:cladegymno == 0 -6.656e-02 2.446e-02 -2.721 0.006514 **
## expoyrs:treatmentopen == 0 -6.362e-03 5.138e-03 -1.238 0.215625
## expoyrs:TempC:PrecC == 0 -6.136e-04 3.882e-04 -1.581 0.113920
## expoyrs:TempC:cladegymno == 0 5.653e-03 1.236e-03 4.572 4.83e-06 ***
## expoyrs:PrecC:cladegymno == 0 -2.504e-04 3.189e-03 -0.079 0.937413
## expoyrs:TempC:treatmentopen == 0 -1.590e-03 6.743e-04 -2.358 0.018379 *
## expoyrs:PrecC:treatmentopen == 0 -3.027e-03 8.088e-04 -3.742 0.000182 ***
## expoyrs:TempC:PrecC:cladegymno == 0 9.022e-05 2.925e-04 0.308 0.757773
## expoyrs:TempC:PrecC:treatmentopen == 0 -4.857e-04 1.020e-04 -4.760 1.94e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Relative effects and confidence intervals
### confidence intervals
ci <- confint(g, calpha = univariate_calpha())$confint
### relative effects
exp(ci)
## Estimate lwr upr
## expoyrs:TempC 0.9901429 0.9844608 0.9958578
## expoyrs:PrecC 0.9956658 0.9893087 1.0020639
## expoyrs:cladeangio 0.8629217 0.8264371 0.9010169
## expoyrs:cladegymno 0.9356092 0.8918077 0.9815620
## expoyrs:treatmentopen 0.9936577 0.9837009 1.0037153
## expoyrs:TempC:PrecC 0.9993865 0.9986265 1.0001472
## expoyrs:TempC:cladegymno 1.0056686 1.0032347 1.0081084
## expoyrs:PrecC:cladegymno 0.9997497 0.9935212 1.0060172
## expoyrs:TempC:treatmentopen 0.9984113 0.9970927 0.9997317
## expoyrs:PrecC:treatmentopen 0.9969780 0.9953989 0.9985596
## expoyrs:TempC:PrecC:cladegymno 1.0000902 0.9995170 1.0006638
## expoyrs:TempC:PrecC:treatmentopen 0.9995144 0.9993145 0.9997144
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
### Model for dowels: open vs. closed (for angio / gymno)
m_d_co <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_dowel_closed_open,
family = gaussian(link = "log"),
offset = log(drymass_start))
ss <- getME(m_d_co, c("theta","fixef"))
m_d_co <- update(m_d_co, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_d_co)$optinfo$conv$opt
## [1] 0
summary(m_d_co)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: gaussian ( log )
## Formula: I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS * PrecS) + (0 + expoyrs |
## PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS * PrecS)
## Data: d_dowel_closed_open
## Offset: log(drymass_start)
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## -1474.3 -1434.1 747.1 -1494.3 401
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3975 -0.2563 0.0623 0.3764 3.7149
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID expoyrs 0.010172 0.10086
## Residual 0.001526 0.03906
## Number of obs: 411, groups: PlotID, 55
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## expoyrs:TempS -0.19370 0.04239 -4.570 4.88e-06 ***
## expoyrs:PrecS -0.02235 0.03816 -0.586 0.558025
## expoyrs:treatmentclosed -0.13694 0.03811 -3.594 0.000326 ***
## expoyrs:treatmentopen -0.17653 0.03746 -4.713 2.44e-06 ***
## expoyrs:TempS:PrecS -0.10548 0.03648 -2.891 0.003838 **
## expoyrs:TempS:treatmentopen -0.03817 0.01312 -2.910 0.003617 **
## expoyrs:PrecS:treatmentopen -0.01947 0.01368 -1.423 0.154849
## expoyrs:TempS:PrecS:treatmentopen -0.02456 0.01559 -1.575 0.115202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## exp:TS exp:PS expyrs:trtmntc expyrs:trtmntp ex:TS:PS ex:TS: ex:PS:
## expyrs:PrcS -0.375
## expyrs:trtmntc -0.197 0.006
## expyrs:trtmntp -0.240 0.011 0.954
## expyr:TS:PS 0.143 -0.355 -0.282 -0.290
## expyrs:TmS: -0.241 0.015 -0.103 0.038 -0.070
## expyrs:PrS: -0.020 -0.200 0.002 0.033 -0.179 0.109
## expy:TS:PS: -0.054 -0.108 0.006 0.026 -0.287 0.284 0.630
### Relative effects for comparison between treatments
### Treatment effect when temperature is 15 ° Celsius and
### when precipitation is 13 dm/a
g <- glht(m_d_co, linfct = c("expoyrs:treatmentopen - expoyrs:treatmentclosed = 0"))
summary(g)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS) + (0 + expoyrs | PlotID) + expoyrs:treatment + expoyrs:treatment:(TempS *
## PrecS), data = d_dowel_closed_open, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:treatmentopen - expoyrs:treatmentclosed == 0 -0.03959 0.01146 -3.454 0.000553 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(g)$confint)
## Estimate lwr upr
## expoyrs:treatmentopen - expoyrs:treatmentclosed 0.9611803 0.939825 0.9830209
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
### Compute fitted values for annual wood decomposition rates in open and
### closed treatment for each experimental site, seperately for natural wood
### (angiosperm and gymnosperm seperately were applicable) and dowels
a <- subset(wm, treatment == "closed" & clade == "angio")
TP_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "Temp", "Prec")]
}))
g <- subset(wm, clade != "angio")
TP_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "Temp", "Prec")]
}))
TPS_a <- do.call("rbind", tapply(1:nrow(a), a$PlotID, function(i) {
a[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_a$expoyrs <- 1
clade <- sort(unique(wm$clade))
TPS_a$clade <- clade[1]
TPS_g <- do.call("rbind", tapply(1:nrow(g), g$PlotID, function(i) {
g[i[1], c("PlotID", "TempS", "PrecS")]
}))
TPS_g$expoyrs <- 1
TPS_g$clade <- clade[2]
trt <- sort(unique(wm$treatment))
TPS_a$treatment <- trt[trt == "closed"]
TPS_a$rate_n_closed_co <- exp(predict(m_n_co, newdata = TPS_a))
TPS_a$treatment <- trt[trt == "open"]
TPS_a$rate_n_open_co <- exp(predict(m_n_co, newdata = TPS_a))
TPS_g$treatment <- trt[trt == "closed"]
TPS_g$rate_n_closed_co <- exp(predict(m_n_co, newdata = TPS_g))
TPS_g$treatment <- trt[trt == "open"]
TPS_g$rate_n_open_co <- exp(predict(m_n_co, newdata = TPS_g))
biome <- do.call("rbind", tapply(1:nrow(wm), wm$PlotID, function(i)
wm[i[1], c("PlotID", "biome")]))
x_a <- merge(merge(TP_a, TPS_a, by = "PlotID"), biome, by = "PlotID")
x_g <- merge(merge(TP_g, TPS_g, by = "PlotID"), biome, by = "PlotID")
### Results for natural wood for closed and open treatments
ret_n_co <- rbind(x_a, x_g)
no_dowel <- levels(wm$PlotID)[table(d_dowel_closed_open$PlotID) == 0]
tmp <- subset(TPS_a, !(PlotID %in% no_dowel))
tmp$rate_n_closed_co <- tmp$rate_n_open_co <- NULL
tmp$treatment <- trt[trt == "closed"]
tmp$rate_d_closed_co <- exp(predict(m_d_co, newdata = tmp))
tmp$treatment <- trt[trt == "open"]
tmp$rate_d_open_co <- exp(predict(m_d_co, newdata = tmp))
### Results for dowels for closed and open treatments
ret_d_co <- merge(merge(TP_a, tmp, by = "PlotID"), biome, by = "PlotID")
Merge fitted values from different models to results data set and save for further analyses
ret_n <- merge(ret_n, ret_n_uo[, c("PlotID", "clade",
"rate_n_uncaged_uo", "rate_n_open_uo")],
by = c("PlotID", "clade"))
ret_n <- merge(ret_n, ret_n_co[, c("PlotID", "clade",
"rate_n_closed_co", "rate_n_open_co")],
by = c("PlotID", "clade"))
ret_d <- merge(ret_d, ret_d_uo[, c("PlotID",
"rate_d_uncaged_uo", "rate_d_open_uo")],
by = "PlotID")
ret_d <- merge(ret_d, ret_d_co[, c("PlotID",
"rate_d_closed_co", "rate_d_open_co")],
by = "PlotID")
save(ret_n, file="data/ret_natural.rda")
save(ret_d, file="data/ret_dowel.rda")
### Output data include annual decomposition rates for natural wood (ret_n) and dowels (ret_d):
#"rate_n_uncaged" = "rate_n_uncagedEstimate": decomposition rate for natural
# wood for uncaged treatment from model uncaged vs. closed
#
#"rate_n_uncagedlwr": lower confidence level of decomposition rate for
# natural wood for uncaged treatment from model uncaged vs. closed
#
#"rate_n_uncagedupr": upper confidence level of decomposition rate for
# natural wood for uncaged treatment from model uncaged vs. closed
#
#"rate_n_closed" = "rate_n_closedEstimate": decomposition rate for natural
# wood for closed treatment from model uncaged vs. closed
#
#"rate_n_closedlwr": lower confidence level of decomposition rate for
# natural wood for closed treatment from model uncaged vs. closed
#
#"rate_n_closedupr": upper confidence level of decomposition rate for
# natural wood for closed treatment from model uncaged vs. closed
#
#"rate_n_uncaged_uo": decomposition rate for natural wood for uncaged
# treatment from model uncaged vs. open
#
#"rate_n_open_uo": decomposition rate for natural wood for open treatment
# from model uncaged vs. open
#
#"rate_n_closed_co": decomposition rate for natural wood for closed
# treatment from model closed vs. open
#
#"rate_n_open_co": decomposition rate for natural wood for open treatment
# from model closed vs. open
#
#"rate_d_uncaged": decomposition rate for dowels for uncaged treatment from
# model uncaged vs. closed
#
#"rate_d_closed": decomposition rate for dowels for closed treatment from
# model uncaged vs. closed
#
#"rate_d_uncaged_uo": decomposition rate for dowels for uncaged treatment
# from model uncaged vs. open
#
#"rate_d_open_uo": decomposition rate for dowels for open treatment from
# model uncaged vs. open
#
#"rate_d_closed_co": decomposition rate for dowels for closed treatment from
# model closed vs. open
#
#"rate_d_open_co": decomposition rate for dowels for open treatment from
# model closed vs. open
Reproduce Figure 1 b and c: Boxplot showing annual decomposition of all decomposers and insect contribution as boxplot for different biomes
x <- ret_n
x$trt <- with(x, rate_n_closed - rate_n_uncaged)
### Median of annual woo decomposition rates per biome
tapply((x$trt*100),x$biome,median)
## Boreal Temperate Tropical
## -0.08758852 0.87671480 3.91539712
tapply(((1 - x$rate_n_uncaged) * 100),x$biome,median)
## Boreal Temperate Tropical
## 3.341168 6.345905 28.200202
### Set graphic parameters
cexlab <- 0.58 # 7 pt
cexax <- 0.58 # 7 pt
cexnum <- 0.67 # 8pt
axline <- 1.5
tlength <- -0.2 # ticklength
posx <- 1.4
posx2 <- 1.9
posy <- 0.7
if (TIFF)
tiff(file = "products/Figure1.tif",width = 4.7, height = 2.5,res=300,units="in")
par(mfrow=c(1,2))
par(mai=c(0.3,0.45,0.3,0.0), mgp=c(1.5,0.4,0))
### Fig. 1 panel b: annual wood decomposition rate with all decomposers involved (uncaged treatment)
boxplot(I((1 - rate_n_uncaged) * 100) ~ biome, data = x, axes=FALSE, ylab="", xlab="",ylim=c(-5,55))
axis(1,c(1, 2, 3),c("Boreal", "Temperate", "Tropical"),tick=TRUE,cex.axis=cexlab,tcl=tlength)
axis(2,seq(0,50,10),seq(0,50,10),las=1,tick=TRUE,cex.axis=cexlab,tcl=tlength)
mtext("Annual mass loss [%]",2,cex=cexnum, line=axline)
mtext("b All decomposers",3,cex=cexnum, line=posy,at=posx)
box()
grid(nx = 1, ny = NULL, col = "lightgray", lty=1 ,equilogs = TRUE)
boxplot(I((1 - rate_n_uncaged) * 100) ~ biome, data = x,axes=FALSE, add=TRUE, col="#F5F5F5")
### Fig. 1 panel c: contribution of insects to annual wood decomposition rate
### (closed minus uncaged treatment)
par(mai=c(0.3,0.4,0.3,0.05))
boxplot(I((x$trt) * 100) ~ biome, data = x, ylab = "",axes=FALSE, xlab="",ylim=c(-5,55), fill="white")
axis(1,c(1, 2, 3),c("Boreal", "Temperate", "Tropical"),tick=TRUE,cex.axis=cexlab,tcl=tlength)
axis(2,seq(0,50,10),seq(0,50,10),las=1,tick=TRUE,cex.axis=cexlab,tcl=tlength)
mtext("c Net insect effect",3,cex=cexnum, line=posy,at=posx2)
box()
grid(nx = 1, ny = NULL, col = "lightgray", lty=1 ,equilogs = TRUE)
boxplot(I((x$trt) * 100) ~ biome, data = x,axes=FALSE, add=TRUE, col="#D8B365")
if (TIFF)
dev.off()
Reproduce Figure 2: Annual decomposition rates per site relative to temperature and precipitation
### Set graphic parameters
cexlab <- 7
cextitle <- 8
ptsize <- 3
### Fig. 2 panel a: annual wood decomposition rate with all decomposers involved
### (uncaged treatment)
x <- ret_n
breaks <- floor(min(1 - x$rate_n_uncaged) * 10):ceiling(max(1 - x$rate_n_uncaged) * 10) / 10
x$ct <- cut(1 - x$rate_n_uncaged, breaks = breaks)
levels(x$ct) <- paste((breaks[-length(breaks)] * 100), "% -", (breaks[-1] * 100), "%")
col <- diverging_hcl(nlevels(x$ct), c = 100, l = c(50, 90))
p1 <- ggplot(x,aes(y = Temp, x = Prec, colour = ct)) + geom_point(size=ptsize) +
scale_color_manual(values = col)+
facet_grid(~ clade,labeller = labeller(clade=c(`angio`="Angiosperm", `gymno`="Gymnosperm")))+
theme_bw()+ guides(color=guide_legend(title="Annual mass loss"))+
theme(panel.border = element_rect(colour = "black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_blank() ,
title =element_text(size=cextitle, face="plain"),plot.title = element_text(hjust=-0.2,face="plain"),
axis.text=element_text(size=cexlab), axis.title=element_text(size=cexlab,face="plain"),
legend.text = element_text(size=cexlab),
strip.text.x = element_text(size = cexlab))+
ggtitle("b All decomposers")+xlab("Precipitation [dm / a]")+ylab("Temperature [°C]")
### Fig. 2 panel b: contribution of insects to annual wood decomposition rate (closed minus uncaged treatment)
x <- ret_n
x$trt <- with(x, rate_n_closed - rate_n_uncaged)
breaks <- floor(min(x$trt) * 20):ceiling(max(x$trt) * 20) / 20
x$ct2 <- cut(x$trt, breaks = breaks)
levels(x$ct2) <- paste((breaks[-length(breaks)] * 100), "% -", (breaks[-1] * 100), "%")
col2 <- diverging_hcl(nlevels(x$ct), c = 100, l = c(50, 90))
p2 <- ggplot(x,aes(y = Temp, x = Prec, colour = ct2)) + geom_point(size=ptsize,aes(shape=termite)) +
scale_color_manual(values = col2)+
facet_grid(~ clade,labeller = labeller(clade=c(`angio`="Angiosperm", `gymno`="Gymnosperm")))+
theme_bw()+
scale_shape_discrete(name="Termites",breaks=c( "17","16", "18"),labels=c( "Yes", "No","Maybe"))+
guides(color=guide_legend(title="Annual mass loss"))+
theme(panel.border = element_rect(colour = "black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_blank() ,
title =element_text(size=cextitle, face='plain'),plot.title = element_text(hjust=-0.2),
axis.text=element_text(size=cexlab), axis.title=element_text(size=cexlab,face="plain"),
legend.text = element_text(size=cexlab),
strip.text.x = element_text(size = cexlab),
legend.spacing.y = unit(0.001, 'cm'))+
ggtitle("c Net insect effect")+xlab("Precipitation [dm / a]")+ylab("Temperature [°C]")
if (TIFF)
tiff(file = "products/Figure2.tif",width = 4.7, height = 5,res=300,units="in")
plot(ggarrange(p1,p2 , ncol = 1, nrow = 2))
if (TIFF)
dev.off()
Reproduce Extended Data Figure 3: Pairwise comparison of decomposition rates among all three treatments
x <- ret_n
x$trt <- with(x, rate_n_closed - rate_n_uncaged)
x$Biome <- ordered(as.character(x$biome),
levels = rev(levels(x$biome)),
labels = rev(levels(x$biome)))
p1 <- ggplot(x,aes(y = (1-rate_n_closed)*100,
x = (1-rate_n_uncaged)*100, colour = Biome)) +
geom_abline(intercept = 0, slope = 1, color = "grey", lty = 2)+ geom_point(size=3)+
xlim(-20, 60) +ylim(-20, 60)+
xlab("Annual mass loss [%] - uncaged")+ylab("Annual mass loss [%] - closed cage")+
theme(legend.position="none")
p2 <- ggplot(x,aes(y = (1-rate_n_open_uo)*100,
x = (1-rate_n_uncaged_uo)*100, colour = Biome)) +
geom_abline(intercept = 0, slope = 1, color = "grey", lty = 2)+ geom_point(size=3)+
xlim(-20, 60) +ylim(-20, 60)+
xlab("Annual mass loss [%] - uncaged")+ylab("Annual mass loss [%] - open cage")+
theme(legend.position="none")
p3 <- ggplot(x,aes(y = (1-rate_n_closed)*100 ,
x = (1-rate_n_open_uo)*100 , colour = Biome)) +
geom_abline(intercept = 0, slope = 1, color = "grey", lty = 2)+ geom_point(size=3)+
xlim(-20, 60) +ylim(-20, 60)+
xlab("Annual mass loss [%] - open cage")+ylab("Annual mass loss [%] - closed cage")+
theme(legend.position = c(1.4, .7))
if (TIFF)
tiff(filename = "products/ExtendedDataFigure3.tif",width = 160, height = 180,units = "mm", res = 300,compression = "lzw")
plot(p1+ p2 +p3+ plot_layout(nrow = 2, byrow = TRUE))
if (TIFF)
dev.off()
Model presented in Supplementary Table S1-1 comparing decomposition rates of natural wood among all three treatments (open, closed, uncaged)
### Subset to natural wood only
d_natural <- subset(wm, dowel == 0)
m_n_all <- glmer(I(drymass_end + .01) ~ 0 + ### intercept zero
expoyrs:(TempS * PrecS * clade) + ### slopes varying in temp and prec
(0 + expoyrs | PlotID) + ### space-varying slopes
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempS * PrecS),
data = d_natural,
family = gaussian(link = "log"),
offset = log(drymass_start))
### increase iterations => model converges
ss <- getME(m_n_all, c("theta","fixef"))
m_n_all <- update(m_n_all, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(m_n_all)$optinfo$conv$opt
## [1] 0
### Transform temperature and precipitation back to degree Celsius and dm per year
Xscale <- model.matrix(m_n_all)
Xorig <- model.matrix(~ 0 + ### intercept zero
expoyrs:(TempC * PrecC * clade) + ### slopes varying in temp and prec
expoyrs:treatment + ### treatment + interactions
expoyrs:treatment:(TempC * PrecC),
data = droplevels(d_natural))
K <- solve(crossprod(Xorig)) %*% crossprod(Xorig, Xscale)
g <- glht(m_n_all, linfct = K)
### Create main model output table (Table S1-1) for natural wood
### Model output
(tab <- summary(g, test = univariate()))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS * clade) + (0 + expoyrs | PlotID) + expoyrs:treatment +
## expoyrs:treatment:(TempS * PrecS), data = d_natural, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:TempC == 0 -1.133e-02 3.018e-03 -3.755 0.000173 ***
## expoyrs:PrecC == 0 -2.981e-03 3.319e-03 -0.898 0.369092
## expoyrs:cladeangio == 0 -1.533e-01 2.262e-02 -6.775 1.24e-11 ***
## expoyrs:cladegymno == 0 -8.402e-02 2.421e-02 -3.470 0.000520 ***
## expoyrs:treatmentopen == 0 -7.224e-03 5.573e-03 -1.296 0.194935
## expoyrs:treatmentuncaged == 0 -2.952e-02 5.584e-03 -5.285 1.25e-07 ***
## expoyrs:TempC:PrecC == 0 -5.403e-04 4.000e-04 -1.351 0.176742
## expoyrs:TempC:cladegymno == 0 5.177e-03 1.014e-03 5.108 3.26e-07 ***
## expoyrs:PrecC:cladegymno == 0 -8.984e-04 2.812e-03 -0.320 0.749329
## expoyrs:TempC:treatmentopen == 0 -1.965e-03 7.453e-04 -2.636 0.008380 **
## expoyrs:TempC:treatmentuncaged == 0 -4.271e-03 7.286e-04 -5.862 4.56e-09 ***
## expoyrs:PrecC:treatmentopen == 0 -3.190e-03 9.047e-04 -3.526 0.000421 ***
## expoyrs:PrecC:treatmentuncaged == 0 -5.162e-03 9.022e-04 -5.722 1.06e-08 ***
## expoyrs:TempC:PrecC:cladegymno == 0 1.719e-05 2.570e-04 0.067 0.946692
## expoyrs:TempC:PrecC:treatmentopen == 0 -4.575e-04 1.133e-04 -4.040 5.35e-05 ***
## expoyrs:TempC:PrecC:treatmentuncaged == 0 -7.255e-04 1.106e-04 -6.561 5.35e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Relative effects and confidence intervals
### confidence intervals
ci <- confint(g, calpha = univariate_calpha())$confint
### relative effects
exp(ci)
## Estimate lwr upr
## expoyrs:TempC 0.9887323 0.9829014 0.9945977
## expoyrs:PrecC 0.9970235 0.9905591 1.0035301
## expoyrs:cladeangio 0.8579151 0.8207133 0.8968031
## expoyrs:cladegymno 0.9194128 0.8768014 0.9640951
## expoyrs:treatmentopen 0.9928023 0.9820164 1.0037067
## expoyrs:treatmentuncaged 0.9709154 0.9603465 0.9816006
## expoyrs:TempC:PrecC 0.9994598 0.9986766 1.0002437
## expoyrs:TempC:cladegymno 1.0051909 1.0031960 1.0071899
## expoyrs:PrecC:cladegymno 0.9991020 0.9936112 1.0046231
## expoyrs:TempC:treatmentopen 0.9980369 0.9965800 0.9994960
## expoyrs:TempC:treatmentuncaged 0.9957378 0.9943169 0.9971608
## expoyrs:PrecC:treatmentopen 0.9968148 0.9950489 0.9985839
## expoyrs:PrecC:treatmentuncaged 0.9948510 0.9930933 0.9966118
## expoyrs:TempC:PrecC:cladegymno 1.0000172 0.9995135 1.0005211
## expoyrs:TempC:PrecC:treatmentopen 0.9995426 0.9993208 0.9997645
## expoyrs:TempC:PrecC:treatmentuncaged 0.9992748 0.9990582 0.9994914
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
### all-pair comparison of main effects only
m_n_all.post<-glht(m_n_all, linfct = c("expoyrs:treatmentopen = 0",
"expoyrs:treatmentuncaged = 0",
"expoyrs:treatmentuncaged - expoyrs:treatmentopen = 0"))
summary(m_n_all.post)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = I(drymass_end + 0.01) ~ 0 + expoyrs:(TempS *
## PrecS * clade) + (0 + expoyrs | PlotID) + expoyrs:treatment +
## expoyrs:treatment:(TempS * PrecS), data = d_natural, family = gaussian(link = "log"),
## control = glmerControl(optCtrl = list(maxfun = 2e+06)), start = ss,
## offset = log(drymass_start))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## expoyrs:treatmentopen == 0 -0.007224 0.005573 -1.296 0.397313
## expoyrs:treatmentuncaged == 0 -0.029516 0.005584 -5.285 < 1e-04 ***
## expoyrs:treatmentuncaged - expoyrs:treatmentopen == 0 -0.022292 0.005505 -4.050 0.000154 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Reproduce Extended Data Figure 2 panel a: differences in decomposition rate between all three treatments
m_n_all.post_ci <- data.frame(confint(m_n_all.post)$confint)
theme_set(theme_classic())
p1 <- ggplot(m_n_all.post_ci, aes(x = Estimate, y = c(1,2,3), xmin = lwr, xmax = upr)) +
geom_point(col="black",size=3) + geom_errorbarh(height=.2)+
geom_vline(xintercept=0, linetype="dashed", color = "black")+
xlab("Linear function")+ylab("")+
ggtitle("a Difference in decomposition rates between treatments")+
scale_y_continuous(breaks=c(1,2,3),labels=c("open vs. closed","uncaged vs. closed","uncaged vs. open"))+
theme(title =element_text(size=cextitle, face='plain'),plot.title = element_text(hjust=0),
axis.text=element_text(size=cexlab, color="black"), axis.title=element_text(size=cexlab,face="plain"),
panel.border = element_rect(colour = "black", fill=NA))
Reproduce Extended Data Figure 2 panel b: Differences in insect colonization (0 = no insect marks, 1 = insect marks) between all three treatments
natural <- subset(wm, wm$dowel == 0)
### Model for pairwise comparisons of insect colonization of natural wood
### between all three treatments (open, closed, uncaged)
mod_marks <- glmer(marks_total ~
expoyrs + (TempS * PrecS * clade) +
treatment * TempS * PrecS + (1| PlotID) ,
data = natural, family = "binomial")
ss <- getME(mod_marks, c("theta","fixef"))
mod_marks <- update(mod_marks, start = ss,
control = glmerControl(optCtrl = list(maxfun = 2e6)))
### check convergence (0 means model converged)
summary(mod_marks)$optinfo$conv$opt
## [1] 0
summary(mod_marks)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: binomial ( logit )
## Formula: marks_total ~ expoyrs + (TempS * PrecS * clade) + treatment *
## TempS * PrecS + (1 | PlotID)
## Data: natural
## Control: glmerControl(optCtrl = list(maxfun = 2e+06))
##
## AIC BIC logLik deviance df.resid
## 2753.2 2863.7 -1358.6 2717.2 3412
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1363 -0.4448 -0.1361 0.4286 12.1063
##
## Random effects:
## Groups Name Variance Std.Dev.
## PlotID (Intercept) 4.678 2.163
## Number of obs: 3430, groups: PlotID, 54
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.485082 0.367628 -9.480 < 2e-16 ***
## expoyrs 0.399333 0.076595 5.214 1.85e-07 ***
## TempS 1.597340 0.348358 4.585 4.53e-06 ***
## PrecS 0.776354 0.334827 2.319 0.020413 *
## cladegymno 0.811922 0.206570 3.930 8.48e-05 ***
## treatmentopen 1.417549 0.132816 10.673 < 2e-16 ***
## treatmentuncaged 1.790294 0.135439 13.218 < 2e-16 ***
## TempS:PrecS 0.161073 0.339403 0.475 0.635088
## TempS:cladegymno -1.107885 0.207705 -5.334 9.61e-08 ***
## PrecS:cladegymno -0.047207 0.399132 -0.118 0.905850
## TempS:treatmentopen 0.011093 0.127430 0.087 0.930628
## TempS:treatmentuncaged -0.008045 0.127883 -0.063 0.949842
## PrecS:treatmentopen -0.457531 0.130009 -3.519 0.000433 ***
## PrecS:treatmentuncaged -0.524621 0.132689 -3.954 7.69e-05 ***
## TempS:PrecS:cladegymno 0.463649 0.367824 1.261 0.207482
## TempS:PrecS:treatmentopen -0.099888 0.142467 -0.701 0.483222
## TempS:PrecS:treatmentuncaged -0.066305 0.142262 -0.466 0.641162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
Transform temperature and precipitation back to degree Celsius and dm per year.
TempC is temperature in °C - 15 and PrecC is precipitation in dm / a - 13.
Effects of TempC can be interpreted as effect of increasing temperature by 1 °C when precipitation is 13 dm/a.
Effects of PreC can be interpreted as effect an additional dm / a when temperature is 15 °C.
Xscale <- model.matrix(mod_marks)
Xorig <- model.matrix(~ expoyrs + (TempC * PrecC * clade) +
treatment * TempC * PrecC,
data = droplevels(natural)[!is.na(natural$marks_total),])
K <- solve(crossprod(Xorig)) %*% crossprod(Xorig, Xscale)
g <- glht(mod_marks, linfct = K)
(tab <- summary(g, test = univariate()))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmer(formula = marks_total ~ expoyrs + (TempS * PrecS * clade) +
## treatment * TempS * PrecS + (1 | PlotID), data = natural,
## family = "binomial", control = glmerControl(optCtrl = list(maxfun = 2e+06)),
## start = ss)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) == 0 -3.485082 0.367628 -9.480 < 2e-16 ***
## expoyrs == 0 0.399333 0.076595 5.214 1.85e-07 ***
## TempC == 0 0.201986 0.044050 4.585 4.53e-06 ***
## PrecC == 0 0.114630 0.049438 2.319 0.020413 *
## cladegymno == 0 0.811922 0.206570 3.930 8.48e-05 ***
## treatmentopen == 0 1.417549 0.132816 10.673 < 2e-16 ***
## treatmentuncaged == 0 1.790294 0.135439 13.218 < 2e-16 ***
## TempC:PrecC == 0 0.003007 0.006337 0.475 0.635088
## TempC:cladegymno == 0 -0.140094 0.026265 -5.334 9.61e-08 ***
## PrecC:cladegymno == 0 -0.006970 0.058933 -0.118 0.905850
## TempC:treatmentopen == 0 0.001403 0.016114 0.087 0.930628
## TempC:treatmentuncaged == 0 -0.001017 0.016171 -0.063 0.949842
## PrecC:treatmentopen == 0 -0.067555 0.019196 -3.519 0.000433 ***
## PrecC:treatmentuncaged == 0 -0.077461 0.019592 -3.954 7.69e-05 ***
## TempC:PrecC:cladegymno == 0 0.008657 0.006868 1.261 0.207482
## TempC:PrecC:treatmentopen == 0 -0.001865 0.002660 -0.701 0.483222
## TempC:PrecC:treatmentuncaged == 0 -0.001238 0.002656 -0.466 0.641162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Post-hoc comparisons
mod_marks.post<-glht(mod_marks, mcp(treatment = "Tukey"))
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found -- default contrast
## might be inappropriate
summary(mod_marks.post)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = marks_total ~ expoyrs + (TempS * PrecS * clade) +
## treatment * TempS * PrecS + (1 | PlotID), data = natural,
## family = "binomial", control = glmerControl(optCtrl = list(maxfun = 2e+06)),
## start = ss)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## open - closed == 0 1.4175 0.1328 10.673 < 0.001 ***
## uncaged - closed == 0 1.7903 0.1354 13.218 < 0.001 ***
## uncaged - open == 0 0.3727 0.1208 3.085 0.00572 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
mod_marks.post <- data.frame(confint(mod_marks.post)$confint)
theme_set(theme_classic())
p2 <- ggplot(mod_marks.post, aes(x = -Estimate, y = c(1,2,3), xmin = -lwr, xmax = -upr)) +
geom_point(col="black",size=3) + geom_errorbarh(height=.2)+
geom_vline(xintercept=0, linetype="dashed", color = "black")+
xlab("Logarithmic effect")+ylab("")+
ggtitle("b Differences in insect colonization rates between treatments")+
scale_y_continuous(breaks=c(1,2,3),labels=c("","",""))+
theme(title =element_text(size=cextitle, face='plain'),plot.title = element_text(hjust=0),
axis.text=element_text(size=cexlab), axis.title=element_text(size=cexlab,face="plain"),
panel.border = element_rect(colour = "black", fill=NA))
if (TIFF)
tiff(filename = "products/ExtendedDataFigure2.tif",width = 9, height = 4.5,res=300,units="in",compression = "lzw")
plot(p1 + p2 +plot_layout(nrow = 1, byrow = TRUE))
if (TIFF)
dev.off()