library("ggplot2")
library("readr")
library("phyloseq")
library("dplyr")
library("tidyr")
library('tibble')
library("rstatix")
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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] mgcv_1.9-4      nlme_3.1-168    car_3.1-5       carData_3.0-6  
##  [5] vegan_2.7-2     permute_0.9-10  rstatix_1.1.0   tibble_3.3.1   
##  [9] tidyr_1.3.2     dplyr_1.2.1     phyloseq_1.52.0 readr_2.1.6    
## [13] ggplot2_4.0.2  
## 
## loaded via a namespace (and not attached):
##  [1] ade4_1.7-23             tidyselect_1.2.1        farver_2.1.2           
##  [4] Biostrings_2.76.0       S7_0.2.1                fastmap_1.2.0          
##  [7] digest_0.6.37           lifecycle_1.0.5         cluster_2.1.8.2        
## [10] survival_3.8-6          magrittr_2.0.3          compiler_4.5.0         
## [13] rlang_1.2.0             sass_0.4.10             tools_4.5.0            
## [16] igraph_2.2.2            yaml_2.3.10             data.table_1.18.2.1    
## [19] knitr_1.51              plyr_1.8.9              RColorBrewer_1.1-3     
## [22] abind_1.4-8             withr_3.0.2             purrr_1.2.1            
## [25] BiocGenerics_0.54.1     grid_4.5.0              stats4_4.5.0           
## [28] multtest_2.64.0         biomformat_1.36.0       colorspace_2.1-2       
## [31] Rhdf5lib_1.30.0         scales_1.4.0            iterators_1.0.14       
## [34] MASS_7.3-65             dichromat_2.0-0.1       cli_3.6.5              
## [37] rmarkdown_2.30          crayon_1.5.3            generics_0.1.4         
## [40] otel_0.2.0              rstudioapi_0.18.0       httr_1.4.8             
## [43] reshape2_1.4.5          tzdb_0.5.0              ape_5.8-1              
## [46] cachem_1.1.0            rhdf5_2.52.1            stringr_1.6.0          
## [49] splines_4.5.0           parallel_4.5.0          XVector_0.48.0         
## [52] vctrs_0.7.3             Matrix_1.7-4            jsonlite_2.0.0         
## [55] IRanges_2.42.0          hms_1.1.4               S4Vectors_0.46.0       
## [58] Formula_1.2-5           foreach_1.5.2           jquerylib_0.1.4        
## [61] glue_1.8.0              codetools_0.2-20        stringi_1.8.7          
## [64] gtable_0.3.6            GenomeInfoDb_1.44.3     UCSC.utils_1.4.0       
## [67] pillar_1.11.1           htmltools_0.5.8.1       rhdf5filters_1.20.0    
## [70] GenomeInfoDbData_1.2.14 R6_2.6.1                evaluate_1.0.5         
## [73] lattice_0.22-9          Biobase_2.68.0          backports_1.5.0        
## [76] broom_1.0.12            bslib_0.10.0            Rcpp_1.1.0             
## [79] xfun_0.61               pkgconfig_2.0.3

Read data

# output of BAI_detBAI_RCS.Rmd: raw BAI, detrended BAI (detBAI) and RCS BAI (rcsBAI)
sample_info_no_R <-read.table("Alaska_info_no_R.txt", row.names=1, header=T, check.names=F)
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, ]

NTC treatment

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
length(blank_OTUs)

## [1] 4

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

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) ]

All OTUs

Data preparation

meta <- as(sample_info_no_R, "data.frame")
keep <- complete.cases(meta[, c("BAI_5y","Habitat","Site")])
meta_clean <- meta[keep, ]

# OTU table + diversity (raw counts)
otu <- otu_table(filt_count_tab_no_r, taxa_are_rows = TRUE)
physeq_tmp <- phyloseq(otu)

div_df <- estimate_richness(physeq_tmp, measures = c("Shannon", "Observed", "InvSimpson"))
rownames(div_df) <- sub("^X", "", rownames(div_df))
rownames(div_df) <- gsub("\\.", "-", rownames(div_df))
div_df$Sample <- rownames(div_df)

meta_clean$Sample <- rownames(meta_clean)

common <- intersect(div_df$Sample, meta_clean$Sample)
div_only <- setdiff(div_df$Sample, meta_clean$Sample)
div_df <- div_df[div_df$Sample %in% common, ]
meta2   <- meta_clean[meta_clean$Sample %in% common, ]
div_df2 <- merge(div_df, meta2, by = "Sample")

div_df2_clean <- div_df2[is.finite(div_df2$BAI_5y), ] # remove samples with NA

Model check: Gaussian vs Gamma (Shannon, BAI 5 years)

#Check for multicollinearity (should be less than 1.1)
vif(lm(log(BAI_5y) ~ Observed + Habitat + Site + log(DBH), data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Observed 1.005313  1        1.002653
## Habitat  1.028227  1        1.014015
## Site     1.040554  2        1.009988
## log(DBH) 1.046153  1        1.022816

vif(lm(log(BAI_5y) ~ Shannon + Habitat + Site + log(DBH), data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Shannon  1.026227  1        1.013029
## Habitat  1.046458  1        1.022965
## Site     1.041340  2        1.010179
## log(DBH) 1.045713  1        1.022601

vif(lm(log(BAI_5y) ~ InvSimpson + Habitat + Site + log(DBH), data = div_df2))

##                GVIF Df GVIF^(1/(2*Df))
## InvSimpson 1.026054  1        1.012943
## Habitat    1.035998  1        1.017840
## Site       1.053360  2        1.013081
## log(DBH)   1.047869  1        1.023655

# same model, two error distributions
m_gaus <- gam(BAI_5y ~ Shannon + Habitat + Site + log(DBH),
  family = gaussian(link = "log"),
  method = "REML",
  data = div_df2_clean)

m_gamma <- gam(BAI_5y ~ Shannon + Habitat + Site + log(DBH),
  family = Gamma(link = "log"),
  method = "REML",
  data = div_df2_clean)

AIC(m_gaus, m_gamma)  # lower = better

##         df      AIC
## m_gaus   7 652.1808
## m_gamma  7 564.3223

gam.check(m_gaus)     # "Resids vs. linear pred." fans out = variance grows with BAI

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-2.060658e-07,-2.060658e-07]
## (score 327.6589 & scale 25.55837).
## 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.0002062811,-0.0002062811]
## (score 282.6833 & scale 0.3686196).
## Hessian positive definite, eigenvalue range [56.2358,56.2358].
## Model rank =  6 / 6

summary(m_gamma)

## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## BAI_5y ~ Shannon + Habitat + Site + log(DBH)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.212916   0.525624  -4.210 5.59e-05 ***
## Shannon          -0.008972   0.127928  -0.070    0.944    
## Habitattreeline   0.491472   0.120995   4.062 9.69e-05 ***
## SiteBrooks_Range -0.120930   0.140410  -0.861    0.391    
## SiteInterior      1.164111   0.149984   7.762 7.37e-12 ***
## log(DBH)          1.257818   0.136637   9.206 5.46e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.776   Deviance explained =   69%
## -REML = 282.68  Scale est. = 0.36862   n = 106

# prediction line at the median DBH of each plot, so it shows the Shannon effect only
nd <- expand.grid(Shannon = seq(min(div_df2_clean$Shannon), max(div_df2_clean$Shannon), length.out = 50),
                  Habitat = unique(div_df2_clean$Habitat),
                  Site = unique(div_df2_clean$Site))
med_dbh <- aggregate(DBH ~ Habitat + Site, data = div_df2_clean, FUN = median)
nd <- merge(nd, med_dbh, by = c("Habitat", "Site"))
nd$pred <- predict(m_gamma, newdata = nd, type = "response")

ggplot(div_df2_clean, aes(x = Shannon, y = BAI_5y)) + 
  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() + labs(x = "Shannon index", y = "Mean BAI over 5 years (cm²/year, log scale)") +
  theme_bw()

RAF diversity 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, for one model per
## growth measure, instead of refitting a Gaussian model beside every model below.

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

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

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

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 71.88924
## 2 detBAI_10y 14.04768
## 3 rcsBAI_10y 37.01982

# Gamma error distribution with a log link (see the family check above)
m_obs_5y <- gam(BAI_5y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(BAI_10y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(BAI_15y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(BAI_30y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(BAI_5y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(BAI_10y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(BAI_15y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(BAI_30y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(BAI_5y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(BAI_10y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(BAI_15y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(BAI_30y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y", "BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y", "BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE  = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj = c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

# BH correction for the 12 tests
res$p_adj <- p.adjust(res$p, method = "BH")   

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST5_All_OTU_BAI.csv", row.names = FALSE)
res

##         Index  Period      Estimate          SE         p    R2_adj  df
## 1    Observed  BAI_5y  0.0010550441 0.001594889 0.5098030 0.7837078 100
## 2    Observed BAI_10y  0.0010115763 0.001527753 0.5094088 0.7935174 100
## 3    Observed BAI_15y  0.0009943306 0.001454333 0.4957444 0.7828433 100
## 4    Observed BAI_30y  0.0010639781 0.001192821 0.3745424 0.7774324 100
## 5     Shannon  BAI_5y -0.0089718891 0.127928448 0.9442286 0.7756628 100
## 6     Shannon BAI_10y -0.0072282573 0.122846243 0.9531971 0.7846923 100
## 7     Shannon BAI_15y -0.0025383921 0.117025170 0.9827377 0.7753536 100
## 8     Shannon BAI_30y  0.0299982033 0.096125951 0.7556357 0.7745415 100
## 9  InvSimpson  BAI_5y -0.0061680741 0.014113365 0.6630262 0.7718928 100
## 10 InvSimpson BAI_10y -0.0056914247 0.013565598 0.6757154 0.7820253 100
## 11 InvSimpson BAI_15y -0.0053783876 0.012929555 0.6783187 0.7732048 100
## 12 InvSimpson BAI_30y -0.0010737663 0.010637913 0.9198020 0.7732729 100
##        p_adj       CI_low     CI_high   SD_index per_SD_pct
## 1  0.9827377 -0.002109170 0.004219258 37.3271602  4.0167543
## 2  0.9827377 -0.002019442 0.004042595 37.3271602  3.8481209
## 3  0.9827377 -0.001891024 0.003879685 37.3271602  3.7812919
## 4  0.9827377 -0.001302544 0.003430501 37.3271602  4.0514478
## 5  0.9827377 -0.262778286 0.244834508  0.4691905 -0.4200677
## 6  0.9827377 -0.250951705 0.236495190  0.4691905 -0.3385685
## 7  0.9827377 -0.234712997 0.229636213  0.4691905 -0.1190280
## 8  0.9827377 -0.160712946 0.220709352  0.4691905  1.4174389
## 9  0.9827377 -0.034168588 0.021832440  4.2451399 -2.5844500
## 10 0.9827377 -0.032605185 0.021222335  4.2451399 -2.3871356
## 11 0.9827377 -0.031030256 0.020273480  4.2451399 -2.2573330
## 12 0.9827377 -0.022179082 0.020031549  4.2451399 -0.4547915

RAF diversity vs detrended BAI

# detrended BAI is already size-free, so no log(DBH) term
m_obs_5y <- gam(detBAI_5y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(detBAI_10y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(detBAI_15y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(detBAI_30y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(detBAI_5y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(detBAI_10y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(detBAI_15y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(detBAI_30y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(detBAI_5y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(detBAI_10y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(detBAI_15y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(detBAI_30y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y", "detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y", "detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj = c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

# BH correction for the 12 tests
res$p_adj <- p.adjust(res$p, method = "BH")   

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST5_All_OTU_detBAI.csv", row.names = FALSE)
res

##         Index     Period      Estimate           SE         p    R2_adj  df
## 1    Observed  detBAI_5y  3.684445e-04 2.436041e-04 0.1335364 0.2441914 101
## 2    Observed detBAI_10y -1.219105e-05 1.360509e-04 0.9287773 0.6836360 101
## 3    Observed detBAI_15y -9.744599e-05 8.273982e-05 0.2416680 0.3783426 101
## 4    Observed detBAI_30y  9.940769e-05 7.023754e-05 0.1600545 0.3242879 101
## 5     Shannon  detBAI_5y  1.982363e-02 1.973390e-02 0.3175159 0.2346485 101
## 6     Shannon detBAI_10y  3.834749e-03 1.093402e-02 0.7265309 0.6836536 101
## 7     Shannon detBAI_15y -3.680093e-03 6.692025e-03 0.5835868 0.3727836 101
## 8     Shannon detBAI_30y  6.223525e-03 5.666559e-03 0.2746892 0.3192438 101
## 9  InvSimpson  detBAI_5y  1.964056e-03 2.180860e-03 0.3699498 0.2335490 101
## 10 InvSimpson detBAI_10y  8.278421e-04 1.205086e-03 0.4936841 0.6841933 101
## 11 InvSimpson detBAI_15y -3.234262e-04 7.392661e-04 0.6626855 0.3720464 101
## 12 InvSimpson detBAI_30y  5.730547e-04 6.268062e-04 0.3627653 0.3168943 101
##        p_adj        CI_low      CI_high   SD_index  per_SD_pct
## 1  0.6341996 -0.0001148004 8.516895e-04 37.3271602  1.38479957
## 2  0.9287773 -0.0002820794 2.576973e-04 37.3271602 -0.04549537
## 3  0.6341996 -0.0002615795 6.668756e-05 37.3271602 -0.36307749
## 4  0.6341996 -0.0000399247 2.387401e-04 37.3271602  0.37174997
## 5  0.6341996 -0.0193231304 5.897038e-02  0.4691905  0.93444458
## 6  0.7925792 -0.0178554036 2.552490e-02  0.4691905  0.18008474
## 7  0.7781157 -0.0169552696 9.595084e-03  0.4691905 -0.17251748
## 8  0.6341996 -0.0050174038 1.746445e-02  0.4691905  0.29242861
## 9  0.6341996 -0.0023621830 6.290296e-03  4.2451399  0.83725497
## 10 0.7405261 -0.0015627237 3.218408e-03  4.2451399  0.35204877
## 11 0.7925792 -0.0017899313 1.143079e-03  4.2451399 -0.13720474
## 12 0.6341996 -0.0006703602 1.816470e-03  4.2451399  0.24356587

RAF diversity vs RCS BAI

# RCS BAI: growth relative to trees of the same cambial age in the same plot (no log(DBH) term).
# 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_5y) ~ Observed + Plot, data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Observed 1.035535  1        1.017613
## Plot     1.035535  5        1.003498

vif(lm(log(rcsBAI_5y) ~ Shannon + Plot, data = div_df2))

##             GVIF Df GVIF^(1/(2*Df))
## Shannon 1.027569  1        1.013691
## Plot    1.027569  5        1.002723

vif(lm(log(rcsBAI_5y) ~ InvSimpson + Plot, data = div_df2))

##                GVIF Df GVIF^(1/(2*Df))
## InvSimpson 1.027238  1        1.013527
## Plot       1.027238  5        1.002691

m_obs_5y <- gam(rcsBAI_5y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(rcsBAI_10y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(rcsBAI_15y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(rcsBAI_30y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(rcsBAI_5y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(rcsBAI_10y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(rcsBAI_15y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(rcsBAI_30y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(rcsBAI_5y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(rcsBAI_10y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(rcsBAI_15y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(rcsBAI_30y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y", "rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y", "rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj = c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

# BH correction for the 12 tests
res$p_adj <- p.adjust(res$p, method = "BH")   

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST5_All_OTU_rcsBAI.csv", row.names = FALSE)
res

##         Index     Period    Estimate          SE         p    R2_adj df
## 1    Observed  rcsBAI_5y 0.001851975 0.002234664 0.4092401 0.3090068 99
## 2    Observed rcsBAI_10y 0.001758124 0.002133209 0.4118250 0.2737775 99
## 3    Observed rcsBAI_15y 0.001784082 0.002023208 0.3800166 0.2288748 99
## 4    Observed rcsBAI_30y 0.002229375 0.001829915 0.2260074 0.1697822 99
## 5     Shannon  rcsBAI_5y 0.112439010 0.175843497 0.5240223 0.2963190 99
## 6     Shannon rcsBAI_10y 0.109415647 0.168323821 0.5171780 0.2598436 99
## 7     Shannon rcsBAI_15y 0.109576621 0.159924061 0.4948321 0.2143226 99
## 8     Shannon rcsBAI_30y 0.155393830 0.144641470 0.2852827 0.1560446 99
## 9  InvSimpson  rcsBAI_5y 0.001620081 0.019510673 0.9339909 0.2916794 99
## 10 InvSimpson rcsBAI_10y 0.002845219 0.018680487 0.8792527 0.2566024 99
## 11 InvSimpson rcsBAI_15y 0.003813671 0.017744326 0.8302692 0.2121456 99
## 12 InvSimpson rcsBAI_30y 0.009553971 0.016108508 0.5544661 0.1505328 99
##        p_adj       CI_low     CI_high   SD_index per_SD_pct
## 1  0.7392881 -0.002582082 0.006286033 37.3271602  7.1574412
## 2  0.7392881 -0.002474625 0.005990873 37.3271602  6.7827028
## 3  0.7392881 -0.002230403 0.005798566 37.3271602  6.8862178
## 4  0.7392881 -0.001401572 0.005860322 37.3271602  8.6776784
## 5  0.7392881 -0.236472637 0.461350658  0.4691905  5.4171673
## 6  0.7392881 -0.224575331 0.443406626  0.4691905  5.2677356
## 7  0.7392881 -0.207747411 0.426900654  0.4691905  5.2756865
## 8  0.7392881 -0.131606228 0.442393887  0.4691905  7.5632981
## 9  0.9339909 -0.037093326 0.040333489  4.2451399  0.6901175
## 10 0.9339909 -0.034220920 0.039911357  4.2451399  1.2151590
## 11 0.9339909 -0.031394921 0.039022263  4.2451399  1.6321328
## 12 0.7392881 -0.022408803 0.041516745  4.2451399  4.1391650

RAF diversity vs DBH and age

# DBH and age as response: tests whether tree size/age differ with RAF diversity (association only)
m_obs_dbh <- gam(DBH ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_dbh <- summary(m_obs_dbh)$p.table["Observed", ]

m_obs_age <- gam(Age ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_age <- summary(m_obs_age)$p.table["Observed", ]

m_sha_dbh <- gam(DBH ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_dbh <- summary(m_sha_dbh)$p.table["Shannon", ]

m_sha_age <- gam(Age ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_age <- summary(m_sha_age)$p.table["Shannon", ]

m_inv_dbh <- gam(DBH ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_dbh <- summary(m_inv_dbh)$p.table["InvSimpson", ]

m_inv_age <- gam(Age ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_age <- summary(m_inv_age)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Shannon", "Shannon", "InvSimpson", "InvSimpson"),
Period = c("DBH", "Age", "DBH", "Age", "DBH", "Age"),
Estimate = c(c_obs_dbh["Estimate"], c_obs_age["Estimate"], c_sha_dbh["Estimate"], c_sha_age["Estimate"], c_inv_dbh["Estimate"], c_inv_age["Estimate"]),
SE = c(c_obs_dbh["Std. Error"], c_obs_age["Std. Error"], c_sha_dbh["Std. Error"], c_sha_age["Std. Error"], c_inv_dbh["Std. Error"], c_inv_age["Std. Error"]),
p  = c(c_obs_dbh["Pr(>|t|)"], c_obs_age["Pr(>|t|)"], c_sha_dbh["Pr(>|t|)"], c_sha_age["Pr(>|t|)"], c_inv_dbh["Pr(>|t|)"], c_inv_age["Pr(>|t|)"]),
R2_adj = c(summary(m_obs_dbh)$r.sq, summary(m_obs_age)$r.sq, summary(m_sha_dbh)$r.sq, summary(m_sha_age)$r.sq, summary(m_inv_dbh)$r.sq, summary(m_inv_age)$r.sq),
df = c(df.residual(m_obs_dbh), df.residual(m_obs_age), df.residual(m_sha_dbh), df.residual(m_sha_age), df.residual(m_inv_dbh), df.residual(m_inv_age)))

# BH correction for the 6 tests
res$p_adj <- p.adjust(res$p, method = "BH")  
res

##        Index Period      Estimate          SE          p     R2_adj  df
## 1   Observed    DBH -0.0007865061 0.001142883 0.49291860 0.02835140 101
## 2   Observed    Age -0.0021546318 0.001086866 0.05014811 0.23364854 101
## 3    Shannon    DBH  0.0600914034 0.092083020 0.51551058 0.02497674 101
## 4    Shannon    Age -0.0537757832 0.088440554 0.54452292 0.21912769 101
## 5 InvSimpson    DBH  0.0069183682 0.010175789 0.49813254 0.02548133 101
## 6 InvSimpson    Age  0.0004162263 0.009791586 0.96617721 0.21985945 101
##       p_adj
## 1 0.6534275
## 2 0.3008887
## 3 0.6534275
## 4 0.6534275
## 5 0.6534275
## 6 0.9661772

ECM OTUs

Data preparation (ECM)

# OTU table + diversity
# Raw count phyloseq object - select ECM
otu <- otu_table(filt_count_tab_no_r, taxa_are_rows = TRUE)
otu_ECM <- tax_tab[tax_tab$Guild %in% c("Ectomycorrhizal"),]
ECM_OTUs_present <- intersect(otu_ECM$OTU, rownames(filt_count_tab_no_r))

# create phyloseq object with raw ECM counts
otu_phy_ECM <- otu_table(filt_count_tab_no_r[ECM_OTUs_present,], taxa_are_rows=T)
otu_phy_ECM_data <- sample_data(sample_info_no_R)
phy_ECM <- phyloseq(otu_phy_ECM, otu_phy_ECM_data)

div_df <- estimate_richness(phy_ECM, measures = c("Shannon", "Observed", "InvSimpson"))
rownames(div_df) <- sub("^X", "", rownames(div_df))
rownames(div_df) <- gsub("\\.", "-", rownames(div_df))
div_df$Sample <- rownames(div_df)

meta_clean$Sample <- rownames(meta_clean)

common <- intersect(div_df$Sample, meta_clean$Sample)
div_only <- setdiff(div_df$Sample, meta_clean$Sample)
div_df <- div_df[div_df$Sample %in% common, ]
meta2   <- meta_clean[meta_clean$Sample %in% common, ]
div_df2 <- merge(div_df, meta2, by = "Sample")

div_df2_clean <- div_df2[is.finite(div_df2$BAI_5y), ] # remove samples with NA

RAF diversity vs raw BAI (ECM)

#Check for multicollinearity
vif(lm(log(BAI_5y) ~ Observed + Habitat + Site + log(DBH), data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Observed 1.073625  1        1.036159
## Habitat  1.030414  1        1.015093
## Site     1.110602  2        1.026573
## log(DBH) 1.042862  1        1.021206

vif(lm(log(BAI_5y) ~ Shannon + Habitat + Site + log(DBH), data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Shannon  1.127566  1        1.061869
## Habitat  1.033749  1        1.016735
## Site     1.164885  2        1.038893
## log(DBH) 1.041368  1        1.020475

vif(lm(log(BAI_5y) ~ InvSimpson + Habitat + Site + log(DBH), data = div_df2))

##                GVIF Df GVIF^(1/(2*Df))
## InvSimpson 1.107889  1        1.052563
## Habitat    1.028180  1        1.013992
## Site       1.149781  2        1.035509
## log(DBH)   1.041178  1        1.020381

m_obs_5y <- gam(BAI_5y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(BAI_10y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(BAI_15y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(BAI_30y ~ Observed + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(BAI_5y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(BAI_10y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(BAI_15y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(BAI_30y ~ Shannon + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(BAI_5y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(BAI_10y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(BAI_15y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(BAI_30y ~ InvSimpson + Habitat + Site + log(DBH), family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y", "BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y", "BAI_5y", "BAI_10y", "BAI_15y", "BAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj = c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

# BH correction for the 12 tests
res$p_adj <- p.adjust(res$p, method = "BH")  

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST5_ECM_OTU_BAI.csv", row.names = FALSE)
res

##         Index  Period      Estimate          SE         p    R2_adj  df
## 1    Observed  BAI_5y  0.0006968489 0.003178169 0.8268935 0.7798278 100
## 2    Observed BAI_10y  0.0007479993 0.003049188 0.8067184 0.7894626 100
## 3    Observed BAI_15y  0.0006895926 0.002904075 0.8127882 0.7790210 100
## 4    Observed BAI_30y  0.0009775249 0.002385423 0.6828358 0.7754553 100
## 5     Shannon  BAI_5y -0.1125461930 0.127377961 0.3790525 0.7657152 100
## 6     Shannon BAI_10y -0.1031942088 0.122650219 0.4021482 0.7751166 100
## 7     Shannon BAI_15y -0.0917943699 0.117025304 0.4346610 0.7676382 100
## 8     Shannon BAI_30y -0.0426459188 0.096551018 0.6596650 0.7715419 100
## 9  InvSimpson  BAI_5y -0.0325433654 0.023770404 0.1740444 0.7687444 100
## 10 InvSimpson BAI_10y -0.0326082601 0.022899322 0.1575641 0.7780343 100
## 11 InvSimpson BAI_15y -0.0316576518 0.021870961 0.1508919 0.7704043 100
## 12 InvSimpson BAI_30y -0.0200934299 0.018116636 0.2700409 0.7733179 100
##        p_adj       CI_low     CI_high   SD_index per_SD_pct
## 1  0.8268935 -0.005608547 0.007002245 19.3371096   1.356624
## 2  0.8268935 -0.005301504 0.006797502 19.3371096   1.456926
## 3  0.8268935 -0.005072009 0.006451194 19.3371096   1.342403
## 4  0.8268935 -0.003755086 0.005710136 19.3371096   1.908229
## 5  0.7451332 -0.365260440 0.140168054  0.4904209  -5.369941
## 6  0.7451332 -0.346528749 0.140140332  0.4904209  -4.934932
## 7  0.7451332 -0.323969240 0.140380501  0.4904209  -4.401961
## 8  0.8268935 -0.234200388 0.148908550  0.4904209  -2.069726
## 9  0.6961776 -0.079703170 0.014616439  2.5836409  -8.064264
## 10 0.6961776 -0.078039863 0.012823343  2.5836409  -8.079677
## 11 0.6961776 -0.075049015 0.011733711  2.5836409  -7.853640
## 12 0.7451332 -0.056036320 0.015849460  2.5836409  -5.058968

RAF diversity vs detrended BAI (ECM)

m_obs_5y <- gam(detBAI_5y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(detBAI_10y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(detBAI_15y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(detBAI_30y ~ Observed + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(detBAI_5y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(detBAI_10y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(detBAI_15y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(detBAI_30y ~ Shannon + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(detBAI_5y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(detBAI_10y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(detBAI_15y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(detBAI_30y ~ InvSimpson + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y", "detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y", "detBAI_5y", "detBAI_10y", "detBAI_15y", "detBAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p  = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj  = c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

res$p_adj <- p.adjust(res$p, method = "BH")   # BH correction for the 12 tests

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST2_ECM_OTU_detBAI.csv", row.names = FALSE)
res

##         Index     Period      Estimate           SE          p    R2_adj  df
## 1    Observed  detBAI_5y  6.440308e-04 0.0004878658 0.18978618 0.2405466 101
## 2    Observed detBAI_10y  5.193072e-05 0.0002718717 0.84889972 0.6837198 101
## 3    Observed detBAI_15y -1.445743e-04 0.0001658867 0.38553472 0.3748175 101
## 4    Observed detBAI_30y  1.976470e-04 0.0001403477 0.16212290 0.3243674 101
## 5     Shannon  detBAI_5y  2.398509e-02 0.0197611119 0.22767303 0.2392999 101
## 6     Shannon detBAI_10y  5.882688e-03 0.0109809999 0.59333318 0.6843740 101
## 7     Shannon detBAI_15y -2.443398e-03 0.0067316466 0.71738449 0.3715468 101
## 8     Shannon detBAI_30y  9.188650e-03 0.0056537536 0.10723012 0.3292083 101
## 9  InvSimpson  detBAI_5y  6.161071e-03 0.0036913509 0.09820328 0.2509308 101
## 10 InvSimpson detBAI_10y  2.415666e-03 0.0020552626 0.24261557 0.6872479 101
## 11 InvSimpson detBAI_15y -3.215264e-04 0.0012672751 0.80023171 0.3710397 101
## 12 InvSimpson detBAI_30y  1.764196e-03 0.0010632451 0.10016640 0.3303066 101
##        p_adj        CI_low      CI_high   SD_index  per_SD_pct
## 1  0.4159124 -3.237636e-04 0.0016118252 19.3371096  1.25315648
## 2  0.8488997 -4.873895e-04 0.0005912510 19.3371096  0.10046944
## 3  0.5783021 -4.736488e-04 0.0001845002 19.3371096 -0.27917452
## 4  0.4159124 -8.076508e-05 0.0004760590 19.3371096  0.38292339
## 5  0.4159124 -1.521564e-02 0.0631858229  0.4904209  1.18322447
## 6  0.7911109 -1.590066e-02 0.0276660381  0.4904209  0.28891589
## 7  0.8488997 -1.579717e-02 0.0109103783  0.4904209 -0.11975757
## 8  0.4159124 -2.026876e-03 0.0204041764  0.4904209  0.45164749
## 9  0.4159124 -1.161576e-03 0.0134837187  2.5836409  1.60453626
## 10 0.4159124 -1.661422e-03 0.0064927537  2.5836409  0.62607295
## 11 0.8488997 -2.835459e-03 0.0021924065  2.5836409 -0.08303638
## 12 0.4159124 -3.449964e-04 0.0038733881  2.5836409  0.45684523

RAF diversity vs rcsBAI (ECM)

# RCS BAI: growth relative to trees of the same cambial age in the same plot (no log(DBH) term).
# 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_5y) ~ Observed + Plot, data = div_df2))

##              GVIF Df GVIF^(1/(2*Df))
## Observed 1.130459  1        1.063230
## Plot     1.130459  5        1.012338

vif(lm(log(rcsBAI_5y) ~ Shannon + Plot, data = div_df2))

##             GVIF Df GVIF^(1/(2*Df))
## Shannon 1.172308  1        1.082732
## Plot    1.172308  5        1.016024

vif(lm(log(rcsBAI_5y) ~ InvSimpson + Plot, data = div_df2))

##                GVIF Df GVIF^(1/(2*Df))
## InvSimpson 1.138875  1        1.067181
## Plot       1.138875  5        1.013089

m_obs_5y <- gam(rcsBAI_5y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_5y <- summary(m_obs_5y)$p.table["Observed", ]

m_obs_10y <- gam(rcsBAI_10y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_10y <- summary(m_obs_10y)$p.table["Observed", ]

m_obs_15y <- gam(rcsBAI_15y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_15y <- summary(m_obs_15y)$p.table["Observed", ]

m_obs_30y <- gam(rcsBAI_30y ~ Observed + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_obs_30y <- summary(m_obs_30y)$p.table["Observed", ]

m_sha_5y <- gam(rcsBAI_5y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_5y <- summary(m_sha_5y)$p.table["Shannon", ]

m_sha_10y <- gam(rcsBAI_10y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_10y <- summary(m_sha_10y)$p.table["Shannon", ]

m_sha_15y <- gam(rcsBAI_15y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_15y <- summary(m_sha_15y)$p.table["Shannon", ]

m_sha_30y <- gam(rcsBAI_30y ~ Shannon + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_sha_30y <- summary(m_sha_30y)$p.table["Shannon", ]

m_inv_5y <- gam(rcsBAI_5y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_5y <- summary(m_inv_5y)$p.table["InvSimpson", ]

m_inv_10y <- gam(rcsBAI_10y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_10y <- summary(m_inv_10y)$p.table["InvSimpson", ]

m_inv_15y <- gam(rcsBAI_15y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_15y <- summary(m_inv_15y)$p.table["InvSimpson", ]

m_inv_30y <- gam(rcsBAI_30y ~ InvSimpson + Plot, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
c_inv_30y <- summary(m_inv_30y)$p.table["InvSimpson", ]

res <- data.frame(
Index = c("Observed", "Observed", "Observed", "Observed", "Shannon", "Shannon", "Shannon", "Shannon", "InvSimpson", "InvSimpson", "InvSimpson", "InvSimpson"),
Period = c("rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y", "rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y", "rcsBAI_5y", "rcsBAI_10y", "rcsBAI_15y", "rcsBAI_30y"),
Estimate = c(c_obs_5y["Estimate"], c_obs_10y["Estimate"], c_obs_15y["Estimate"], c_obs_30y["Estimate"], c_sha_5y["Estimate"], c_sha_10y["Estimate"], c_sha_15y["Estimate"], c_sha_30y["Estimate"], c_inv_5y["Estimate"], c_inv_10y["Estimate"], c_inv_15y["Estimate"], c_inv_30y["Estimate"]),
SE = c(c_obs_5y["Std. Error"], c_obs_10y["Std. Error"], c_obs_15y["Std. Error"], c_obs_30y["Std. Error"], c_sha_5y["Std. Error"], c_sha_10y["Std. Error"], c_sha_15y["Std. Error"], c_sha_30y["Std. Error"], c_inv_5y["Std. Error"], c_inv_10y["Std. Error"], c_inv_15y["Std. Error"], c_inv_30y["Std. Error"]),
p = c(c_obs_5y["Pr(>|t|)"], c_obs_10y["Pr(>|t|)"], c_obs_15y["Pr(>|t|)"], c_obs_30y["Pr(>|t|)"], c_sha_5y["Pr(>|t|)"], c_sha_10y["Pr(>|t|)"], c_sha_15y["Pr(>|t|)"], c_sha_30y["Pr(>|t|)"], c_inv_5y["Pr(>|t|)"], c_inv_10y["Pr(>|t|)"], c_inv_15y["Pr(>|t|)"], c_inv_30y["Pr(>|t|)"]),
R2_adj= c(summary(m_obs_5y)$r.sq, summary(m_obs_10y)$r.sq, summary(m_obs_15y)$r.sq, summary(m_obs_30y)$r.sq, summary(m_sha_5y)$r.sq, summary(m_sha_10y)$r.sq, summary(m_sha_15y)$r.sq, summary(m_sha_30y)$r.sq, summary(m_inv_5y)$r.sq, summary(m_inv_10y)$r.sq, summary(m_inv_15y)$r.sq, summary(m_inv_30y)$r.sq),
df = c(df.residual(m_obs_5y), df.residual(m_obs_10y), df.residual(m_obs_15y), df.residual(m_obs_30y), df.residual(m_sha_5y), df.residual(m_sha_10y), df.residual(m_sha_15y), df.residual(m_sha_30y), df.residual(m_inv_5y), df.residual(m_inv_10y), df.residual(m_inv_15y), df.residual(m_inv_30y)))

# BH correction for the 12 tests
res$p_adj <- p.adjust(res$p, method = "BH")  

## Writing the table of results
sd_obs <- sd(div_df2_clean$Observed)
sd_sha <- sd(div_df2_clean$Shannon)
sd_inv <- sd(div_df2_clean$InvSimpson)

res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
res$SD_index   <- c(sd_obs, sd_obs, sd_obs, sd_obs, sd_sha, sd_sha, sd_sha, sd_sha, sd_inv, sd_inv, sd_inv, sd_inv)
res$per_SD_pct <- 100 * (exp(res$Estimate * res$SD_index) - 1)
write.csv(res, "TableST5_ECM_OTU_rcsBAI.csv", row.names = FALSE)
res

##         Index     Period     Estimate          SE         p    R2_adj df
## 1    Observed  rcsBAI_5y  0.002945501 0.004484028 0.5127776 0.3043347 99
## 2    Observed rcsBAI_10y  0.002875807 0.004285419 0.5037401 0.2701531 99
## 3    Observed rcsBAI_15y  0.002815023 0.004069542 0.4907262 0.2243564 99
## 4    Observed rcsBAI_30y  0.003546317 0.003693550 0.3393251 0.1651666 99
## 5     Shannon  rcsBAI_5y -0.051521601 0.179899221 0.7751768 0.2890125 99
## 6     Shannon rcsBAI_10y -0.037521047 0.172489049 0.8282453 0.2560762 99
## 7     Shannon rcsBAI_15y -0.021596927 0.164010675 0.8955046 0.2143644 99
## 8     Shannon rcsBAI_30y  0.045344691 0.149004247 0.7615243 0.1575800 99
## 9  InvSimpson  rcsBAI_5y -0.029724208 0.033311630 0.3743921 0.2934624 99
## 10 InvSimpson rcsBAI_10y -0.027782965 0.032026094 0.3877604 0.2614087 99
## 11 InvSimpson rcsBAI_15y -0.025159004 0.030558612 0.4123144 0.2217315 99
## 12 InvSimpson rcsBAI_30y -0.012059048 0.028018727 0.6678457 0.1633334 99
##        p_adj       CI_low    CI_high   SD_index per_SD_pct
## 1  0.8790474 -0.005951784 0.01184279 19.3371096   5.861079
## 2  0.8790474 -0.005627395 0.01137901 19.3371096   5.718508
## 3  0.8790474 -0.005259831 0.01088988 19.3371096   5.594321
## 4  0.8790474 -0.003782489 0.01087512 19.3371096   7.098150
## 5  0.8955046 -0.408480685 0.30543748  0.4904209  -2.495072
## 6  0.8955046 -0.379776742 0.30473465  0.4904209  -1.823284
## 7  0.8955046 -0.347029689 0.30383584  0.4904209  -1.053569
## 8  0.8955046 -0.250312062 0.34100144  0.4904209   2.248709
## 9  0.8790474 -0.095821710 0.03637329  2.5836409  -7.392188
## 10 0.8790474 -0.091329684 0.03576375  2.5836409  -6.926549
## 11 0.8790474 -0.085793921 0.03547591  2.5836409  -6.293425
## 12 0.8955046 -0.067654280 0.04353618  2.5836409  -3.067589

ECM diversity vs tree age

Distributions and transformations

par(mfrow = c(2, 3))
hist(div_df2_clean$Age, main = "Age")
hist(log(div_df2_clean$Age), main = "log(Age)")
hist(div_df2_clean$Observed, main = "OTU richness")   # counts, right-skewed
hist(div_df2_clean$Shannon, main = "Shannon")         # close to symmetric
hist(div_df2_clean$InvSimpson, main = "InvSimpson")   # positive, right-skewed
hist(log(div_df2_clean$InvSimpson), main = "log(InvSimpson)")

par(mfrow = c(1, 1))

# tree age is right-skewed (skewness), log(Age) is more symmetric

OTU richness vs age

# counts and overdispersed: negative binomial, or Gaussian on log(richness)
m_nb <- gam(Observed ~ log(Age) + Habitat + Site, family = nb(), method = "REML", data = div_df2_clean)
m_ln <- gam(log(Observed) ~ log(Age) + Habitat + Site, family = gaussian(), method = "REML", data = div_df2_clean)

# the two AICs are on different response scales; the Jacobian term makes them comparable
AIC(m_nb)

## [1] 891.0162

AIC(m_ln) + 2 * sum(log(div_df2_clean$Observed))

## [1] 878.1596

# is the age effect linear? edf close to 1 = yes, a linear term is enough (we have 1.001)
m_s <- gam(Observed ~ s(log(Age), k = 5) + Habitat + Site, family = nb(), method = "REML", data = div_df2_clean)
summary(m_s)

## 
## Family: Negative Binomial(5.471) 
## Link function: log 
## 
## Formula:
## Observed ~ s(log(Age), k = 5) + Habitat + Site
## 
## Parametric coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.429355   0.091751  37.377  < 2e-16 ***
## Habitattreeline  -0.009272   0.095746  -0.097  0.92285    
## SiteBrooks_Range  0.308181   0.105625   2.918  0.00353 ** 
## SiteInterior      0.212946   0.118570   1.796  0.07250 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value
## s(log(Age)) 1.001  1.001  2.638   0.104
## 
## R-sq.(adj) =  0.0464   Deviance explained = 10.7%
## -REML =    448  Scale est. = 1         n = 106

summary(m_ln)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(Observed) ~ log(Age) + Habitat + Site
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.8486270  0.4434779   8.678 7.23e-14 ***
## log(Age)         -0.1177362  0.0922555  -1.276  0.20481    
## Habitattreeline   0.0000311  0.0933058   0.000  0.99973    
## SiteBrooks_Range  0.3110503  0.1024934   3.035  0.00306 ** 
## SiteInterior      0.2038428  0.1151967   1.770  0.07983 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.0709   Deviance explained = 10.6%
## -REML = 70.387  Scale est. = 0.19924   n = 106

gam.check(m_ln)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-8.845829e-08,-8.845829e-08]
## (score 70.38683 & scale 0.1992414).
## Hessian positive definite, eigenvalue range [50.5,50.5].
## Model rank =  5 / 5

Shannon index vs age

m_sh <- gam(Shannon ~ log(Age) + Habitat + Site, family = gaussian(), method = "REML", data = div_df2_clean)
m_sh_alt <- gam(Shannon ~ log(Age) + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
AIC(m_sh, m_sh_alt) # same response scale, directly comparable

##          df      AIC
## m_sh      6 147.5350
## m_sh_alt  6 150.0633

m_sh_s <- gam(Shannon ~ s(log(Age), k = 5) + Habitat + Site, family = gaussian(), method = "REML", data = div_df2_clean)
summary(m_sh_s)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Shannon ~ s(log(Age), k = 5) + Habitat + Site
## 
## Parametric coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.80619    0.09319  19.382  < 2e-16 ***
## Habitattreeline  -0.09119    0.09820  -0.929 0.355293    
## SiteBrooks_Range  0.37633    0.10787   3.489 0.000721 ***
## SiteInterior      0.20534    0.12124   1.694 0.093426 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df     F p-value
## s(log(Age))   1      1 0.504   0.479
## 
## R-sq.(adj) =  0.0824   Deviance explained = 11.7%
## -REML = 76.209  Scale est. = 0.22071   n = 106

summary(m_sh)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Shannon ~ log(Age) + Habitat + Site
## 
## Parametric coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.11271    0.46675   4.526 1.64e-05 ***
## log(Age)         -0.06896    0.09710  -0.710 0.479180    
## Habitattreeline  -0.09119    0.09820  -0.929 0.355293    
## SiteBrooks_Range  0.37633    0.10787   3.489 0.000721 ***
## SiteInterior      0.20534    0.12124   1.694 0.093426 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.0824   Deviance explained = 11.7%
## -REML = 75.554  Scale est. = 0.22071   n = 106

gam.check(m_sh)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-7.294381e-08,-7.294381e-08]
## (score 75.55356 & scale 0.2207054).
## Hessian positive definite, eigenvalue range [50.5,50.5].
## Model rank =  5 / 5

Inverted Simpson index vs age

m_is <- gam(InvSimpson ~ log(Age) + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
m_is_ln <- gam(log(InvSimpson) ~ log(Age) + Habitat + Site, family = gaussian(), method = "REML", data = div_df2_clean)

AIC(m_is)

## [1] 471.5441

AIC(m_is_ln) + 2 * sum(log(div_df2_clean$InvSimpson)) # Jacobian term

## [1] 469.3961

m_is_s <- gam(InvSimpson ~ s(log(Age), k = 5) + Habitat + Site, family = Gamma(link = "log"), method = "REML", data = div_df2_clean)
summary(m_is_s)

## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## InvSimpson ~ s(log(Age), k = 5) + Habitat + Site
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.465633   0.090302  16.230  < 2e-16 ***
## Habitattreeline  -0.003046   0.095209  -0.032 0.974541    
## SiteBrooks_Range  0.365310   0.104801   3.486 0.000728 ***
## SiteInterior      0.190485   0.118247   1.611 0.110323    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df     F p-value
## s(log(Age)) 1.106  1.205 0.015   0.973
## 
## R-sq.(adj) =  0.0629   Deviance explained = 11.1%
## -REML = 238.44  Scale est. = 0.20698   n = 106

summary(m_is)

## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## InvSimpson ~ log(Age) + Habitat + Site
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.452857   0.452354   3.212 0.001770 ** 
## log(Age)          0.002665   0.094102   0.028 0.977463    
## Habitattreeline  -0.001285   0.095173  -0.014 0.989253    
## SiteBrooks_Range  0.363344   0.104545   3.475 0.000753 ***
## SiteInterior      0.194176   0.117502   1.653 0.101532    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.0616   Deviance explained = 10.9%
## -REML = 237.78  Scale est. = 0.2073    n = 106

gam.check(m_is)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-2.488994e-08,-2.488994e-08]
## (score 237.782 & scale 0.2072966).
## Hessian positive definite, eigenvalue range [53.92398,53.92398].
## Model rank =  5 / 5

ECM relative abundance vs age

# ECM relative abundance: share of ECM reads per sample
ecm_reads <- colSums(filt_count_tab_no_r[ECM_OTUs_present, ])
all_reads <- colSums(filt_count_tab_no_r)
ecm_abund <- data.frame(Sample = names(ecm_reads), ECM_prop = ecm_reads / all_reads)
div_df2_clean <- merge(div_df2_clean, ecm_abund, by = "Sample")
summary(div_df2_clean$ECM_prop)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.08556 0.59175 0.72107 0.68575 0.80287 0.94489

par(mfrow = c(1, 2))
hist(div_df2_clean$ECM_prop, main = "ECM share of reads")                     # left-skewed
hist(qlogis(div_df2_clean$ECM_prop), main = "logit(ECM share)")               # more symmetric

par(mfrow = c(1, 1))

# proportions between 0 and 1: beta regression, or Gaussian on the logit
m_ab <- gam(ECM_prop ~ log(Age) + Habitat + Site, family = betar(link = "logit"), method = "REML", data = div_df2_clean)
m_ab_lg <- gam(qlogis(ECM_prop) ~ log(Age) + Habitat + Site, family = gaussian(), method = "REML", data = div_df2_clean)

# comparable AICs: the logit model needs the Jacobian term
AIC(m_ab)

## [1] -103.3698

AIC(m_ab_lg) + 2 * sum(log(div_df2_clean$ECM_prop * (1 - div_df2_clean$ECM_prop)))

## [1] -108.2494

# is the age effect linear? edf close to 1 = yes
m_ab_s <- gam(ECM_prop ~ s(log(Age), k = 5) + Habitat + Site, family = betar(link = "logit"), method = "REML", data = div_df2_clean)
summary(m_ab_s)

## 
## Family: Beta regression(8.465) 
## Link function: logit 
## 
## Formula:
## ECM_prop ~ s(log(Age), k = 5) + Habitat + Site
## 
## Parametric coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.3414     0.1327   2.572 0.010102 *  
## Habitattreeline    0.2798     0.1445   1.936 0.052919 .  
## SiteBrooks_Range   0.5863     0.1578   3.716 0.000202 ***
## SiteInterior       0.2926     0.1755   1.667 0.095549 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df Chi.sq p-value
## s(log(Age))   1      1  2.703     0.1
## 
## R-sq.(adj) =  0.166   Deviance explained = 20.9%
## -REML = -51.232  Scale est. = 1         n = 106

summary(m_ab)

## 
## Family: Beta regression(8.465) 
## Link function: logit 
## 
## Formula:
## ECM_prop ~ log(Age) + Habitat + Site
## 
## Parametric coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.4050     0.7000   2.007 0.044725 *  
## log(Age)          -0.2393     0.1456  -1.644 0.100179    
## Habitattreeline    0.2798     0.1445   1.936 0.052919 .  
## SiteBrooks_Range   0.5863     0.1578   3.716 0.000202 ***
## SiteInterior       0.2926     0.1755   1.667 0.095548 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.166   Deviance explained = 20.9%
## -REML = -51.887  Scale est. = 1         n = 106

gam.check(m_ab)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [3.609446e-05,3.609446e-05]
## (score -51.88722 & scale 1).
## Hessian positive definite, eigenvalue range [56.24001,56.24001].
## Model rank =  5 / 5

# age effect at the plot level
nd_ab <- expand.grid(Age = seq(min(div_df2_clean$Age), max(div_df2_clean$Age), length.out = 50),
                     Habitat = unique(div_df2_clean$Habitat),
                     Site = unique(div_df2_clean$Site))
nd_ab$pred <- predict(m_ab, newdata = nd_ab, type = "response")

ggplot(div_df2_clean, aes(x = Age, y = ECM_prop)) +
  geom_point(alpha = 0.6, color = "#EE7733") +
  geom_line(data = nd_ab, aes(y = pred), color = "#009988", linewidth = 1) +
  facet_grid(Habitat ~ Site) +
  scale_x_log10() +
  theme_bw() + labs(x = "Tree age (years, log scale)", y = "ECM share of reads")

Summary of the four models

# estimates are on the scale of each model's link: log for richness, identity for Shannon,
# log for InvSimpson and logit for the ECM share

res <- data.frame(Response = c("Observed (log)", "Shannon", "InvSimpson", "ECM share (logit)"),
                  Estimate = c(summary(m_ln)$p.table["log(Age)", 1],
                               summary(m_sh)$p.table["log(Age)", 1],
                               summary(m_is)$p.table["log(Age)", 1],
                               summary(m_ab)$p.table["log(Age)", 1]),
                  SE = c(summary(m_ln)$p.table["log(Age)", 2],
                         summary(m_sh)$p.table["log(Age)", 2],
                         summary(m_is)$p.table["log(Age)", 2],
                         summary(m_ab)$p.table["log(Age)", 2]),
                  p = c(summary(m_ln)$p.table["log(Age)", 4],
                        summary(m_sh)$p.table["log(Age)", 4],
                        summary(m_is)$p.table["log(Age)", 4],
                        summary(m_ab)$p.table["log(Age)", 4]),
                  R2_adj = c(summary(m_ln)$r.sq, summary(m_sh)$r.sq, summary(m_is)$r.sq, summary(m_ab)$r.sq))

res$CI_low  <- res$Estimate - 1.96 * res$SE
res$CI_high <- res$Estimate + 1.96 * res$SE
res$p_adj   <- p.adjust(res$p, method = "BH")   # four tests, as elsewhere
res <- res[, c("Response", "Estimate", "SE", "CI_low", "CI_high", "R2_adj", "p", "p_adj")]

# age effect for OTU richness, at the plot level
af <- div_df2_clean[div_df2_clean$Site == "Alaska_Range" & div_df2_clean$Habitat == "forest", ]
at <- div_df2_clean[div_df2_clean$Site == "Alaska_Range" & div_df2_clean$Habitat == "treeline", ]
bf <- div_df2_clean[div_df2_clean$Site == "Brooks_Range" & div_df2_clean$Habitat == "forest", ]
bt <- div_df2_clean[div_df2_clean$Site == "Brooks_Range" & div_df2_clean$Habitat == "treeline", ]
iff <- div_df2_clean[div_df2_clean$Site == "Interior" & div_df2_clean$Habitat == "forest", ]
it <- div_df2_clean[div_df2_clean$Site == "Interior" & div_df2_clean$Habitat == "treeline", ]

# prediction range per panel, so the line is not extrapolated beyond the ages sampled there
nd <- rbind(
  data.frame(Age = seq(min(af$Age), max(af$Age), length.out = 50), Habitat = "forest",   Site = "Alaska_Range"),
  data.frame(Age = seq(min(at$Age), max(at$Age), length.out = 50), Habitat = "treeline", Site = "Alaska_Range"),
  data.frame(Age = seq(min(bf$Age), max(bf$Age), length.out = 50), Habitat = "forest",   Site = "Brooks_Range"),
  data.frame(Age = seq(min(bt$Age), max(bt$Age), length.out = 50), Habitat = "treeline", Site = "Brooks_Range"),
  data.frame(Age = seq(min(iff$Age), max(iff$Age), length.out = 50), Habitat = "forest",   Site = "Interior"),
  data.frame(Age = seq(min(it$Age), max(it$Age), length.out = 50), Habitat = "treeline", Site = "Interior"))
nd$pred <- exp(predict(m_ln, newdata = nd)) # log(richness) -> richness

# one prediction grid per plot, so the line is not drawn beyond the ages sampled there
ggplot(div_df2_clean, aes(x = Age, y = Observed)) +
  geom_point(alpha = 0.6, color = "#EE7733") +
  geom_line(data = nd, aes(y = pred), color = "#009988", linewidth = 1) +
  facet_grid(Habitat ~ Site) +
  scale_x_log10() +
  theme_bw() +
  labs(x = "Tree age (years, log scale)", y = "ECM OTU richness")

res

##            Response     Estimate         SE     CI_low    CI_high     R2_adj
## 1    Observed (log) -0.117736190 0.09225550 -0.2985570 0.06308458 0.07088234
## 2           Shannon -0.068964550 0.09709770 -0.2592760 0.12134694 0.08235453
## 3        InvSimpson  0.002664906 0.09410195 -0.1817749 0.18710472 0.06163536
## 4 ECM share (logit) -0.239315890 0.14557057 -0.5246342 0.04600243 0.16647056
##           p     p_adj
## 1 0.2048128 0.4096255
## 2 0.4791801 0.6389068
## 3 0.9774633 0.9774633
## 4 0.1001792 0.4007169

## Writing the table of results
res$df      <- c(df.residual(m_ln), df.residual(m_sh), df.residual(m_is), df.residual(m_ab))
res$CI_low  <- res$Estimate - qt(0.975, res$df) * res$SE
res$CI_high <- res$Estimate + qt(0.975, res$df) * res$SE
write.csv(res, "Table-ST6_div_vs_age.csv", row.names = FALSE)

GAM: Tree growth vs RAF diversity

Kristina Kuprina

2026-09-30

Read data

NTC treatment

All OTUs

Data preparation

Model check: Gaussian vs Gamma (Shannon, BAI 5 years)

RAF diversity vs raw BAI

RAF diversity vs detrended BAI

RAF diversity vs RCS BAI

RAF diversity vs DBH and age

ECM OTUs

Data preparation (ECM)

RAF diversity vs raw BAI (ECM)

RAF diversity vs detrended BAI (ECM)

RAF diversity vs rcsBAI (ECM)

ECM diversity vs tree age

Distributions and transformations

OTU richness vs age

Shannon index vs age

Inverted Simpson index vs age

ECM relative abundance vs age

Summary of the four models