library("ggplot2")
library("readr")
library("DESeq2")
library("phyloseq")
library('ggrepel')
library('ggalluvial')
library("microeco")
library("dplyr")
library("tidyr")
library('tibble')
library('pals')
library("rstatix")
library("cowplot")
library("vegan")
library("car")
library("mgcv")
sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Berlin
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] mgcv_1.9-4                  nlme_3.1-168               
##  [3] car_3.1-5                   carData_3.0-6              
##  [5] vegan_2.7-2                 permute_0.9-10             
##  [7] cowplot_1.2.0               rstatix_1.1.0              
##  [9] pals_1.10                   tibble_3.3.1               
## [11] tidyr_1.3.2                 dplyr_1.2.1                
## [13] microeco_2.0.0              ggalluvial_0.12.5          
## [15] ggrepel_0.9.6               phyloseq_1.52.0            
## [17] DESeq2_1.48.2               SummarizedExperiment_1.38.1
## [19] Biobase_2.68.0              MatrixGenerics_1.20.0      
## [21] matrixStats_1.5.0           GenomicRanges_1.60.0       
## [23] GenomeInfoDb_1.44.3         IRanges_2.42.0             
## [25] S4Vectors_0.46.0            BiocGenerics_0.54.1        
## [27] generics_0.1.4              readr_2.1.6                
## [29] ggplot2_4.0.2              
## 
## loaded via a namespace (and not attached):
##  [1] rlang_1.2.0             magrittr_2.0.3          ade4_1.7-23            
##  [4] otel_0.2.0              compiler_4.5.0          vctrs_0.7.3            
##  [7] maps_3.4.3              reshape2_1.4.5          stringr_1.6.0          
## [10] pkgconfig_2.0.3         crayon_1.5.3            fastmap_1.2.0          
## [13] backports_1.5.0         XVector_0.48.0          rmarkdown_2.30         
## [16] tzdb_0.5.0              UCSC.utils_1.4.0        purrr_1.2.1            
## [19] xfun_0.61               cachem_1.1.0            jsonlite_2.0.0         
## [22] biomformat_1.36.0       rhdf5filters_1.20.0     DelayedArray_0.34.1    
## [25] Rhdf5lib_1.30.0         BiocParallel_1.42.2     broom_1.0.12           
## [28] parallel_4.5.0          cluster_2.1.8.2         R6_2.6.1               
## [31] bslib_0.10.0            stringi_1.8.7           RColorBrewer_1.1-3     
## [34] jquerylib_0.1.4         Rcpp_1.1.0              iterators_1.0.14       
## [37] knitr_1.51              Matrix_1.7-4            splines_4.5.0          
## [40] igraph_2.2.2            tidyselect_1.2.1        rstudioapi_0.18.0      
## [43] dichromat_2.0-0.1       abind_1.4-8             yaml_2.3.10            
## [46] codetools_0.2-20        lattice_0.22-9          plyr_1.8.9             
## [49] withr_3.0.2             S7_0.2.1                evaluate_1.0.5         
## [52] survival_3.8-6          Biostrings_2.76.0       pillar_1.11.1          
## [55] foreach_1.5.2           hms_1.1.4               scales_1.4.0           
## [58] glue_1.8.0              mapproj_1.2.12          tools_4.5.0            
## [61] data.table_1.18.2.1     locfit_1.5-9.12         rhdf5_2.52.1           
## [64] grid_4.5.0              ape_5.8-1               colorspace_2.1-2       
## [67] GenomeInfoDbData_1.2.14 Formula_1.2-5           cli_3.6.5              
## [70] S4Arrays_1.8.1          gtable_0.3.6            sass_0.4.10            
## [73] digest_0.6.37           SparseArray_1.8.1       farver_2.1.2           
## [76] htmltools_0.5.8.1       multtest_2.64.0         lifecycle_1.0.5        
## [79] httr_1.4.8              MASS_7.3-65

Data preparation

sample_info_no_R <- read.table("Alaska_info_no_R.txt", row.names=1, header=T, check.names=F) # output of BAI_detBAI_RCS.Rmd: raw BAI, detrended BAI (detBAI) and RCS BAI (rcsBAI)
tax_tab <- read_delim("Alaska_taxonomy.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)
## New names:
## Rows: 467 Columns: 18
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "\t" chr
## (16): OTU, Species...4, ID1, ID2, reps/refs, Kingdom, Phylum, Class, Ord... dbl
## (2): size, match
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `Species` -> `Species...4`
## • `Species` -> `Species...14`
count_tab <-read.table("Alaska_counts.txt", header=T, row.names=1, check.names=F)
count_tab <- count_tab[rownames(count_tab) %in% tax_tab$OTU, ]

Treatment of NTC

NTC <-rowSums(count_tab[, 1:8]) # sum of NTC counts for each OTU
sample_counts <- rowSums(count_tab[, 9:143]) # sum of sample counts for each OTU

# normalize sample counts by dividing the samples' total by 17 – as there are 17 as many samples (135) as there are NTCs (8) (135/8=17)
norm_sample_OTU_counts <- sample_counts/17

# which OTUs are deemed likely contaminants based on the threshold noted above:
blank_OTUs <- names(NTC[NTC *10 > norm_sample_OTU_counts]) # OTUs which have more counts than 1/10 of the sample counts
filt_count_tab <- count_tab[!rownames(count_tab) %in% blank_OTUs, -c(1:8)] # removing blank_OTUs and the blank samples from further analysis
filt_count_tab <- filt_count_tab[rowSums(filt_count_tab == 0) != ncol(filt_count_tab), ] # delete OTUs with 0 counts

Normalizing for sampling depth (variance stabilizing transformation)

#reorder the row and column names
sample_info_no_R<-sample_info_no_R[order(rownames(sample_info_no_R)), ]
filt_count_tab<- filt_count_tab[,order(names(filt_count_tab))]
filt_count_tab_no_r <- filt_count_tab[colnames(filt_count_tab) %in% rownames(sample_info_no_R) ]

# create a DESeq2 object
deseq_counts <- DESeqDataSetFromMatrix(filt_count_tab_no_r, colData = sample_info_no_R, design = ~ Habitat) 
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
deseq_counts <- estimateSizeFactors(deseq_counts,type = "poscounts")
deseq_counts_vst <- varianceStabilizingTransformation(deseq_counts,blind = TRUE)
vst_trans_count_tab <- assay(deseq_counts_vst) #save as matrix
vst_trans_count_tab_df0 <- as.data.frame(vst_trans_count_tab)
#transform lowest values (absent counts) to 0
vst_trans_count_tab_df <- sweep(vst_trans_count_tab_df0, 2, apply(vst_trans_count_tab, 2, min), "-")

Functional guilds relative abundancy

tax_table<-as.data.frame(tax_tab)
row.names(tax_table) <- tax_table$OTU
tax_table <- tax_table[, -1]
tax_table<- tax_table[order(rownames(tax_table)),]

# create microeco file
meco_fungi <- microtable$new(sample_table = sample_info_no_R, otu_table = filt_count_tab_no_r, tax_table = tax_table)
mycol <- c("#A6CEE3" ,"#1F78B4" ,"#B2DF8A" ,"#33A02C" ,"#FB9A99", "#E31A1C" ,"#FDBF6F" ,"#FF7F00" ,"#CAB2D6" ,"#6A3D9A" ,"#FFFF99" ,"#B15928" , "#6fdc8c", "green4", "#004144")

t1 <- trans_abund$new(dataset = meco_fungi, taxrank = "Guild", ntaxa = 15)
## No taxa_abund list found. Calculate relative abundance with cal_abund function ...
## The result is stored in object$taxa_abund ...
## The transformed abundance data is stored in object$data_abund ...
t1$plot_bar(bar_type = "notfull", use_alluvium = TRUE,xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, color_values = mycol)
## Warning: The `bar_type` argument of `plot_bar()` is deprecated as of microeco 1.7.0.
## ℹ Please use the `bar_full` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

palette <- c("#EE7733","#CC6677", "#009988","#337538","#33BBEE", "#5E4FA2")
t1$plot_box(group = "Plot",  xtext_angle = 90, color_values = palette) + ylab("Relative abundance (%)")

Setup for the guild comparisons

abund_rel <- t1$data_abund  # relative abundance data

sapro_guilds <- c("Wood Saprotroph","Plant Saprotroph","Undefined Saprotroph")

# sum the saprotroph guilds within each tree -> one row per tree
sapro_rel_hab <- aggregate(Abundance ~ Sample + Site + Habitat, data = abund_rel[abund_rel$Taxonomy %in% sapro_guilds,], FUN = sum)

anyDuplicated(sapro_rel_hab$Sample) # should be 0
## [1] 0

Treeline vs forest - relative abundance data

# ECM
ECM_hab_IN <- wilcox_test(abund_rel[abund_rel$Taxonomy == "Ectomycorrhizal" & abund_rel$Site == "Interior",], Abundance ~ Habitat, detailed = TRUE)
ECM_hab_AR <- wilcox_test(abund_rel[abund_rel$Taxonomy == "Ectomycorrhizal" & abund_rel$Site == "Alaska_Range",], Abundance ~ Habitat, detailed = TRUE)
ECM_hab_BR <- wilcox_test(abund_rel[abund_rel$Taxonomy == "Ectomycorrhizal" & abund_rel$Site == "Brooks_Range",], Abundance ~ Habitat, detailed = TRUE)

ECM_hab <- rbind(ECM_hab_IN, ECM_hab_AR, ECM_hab_BR)
ECM_hab$Site <- c("Interior","Alaska_Range","Brooks_Range")
ECM_hab$p.adj <- p.adjust(ECM_hab$p, method = "BH")
ECM_hab
## # A tibble: 3 × 14
##   estimate .y.    group1 group2    n1    n2 statistic       p conf.low conf.high
##      <dbl> <chr>  <chr>  <chr>  <int> <int>     <dbl>   <dbl>    <dbl>     <dbl>
## 1    -1.40 Abund… forest treel…    19    11        94 0.672      -12.1     11.2 
## 2    -8.09 Abund… forest treel…    19    21       143 0.130      -19.9      2.40
## 3   -13.1  Abund… forest treel…    21    20       111 0.00914    -21.3     -3.05
## # ℹ 4 more variables: method <chr>, alternative <chr>, Site <chr>, p.adj <dbl>
# Saprotrophs
SAP_hab_IN <- wilcox_test(sapro_rel_hab[sapro_rel_hab$Site == "Interior",], Abundance ~ Habitat, detailed = TRUE)
SAP_hab_AR <- wilcox_test(sapro_rel_hab[sapro_rel_hab$Site == "Alaska_Range",], Abundance ~ Habitat, detailed = TRUE)
SAP_hab_BR <- wilcox_test(sapro_rel_hab[sapro_rel_hab$Site == "Brooks_Range",], Abundance ~ Habitat, detailed = TRUE)

SAP_hab <- rbind(SAP_hab_IN, SAP_hab_AR, SAP_hab_BR)
SAP_hab$Site <- c("Interior","Alaska_Range","Brooks_Range")
SAP_hab$p.adj <- p.adjust(SAP_hab$p, method = "BH")
SAP_hab
## # A tibble: 3 × 14
##   estimate .y.    group1 group2    n1    n2 statistic       p conf.low conf.high
##      <dbl> <chr>  <chr>  <chr>  <int> <int>     <dbl>   <dbl>    <dbl>     <dbl>
## 1    -2.40 Abund… forest treel…    19    11        88 4.97e-1   -10.5       4.26
## 2     2.97 Abund… forest treel…    19    21       234 3.61e-1    -5.36     14.3 
## 3    14.5  Abund… forest treel…    21    20       343 3.25e-4     5.99     21.3 
## # ℹ 4 more variables: method <chr>, alternative <chr>, Site <chr>, p.adj <dbl>

Guild abundance vs tree growth

Data preparation

# one value per tree, from the relative abundance data
ECM_tree <- abund_rel[abund_rel$Taxonomy == "Ectomycorrhizal", c("Sample","Abundance")]
names(ECM_tree)[2] <- "ECM"

SAP_tree <- aggregate(Abundance ~ Sample, data = abund_rel[abund_rel$Taxonomy %in% sapro_guilds,], FUN = sum)
names(SAP_tree)[2] <- "SAP"

meta_bai <- sample_info_no_R
meta_bai$Sample <- rownames(meta_bai)

guild_bai <- merge(ECM_tree, SAP_tree, by = "Sample")
guild_bai <- merge(guild_bai, meta_bai, by = "Sample")
guild_bai <- guild_bai[is.finite(guild_bai$BAI_5y), ] # remove samples with NA (same 5 trees for all BAI variables)
nrow(guild_bai) # should be 106
## [1] 106

Model check: Gaussian vs Gamma (ECM, BAI 10 years)

#Check for multicollinearity
vif(lm(log(BAI_10y) ~ ECM + Habitat + Site + log(DBH), data = guild_bai))
##              GVIF Df GVIF^(1/(2*Df))
## ECM      1.325681  1        1.151382
## Habitat  1.089585  1        1.043832
## Site     1.167916  2        1.039568
## log(DBH) 1.147113  1        1.071034
vif(lm(log(BAI_10y) ~ SAP + Habitat + Site + log(DBH), data = guild_bai))
##              GVIF Df GVIF^(1/(2*Df))
## SAP      1.280405  1        1.131550
## Habitat  1.107038  1        1.052159
## Site     1.186463  2        1.043671
## log(DBH) 1.096586  1        1.047180
# same model, two error distributions
m_gaus <- gam(BAI_10y ~ ECM + Habitat + Site + log(DBH),
  family = gaussian(link = "log"),
  method = "REML",
  data = guild_bai)

m_gamma <- gam(BAI_10y ~ ECM + Habitat + Site + log(DBH),
  family = Gamma(link = "log"),
  method = "REML",
  data = guild_bai)

AIC(m_gaus, m_gamma)  # lower = better
##         df      AIC
## m_gaus   7 592.4651
## m_gamma  7 527.4465
gam.check(m_gaus)     # "Resids vs. linear pred." fans out 

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-2.353623e-09,-2.353623e-09]
## (score 302.3649 & scale 14.55027).
## Hessian positive definite, eigenvalue range [50,50].
## Model rank =  6 / 6
gam.check(m_gamma)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-0.000190216,-0.000190216]
## (score 267.897 & scale 0.3361154).
## Hessian positive definite, eigenvalue range [55.74499,55.74499].
## Model rank =  6 / 6
summary(m_gamma)
## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## BAI_10y ~ ECM + Habitat + Site + log(DBH)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.871952   0.523894  -5.482 3.16e-07 ***
## ECM               0.004596   0.003890   1.181 0.240251    
## Habitattreeline   0.463872   0.117894   3.935 0.000154 ***
## SiteBrooks_Range -0.093960   0.141154  -0.666 0.507162    
## SiteInterior      1.032363   0.147818   6.984 3.22e-10 ***
## log(DBH)          1.334744   0.136653   9.767 3.20e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.801   Deviance explained = 69.5%
## -REML =  267.9  Scale est. = 0.33612   n = 106
# prediction line at the median DBH of each plot, so it shows the ECM effect only
nd <- expand.grid(ECM = seq(min(guild_bai$ECM), max(guild_bai$ECM), length.out = 50),
                  Habitat = unique(guild_bai$Habitat),
                  Site = unique(guild_bai$Site))
med_dbh <- aggregate(DBH ~ Habitat + Site, data = guild_bai, FUN = median)
nd <- merge(nd, med_dbh, by = c("Habitat", "Site"))
nd$pred <- predict(m_gamma, newdata = nd, type = "response")

ggplot(guild_bai, aes(x = ECM, y = BAI_10y)) +
  geom_point(alpha = 0.6, color = "#EE7733") +
  geom_line(data = nd, aes(y = pred), color = "#009988", linewidth = 1) +
  facet_grid(Habitat ~ Site) +
  scale_y_log10() +
  theme_bw() + labs(x = "Relative abundance of ECM fungi (%)", y = "Mean BAI over 10 years (cm²/year, log scale)")

Does the ECM effect differ between sites? -> no

# both models: same family and same covariates, only the interaction differs
m_ecm <- gam(rcsBAI_10y ~ ECM + Habitat + Site , family = Gamma(link = "log"), method = "REML", data = guild_bai)
m_ecm_int <- gam(rcsBAI_10y ~ ECM * Site + Habitat , family = Gamma(link = "log"), method = "REML", data = guild_bai) #-> no
anova(m_ecm, m_ecm_int, test = "F")
## Analysis of Deviance Table
## 
## Model 1: rcsBAI_10y ~ ECM + Habitat + Site
## Model 2: rcsBAI_10y ~ ECM * Site + Habitat
##   Resid. Df Resid. Dev Df Deviance      F Pr(>F)
## 1       101     62.905                          
## 2        99     62.218  2  0.68694 0.5239 0.5938
# slope within each site
guild_bai$Site_IN <- relevel(factor(guild_bai$Site), ref = "Interior")
guild_bai$Site_AR <- relevel(factor(guild_bai$Site), ref = "Alaska_Range")
guild_bai$Site_BR <- relevel(factor(guild_bai$Site), ref = "Brooks_Range")

summary(gam(rcsBAI_10y ~ ECM * Site_IN + Habitat , family = Gamma(link = "log"), method = "REML", data = guild_bai))$p.table["ECM",] #-> ns
##     Estimate   Std. Error      t value     Pr(>|t|) 
## -0.008059325  0.009446967 -0.853112507  0.395656259
summary(gam(rcsBAI_10y ~ ECM * Site_AR + Habitat , family = Gamma(link = "log"), method = "REML", data = guild_bai))$p.table["ECM",] #-> ns
##    Estimate  Std. Error     t value    Pr(>|t|) 
## 0.005547571 0.008726679 0.635702376 0.526437215
summary(gam(rcsBAI_10y ~ ECM * Site_BR + Habitat , family = Gamma(link = "log"), method = "REML", data = guild_bai))$p.table["ECM",] #-> ns
##     Estimate   Std. Error      t value     Pr(>|t|) 
## -0.005489628  0.008486473 -0.646868038  0.519214337

Guild abundance vs raw BAI

## Choice of error distribution: Gamma vs Gaussian, both with a log link.
## Positive dAIC = the Gamma fit is better. Checked once here instead of refitting
## a Gaussian model beside every model below.

m_chk_bai <- gam(BAI_10y ~ ECM + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
g_chk_bai <- gam(BAI_10y ~ ECM + Habitat + Site + log(DBH), family = gaussian(link = "log"), method = "REML", data = guild_bai)

m_chk_det <- gam(detBAI_10y ~ ECM + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
g_chk_det <- gam(detBAI_10y ~ ECM + Habitat + Site, family = gaussian(link = "log"), method = "REML", data = guild_bai)

m_chk_rcs <- gam(rcsBAI_10y ~ ECM + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
g_chk_rcs <- gam(rcsBAI_10y ~ ECM + Plot, family = gaussian(link = "log"), method = "REML", data = guild_bai)

data.frame(Model = c("BAI_10y", "detBAI_10y", "rcsBAI_10y"),
           dAIC  = c(AIC(g_chk_bai) - AIC(m_chk_bai),
                     AIC(g_chk_det) - AIC(m_chk_det),
                     AIC(g_chk_rcs) - AIC(m_chk_rcs)))
##        Model     dAIC
## 1    BAI_10y 65.01858
## 2 detBAI_10y 13.56201
## 3 rcsBAI_10y 37.97840
# Gamma error distribution with a log link
m_ecm_5y <- gam(BAI_5y ~ ECM + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_5y <- summary(m_ecm_5y)$p.table["ECM", ]

m_ecm_10y <- gam(BAI_10y ~ ECM + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_10y <- summary(m_ecm_10y)$p.table["ECM", ]

m_ecm_15y <- gam(BAI_15y ~ ECM + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_15y <- summary(m_ecm_15y)$p.table["ECM", ]

m_ecm_30y <- gam(BAI_30y ~ ECM + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_30y <- summary(m_ecm_30y)$p.table["ECM", ]

m_sap_5y <- gam(BAI_5y ~ SAP + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_5y <- summary(m_sap_5y)$p.table["SAP", ]

m_sap_10y <- gam(BAI_10y ~ SAP + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_10y <- summary(m_sap_10y)$p.table["SAP", ]

m_sap_15y <- gam(BAI_15y ~ SAP + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_15y <- summary(m_sap_15y)$p.table["SAP", ]

m_sap_30y <- gam(BAI_30y ~ SAP + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_30y <- summary(m_sap_30y)$p.table["SAP", ]

res <- data.frame(
Guild = c("ECM", "ECM", "ECM", "ECM", "SAP", "SAP", "SAP", "SAP"),
Period = c("BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y", "BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y"),
Estimate = c(c_ecm_5y["Estimate"], c_ecm_10y["Estimate"], c_ecm_15y["Estimate"], c_ecm_30y["Estimate"], c_sap_5y["Estimate"], c_sap_10y["Estimate"], c_sap_15y["Estimate"], c_sap_30y["Estimate"]),
SE = c(c_ecm_5y["Std. Error"], c_ecm_10y["Std. Error"], c_ecm_15y["Std. Error"], c_ecm_30y["Std. Error"], c_sap_5y["Std. Error"], c_sap_10y["Std. Error"], c_sap_15y["Std. Error"], c_sap_30y["Std. Error"]),
p = c(c_ecm_5y["Pr(>|t|)"], c_ecm_10y["Pr(>|t|)"], c_ecm_15y["Pr(>|t|)"], c_ecm_30y["Pr(>|t|)"], c_sap_5y["Pr(>|t|)"], c_sap_10y["Pr(>|t|)"], c_sap_15y["Pr(>|t|)"], c_sap_30y["Pr(>|t|)"]),
R2_adj= c(summary(m_ecm_5y)$r.sq, summary(m_ecm_10y)$r.sq, summary(m_ecm_15y)$r.sq, summary(m_ecm_30y)$r.sq, summary(m_sap_5y)$r.sq, summary(m_sap_10y)$r.sq, summary(m_sap_15y)$r.sq, summary(m_sap_30y)$r.sq),
df= c(df.residual(m_ecm_5y), df.residual(m_ecm_10y), df.residual(m_ecm_15y), df.residual(m_ecm_30y), df.residual(m_sap_5y), df.residual(m_sap_10y), df.residual(m_sap_15y), df.residual(m_sap_30y)))

# BH correction for the 8 tests
res$p_adj <- p.adjust(res$p, method = "BH")   
res$Effect_10pct <- exp(10 * res$Estimate)   # change in BAI per 10 percentage points of guild abundance
res
##   Guild  Period     Estimate          SE         p    R2_adj  df     p_adj
## 1   ECM  BAI_5y  0.004133358 0.004058938 0.3109763 0.7905557 100 0.4621453
## 2   ECM BAI_10y  0.004595992 0.003890339 0.2402514 0.8006379 100 0.4621453
## 3   ECM BAI_15y  0.004492167 0.003701337 0.2277360 0.7868249 100 0.4621453
## 4   ECM BAI_30y  0.003268572 0.003042829 0.2853233 0.7716707 100 0.4621453
## 5   SAP  BAI_5y -0.003655450 0.004952102 0.4621453 0.7803871 100 0.4621453
## 6   SAP BAI_10y -0.003955726 0.004746442 0.4066004 0.7901313 100 0.4621453
## 7   SAP BAI_15y -0.003842385 0.004519367 0.3972425 0.7795364 100 0.4621453
## 8   SAP BAI_30y -0.003063235 0.003709709 0.4109203 0.7747531 100 0.4621453
##   Effect_10pct
## 1    1.0421997
## 2    1.0470324
## 3    1.0459459
## 4    1.0332258
## 5    0.9641056
## 6    0.9612149
## 7    0.9623050
## 8    0.9698321

Guild abundance vs detrended BAI

# detrended BAI is already size-free, so no log(DBH) needed
m_ecm_5y <- gam(detBAI_5y ~ ECM + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_5y <- summary(m_ecm_5y)$p.table["ECM", ]

m_ecm_10y <- gam(detBAI_10y ~ ECM + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_10y <- summary(m_ecm_10y)$p.table["ECM", ]

m_ecm_15y <- gam(detBAI_15y ~ ECM + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_15y <- summary(m_ecm_15y)$p.table["ECM", ]

m_ecm_30y <- gam(detBAI_30y ~ ECM + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_30y <- summary(m_ecm_30y)$p.table["ECM", ]

m_sap_5y <- gam(detBAI_5y ~ SAP + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_5y <- summary(m_sap_5y)$p.table["SAP", ]

m_sap_10y <- gam(detBAI_10y ~ SAP + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_10y <- summary(m_sap_10y)$p.table["SAP", ]

m_sap_15y <- gam(detBAI_15y ~ SAP + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_15y <- summary(m_sap_15y)$p.table["SAP", ]

m_sap_30y <- gam(detBAI_30y ~ SAP + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_30y <- summary(m_sap_30y)$p.table["SAP", ]

res <- data.frame(
Guild = c("ECM", "ECM", "ECM", "ECM", "SAP", "SAP", "SAP", "SAP"),
Period = c("detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y", "detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y"),
Estimate = c(c_ecm_5y["Estimate"], c_ecm_10y["Estimate"], c_ecm_15y["Estimate"], c_ecm_30y["Estimate"], c_sap_5y["Estimate"], c_sap_10y["Estimate"], c_sap_15y["Estimate"], c_sap_30y["Estimate"]),
SE = c(c_ecm_5y["Std. Error"], c_ecm_10y["Std. Error"], c_ecm_15y["Std. Error"], c_ecm_30y["Std. Error"], c_sap_5y["Std. Error"], c_sap_10y["Std. Error"], c_sap_15y["Std. Error"], c_sap_30y["Std. Error"]),
p  = c(c_ecm_5y["Pr(>|t|)"], c_ecm_10y["Pr(>|t|)"], c_ecm_15y["Pr(>|t|)"], c_ecm_30y["Pr(>|t|)"], c_sap_5y["Pr(>|t|)"], c_sap_10y["Pr(>|t|)"], c_sap_15y["Pr(>|t|)"], c_sap_30y["Pr(>|t|)"]),
R2_adj   = c(summary(m_ecm_5y)$r.sq, summary(m_ecm_10y)$r.sq, summary(m_ecm_15y)$r.sq, summary(m_ecm_30y)$r.sq, summary(m_sap_5y)$r.sq, summary(m_sap_10y)$r.sq, summary(m_sap_15y)$r.sq, summary(m_sap_30y)$r.sq),
df = c(df.residual(m_ecm_5y), df.residual(m_ecm_10y), df.residual(m_ecm_15y), df.residual(m_ecm_30y), df.residual(m_sap_5y), df.residual(m_sap_10y), df.residual(m_sap_15y), df.residual(m_sap_30y)))

# BH correction for the 8 tests
res$p_adj <- p.adjust(res$p, method = "BH")   
res$Effect_10pct <- exp(10 * res$Estimate)   # change in BAI per 10 percentage points of guild abundance
res
##   Guild     Period      Estimate           SE         p    R2_adj  df     p_adj
## 1   ECM  detBAI_5y -1.400317e-04 0.0006027762 0.8167664 0.2283103 101 0.8167664
## 2   ECM detBAI_10y  3.198858e-04 0.0003311044 0.3362931 0.6879030 101 0.6591815
## 3   ECM detBAI_15y  6.522088e-05 0.0002037305 0.7495288 0.3714417 101 0.8167664
## 4   ECM detBAI_30y  1.492564e-04 0.0001726922 0.3894758 0.3166561 101 0.6591815
## 5   SAP  detBAI_5y -6.149961e-04 0.0007465263 0.4119884 0.2332570 101 0.6591815
## 6   SAP detBAI_10y -4.219750e-04 0.0004122617 0.3084880 0.6881250 101 0.6591815
## 7   SAP detBAI_15y -8.658286e-05 0.0002536714 0.7335726 0.3712627 101 0.8167664
## 8   SAP detBAI_30y -2.733358e-04 0.0002140979 0.2046407 0.3231283 101 0.6591815
##   Effect_10pct
## 1    0.9986007
## 2    1.0032040
## 3    1.0006524
## 4    1.0014937
## 5    0.9938689
## 6    0.9957891
## 7    0.9991345
## 8    0.9972704

Guild abundance vs RCS BAI

# RCS BAI: growth relative to trees of the same cambial age in the same plot (no log(DBH)).
# Each plot has its own regional curve, so the index level differs between plots:
# Plot (= Site x Habitat) is used instead of Habitat + Site

#Check for multicollinearity
vif(lm(log(rcsBAI_10y) ~ ECM + Plot, data = guild_bai))
##          GVIF Df GVIF^(1/(2*Df))
## ECM  1.269218  1        1.126596
## Plot 1.269218  5        1.024127
vif(lm(log(rcsBAI_10y) ~ SAP + Plot, data = guild_bai))
##          GVIF Df GVIF^(1/(2*Df))
## SAP  1.326615  1        1.151788
## Plot 1.326615  5        1.028666
m_ecm_5y <- gam(rcsBAI_5y ~ ECM + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_5y <- summary(m_ecm_5y)$p.table["ECM", ]

m_ecm_10y <- gam(rcsBAI_10y ~ ECM + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_10y <- summary(m_ecm_10y)$p.table["ECM", ]

m_ecm_15y <- gam(rcsBAI_15y ~ ECM + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_15y <- summary(m_ecm_15y)$p.table["ECM", ]

m_ecm_30y <- gam(rcsBAI_30y ~ ECM + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_ecm_30y <- summary(m_ecm_30y)$p.table["ECM", ]

m_sap_5y <- gam(rcsBAI_5y ~ SAP + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_5y <- summary(m_sap_5y)$p.table["SAP", ]

m_sap_10y <- gam(rcsBAI_10y ~ SAP + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_10y <- summary(m_sap_10y)$p.table["SAP", ]

m_sap_15y <- gam(rcsBAI_15y ~ SAP + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_15y <- summary(m_sap_15y)$p.table["SAP", ]

m_sap_30y <- gam(rcsBAI_30y ~ SAP + Plot, family = Gamma(link = "log"), method = "REML", data = guild_bai)
c_sap_30y <- summary(m_sap_30y)$p.table["SAP", ]

res <- data.frame(
Guild = c("ECM", "ECM", "ECM", "ECM", "SAP", "SAP", "SAP", "SAP"),
Period = c("rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y", "rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y"),
Estimate = c(c_ecm_5y["Estimate"], c_ecm_10y["Estimate"], c_ecm_15y["Estimate"], c_ecm_30y["Estimate"], c_sap_5y["Estimate"], c_sap_10y["Estimate"], c_sap_15y["Estimate"], c_sap_30y["Estimate"]),
SE  = c(c_ecm_5y["Std. Error"], c_ecm_10y["Std. Error"], c_ecm_15y["Std. Error"], c_ecm_30y["Std. Error"], c_sap_5y["Std. Error"], c_sap_10y["Std. Error"], c_sap_15y["Std. Error"], c_sap_30y["Std. Error"]),
p = c(c_ecm_5y["Pr(>|t|)"], c_ecm_10y["Pr(>|t|)"], c_ecm_15y["Pr(>|t|)"], c_ecm_30y["Pr(>|t|)"], c_sap_5y["Pr(>|t|)"], c_sap_10y["Pr(>|t|)"], c_sap_15y["Pr(>|t|)"], c_sap_30y["Pr(>|t|)"]),
R2_adj = c(summary(m_ecm_5y)$r.sq, summary(m_ecm_10y)$r.sq, summary(m_ecm_15y)$r.sq, summary(m_ecm_30y)$r.sq, summary(m_sap_5y)$r.sq, summary(m_sap_10y)$r.sq, summary(m_sap_15y)$r.sq, summary(m_sap_30y)$r.sq),
df = c(df.residual(m_ecm_5y), df.residual(m_ecm_10y), df.residual(m_ecm_15y), df.residual(m_ecm_30y), df.residual(m_sap_5y), df.residual(m_sap_10y), df.residual(m_sap_15y), df.residual(m_sap_30y)))

# BH correction for the 8 tests
res$p_adj <- p.adjust(res$p, method = "BH")   
res$Effect_10pct <- exp(10 * res$Estimate)   # change in BAI per 10 percentage points of guild abundance
res
##   Guild     Period      Estimate          SE         p    R2_adj df     p_adj
## 1   ECM  rcsBAI_5y -0.0026286614 0.005476495 0.6322935 0.2996420 99 0.9122964
## 2   ECM rcsBAI_10y -0.0020507612 0.005252767 0.6970682 0.2643650 99 0.9122964
## 3   ECM rcsBAI_15y -0.0016322076 0.004997894 0.7446757 0.2198988 99 0.9122964
## 4   ECM rcsBAI_30y -0.0019151393 0.004563198 0.6756196 0.1645476 99 0.9122964
## 5   SAP  rcsBAI_5y  0.0030183171 0.006946127 0.6648475 0.2947995 99 0.9122964
## 6   SAP rcsBAI_10y  0.0020979103 0.006661072 0.7534618 0.2596935 99 0.9122964
## 7   SAP rcsBAI_15y  0.0013469161 0.006333303 0.8320201 0.2157366 99 0.9122964
## 8   SAP rcsBAI_30y  0.0006375491 0.005773630 0.9122964 0.1588023 99 0.9122964
##   Effect_10pct
## 1    0.9740559
## 2    0.9797012
## 3    0.9838104
## 4    0.9810308
## 5    1.0306433
## 6    1.0212007
## 7    1.0135603
## 8    1.0063959