Tomicus piniperda, wildfire, phloem nutrient chemistry and fungal communities
Kerri Kluting
2021-12-02
- The data
- Analysis Part I: fungal communities
- Analysis Part II: phloem chemistry
- Quick look at the phloem chemistry data
- Burnt vs unburnt C:N, C:P and N:P pre-colonization
- Fig. 6 - nutrient ratio values for all samples
- Fig. S2 - nutrient ratio values for all samples organized by tree within burnt vs unburnt sites
- Linear mixed-effects models, partially Bayesian implementation
- Fig. S3 - model-predicted values
- Comparison of model-predicted vs mean nutrient ratio values (subsets of Fig. 5 and Fig. S3)
- Checking model assumptions
- \(R^{2}\), ICC, Wald \(\chi^{2}\) test, 95% CI, and summary for each model
- Linear mixed-effects models, partially Bayesian implementation, burnt and unburnt subsets
- Session Info
The data
Sampling sites (Fig. 1)
Phloem samples were taken from 12 trees, three from each of four forest sites in southern Sweden: Hälleskogsbrännan (HSB; 59°54’35.5“N 16°08’42.2”E), Ecopark Öjesjön (OJ; 59°51’04.2“N 16°14’57.8”E), Knutby (KB; 59°56’46.1“N 18°18’00.4”E and Skyttorp (SK; 60°06’32.4“N 17°49’34.0”E). In addition, 48 adult male and 48 female beetles were sampled from each site.
HSB and OJ were hit by wildfire during summer 2014, while KB and SK were not.
Sampling scheme
Beetles were paired together (due to such small amount of tissue), resulting in 24 female and 24 male beetle samples from each site.
Five trees from each site were felled in the spring, prior to potential bark beetle infestation, and a single phloem sample was taken from each tree (‘t0’/pre-colonization phloem). Three trees from each site were then chosen for inclusion in this study based on successful beetle colonization.
Phloem samples were taken from each tree at two time points: first when the bark beetle larvae where in the first instar stage of development (‘t1’/ 1st instar phloem), and again when they were in the 3rd instar of larval development (‘t2’/3rd instar phloem). At each of these time points, ten phloem samples were taken from each tree: two from uncolonized phloem tissue, and eight from the edges of bark beetle galleries.
In total, 192 beetle samples and 252 phloem samples were taken for inclusion in this study. DNA was sampled from a body wash and from the guts from each of the beetle samples and from all phloem samples, and sequenced using PacBio.
Carbon (C), nitrogen (N), and phosphorus (P) concentrations were measured for each phloem sample, and molar C:N, C:P, and N:P ratios were calculated.
Data files
Statistical analyses conducted in publication can be reproduced using the following files:
- metadata file: ‘tp_metadata_20211202.txt’
- OTU / sample matrix: ‘tp_initial_otu_table.txt’
- SINTAX (Edgar, 2016) taxonomy predictions for each OTU sequence from which the ITS2 region could be successfully extracted: ‘tp_sintax_results.txt’
Additional data files associated with this publication:
OTU representative sequences for the OTUs included in the initial OTU table (‘tp_initial_otu_table.txt’): ‘tp_otu_seqs.fasta’
output files from OTU clustering (‘scata3997_JUVgMBeZyA.zip’) generated via the SCATA pipeline (https://scata.mykopat.slu.se/), which contains:
- unfiltered OTU representative sequences (‘all_clusters_scata3997.fas’)
- unfiltered OTU table (‘all_tag_by_cluster_counts.txt’)tagged primer sequences (sample primer tag combinations can be found in the metadata file)
raw sequence data
Analysis Part I: fungal communities
Overview of community sequencing dataset
Initial OTU table
# otutab_in is file 'tp_initial_otu_table.txt'
initial_otutab <- otutab_in
# get list of fungi only OTUs
sintax_results_fungi <-
sintax_results[startsWith(sintax_results$cutoff_0.8, "d:Fungi"),]
# subset the initial_otutab to include only the fungal OTUs
fungi_otus <- row.names(sintax_results_fungi)
# filter initial_otutab
otutab <- initial_otutab[, names(initial_otutab) %in% fungi_otus]There are 412 samples, 362 OTUs, and 15854 sequence reads in OTU table (SCATA output, samples with primer tag combinations not included in this study plus singleton OTUs filtered out).
fungal OTUs only
There are 412 samples, 290 OTUs, and 13397 sequence reads in the OTU table after filtering non-fungal OTUs out.
Filtered OTU table
Remove positive control sequence
One OTU (scata3997_1; 1084 reads) that corresponds to positive control DNA used in PCR, so it should be removed (DNA extracted from Agaricus bisporus was used as a positive control, and scata3997_1 matches A. bisporus)
# remove OTU_1 from otutab
otutab <- otutab[, names(otutab) != "scata3997_1"]
# remove samples with 0 reads left
otutab <- otutab[rowSums(otutab) > 0,] There are 401 samples, 289 OTUs, and 12313 sequence reads in the fungal community dataset after removing scata3997_1.
Remove samples with less than 10 reads
Some of the samples have sequencing depths of less than 10 reads:
Total samples with less than 10 reads: 70
Minimum sequencing depth: 1
Maximum sequencing depth: 156 (sample HBTPFB11)
# remove samples with less than 10 reads for downstream analysis
otutab_min10 <- subset(otutab, rowSums(otutab) > 9)
# remove OTUs that have zero reads after removing samples with less than 10 reads
otutab_min10 <- otutab_min10[, colSums(otutab_min10) > 0]After removing samples with less than 10 reads, the community dataset used for downstream analysis includes:
Total reads: 11962
Total samples: 331
Total OTUs: 286
Fig. 2 nMDS beetle samples
nMDS results
Call:
metaMDS(comm = otutab_beetle_ra, distance = "bray", k = 2, trymax = 200)
global Multidimensional Scaling using monoMDS
Data: otutab_beetle_ra
Distance: bray
Dimensions: 2
Stress: 0.2405601
Stress type 1, weak ties
No convergent solutions - best solution after 200 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'otutab_beetle_ra'nMDS plot
plot saved as ‘tp_Fig2_nMDS_beetles.pdf’
Figure 2. Non-metric multidimensional scaling (nMDS) ordination showing fungal community differentiation in bark beetle samples associated with fire history of sampling locality. Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 167 samples, 258 OTUs and 7155 sequence reads represented in the plot. Plot stress = 0.241.
This is the same plot as above, just highlighting sex and body region:
PERMANOVA
Fire history (burnt vs unburnt)
The most conspicuous pattern for beetle samples is a difference based on burnt vs unburnt collection locality. This is supported by PERMANOVA. A test for homogeneity of group dispersion tests for a significant difference in dispersion, which can lead to a false positive PERMANOVA test. However, a clear distinction between samples from burnt and unburnt sites can be observed.
# PERMANOVA test for significant difference between beetle sample fungal communities based on fire history (burnt vs unburnt)
adonis2(
otutab_beetle_ra ~ fire_hist,
data = meta_seq_beetle,
by = 'terms',
method = 'bray'
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_beetle_ra ~ fire_hist, data = meta_seq_beetle, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 4.473 0.07355 13.099 0.001 ***
Residual 165 56.346 0.92645
Total 166 60.819 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$fire_hist)) Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 1 0.27326 0.273258 22.387 0.000004764 ***
Residuals 165 2.01402 0.012206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sex (female vs male) and body region (body vs gut)
Grouping of samples based on beetle sex or body region is not apparent in ordination. PERMANOVA indicates that a small amount of variation can be explained by sex and body region, but this is a negligible amount. PERMANOVA indicates that a slightly larger amount of variation can be attributed to beetle body region, however, this may be due to the significantly different group dispersions.
sex
# PERMANOVA test for significant difference between beetle sample fungal communities based on sex (female vs male)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_beetle, site)
adonis2(
otutab_beetle_ra ~ sex,
data = meta_seq_beetle,
by = 'terms',
method = 'bray',
permutations = perm
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_beetle, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_beetle_ra ~ sex, data = meta_seq_beetle, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
sex 1 0.580 0.00954 1.5885 0.025 *
Residual 165 60.239 0.99046
Total 166 60.819 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$sex)) Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 1 0.02047 0.020472 1.6313 0.2033
Residuals 165 2.07062 0.012549body region
# PERMANOVA test for significant difference between beetle sample fungal communities based on body region (body vs gut)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_beetle, site)
adonis2(
otutab_beetle_ra ~ region,
data = meta_seq_beetle,
by = 'terms',
method = 'bray',
permutations = perm
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_beetle, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_beetle_ra ~ region, data = meta_seq_beetle, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
region 1 0.953 0.01567 2.6264 0.005 **
Residual 165 59.866 0.98433
Total 166 60.819 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$region)) Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 1 0.04267 0.042667 3.5917 0.05982 .
Residuals 165 1.96007 0.011879
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fig. S1 displaying body and gut samples from the same individuals
plot saved as ‘tp_FigS1_nMDS_beetles_fig2_with_lines.pdf’
Figure S1. Plot from Fig. 2 with lines added to connect body and gut samples that originate from the same individual beetle sample, demonstrating that there is not a pattern driven by the non-independence of the samples.
IDs for beetle samples from which both body and gut samples are included in analysis:
$beetle_individual
[1] "HBTPF10" "HBTPF13" "HBTPF8" "HBTPM4" "KBTPF1" "KBTPF10" "KBTPF11" "KBTPF13" "KBTPF14" "KBTPF16" "KBTPF2" "KBTPF21" "KBTPF23" "KBTPF5" "KBTPM1"
[16] "KBTPM23" "OJTPF14" "OJTPF16" "OJTPF20" "OJTPF21" "OJTPF22" "OJTPF24" "OJTPF25" "OJTPF26" "OJTPM27" "OJTPM30" "OJTPM31" "OJTPM32" "OJTPM33" "OJTPM35"
[31] "SKTPF1" "SKTPF13" "SKTPF16" "SKTPF18" "SKTPF20" "SKTPF21" "SKTPF22" "SKTPF23" "SKTPF5" "SKTPF7" "SKTPM1" "SKTPM2" "SKTPM3" "SKTPM4" "SKTPM8"
[46] "SKTPM9"Fig. 3 - nMDS phloem samples
nMDS results
Call:
metaMDS(comm = otutab_phloem_ra, distance = "bray", k = 2, trymax = 200)
global Multidimensional Scaling using monoMDS
Data: otutab_phloem_ra
Distance: bray
Dimensions: 2
Stress: 0.1859736
Stress type 1, weak ties
No convergent solutions - best solution after 200 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'otutab_phloem_ra'nMDS plot
plot saved as ‘tp_Fig3_nMDS_phloem.pdf’
Figure 3. Non-metric multidimensional scaling (nMDS) ordination displaying differentiation of fungal communities in phloem samples based on fire history of sampling locality. Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 164 samples, 96 OTUs and 4807 sequence reads represented in the plot. Plot stress = 0.186.
PERMANOVA
Fire history (burnt vs unburnt)
# PERMANOVA test for significant difference between phloem sample fungal communities based fire history (burnt vs unburnt)
adonis2(
otutab_phloem_ra ~ fire_hist,
data = meta_seq_phloem,
by = 'terms',
method = 'bray'
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_phloem_ra ~ fire_hist, data = meta_seq_phloem, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 2.033 0.05439 9.3182 0.001 ***
Residual 162 35.346 0.94561
Total 163 37.379 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_phloem_ra, method = "bray"), meta_seq_phloem$fire_hist)) Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 1 0.0003 0.0003371 0.0149 0.9029
Residuals 162 3.6593 0.0225883Uncolonized phloem vs colonized first instar vs colonized 3rd instar
# PERMAVONA test of significance for samples grouped (into 3 groups) as 1) uncolonized samples (pre-col, 1st instar, and 3rd instar collapsed), 2) 1st instar colonized, and 3) 3rd instar colonized)
# 'hm_grp' is grouping in Fig 3 heatmap and nMDS below
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem, site)
adonis2(
otutab_phloem_ra ~ hm_grp,
data = meta_seq_phloem,
by = 'terms',
method = 'bray',
permutations = perm
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_phloem, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_phloem_ra ~ hm_grp, data = meta_seq_phloem, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
hm_grp 5 7.211 0.19293 7.5539 0.005 **
Residual 158 30.167 0.80707
Total 163 37.379 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwise PERMAVONA
Uncolonized phloem vs colonized first instar, and colonized first instar vs colonized 3rd instar
test of significance for uncolonized samples (pre-colonization, 1st instar, and 3rd instar combined) vs 1st instar colonized samples, and for 1st instar vs 3rd instar colonized samples)
# PERMANOVA for uncolonized vs 1st instar colonized samples
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem_uncol_v_1st, site)
permanova_uv1 <- adonis2(
otutab_phloem_uncol_vs_1st_ra ~ hm_grp,
data = meta_seq_phloem_uncol_v_1st,
by = 'terms',
method = 'bray',
permutations = perm
)
permanova_uv1 Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_phloem_uncol_v_1st, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_phloem_uncol_vs_1st_ra ~ hm_grp, data = meta_seq_phloem_uncol_v_1st, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
hm_grp 3 2.8119 0.14808 4.5772 0.005 **
Residual 79 16.1776 0.85192
Total 82 18.9895 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem_1st_vs_3rd, site)
# PERMANOVA for 1st vs 3rd instar colonized samples
permanova_1v3 <- adonis2(
otutab_phloem_1st_vs_3rd_ra ~ hm_grp,
data = meta_seq_phloem_1st_vs_3rd,
by = 'terms',
method = 'bray',
permutations = perm
)
permanova_1v3 Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_phloem_1st_vs_3rd, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_phloem_1st_vs_3rd_ra ~ hm_grp, data = meta_seq_phloem_1st_vs_3rd, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
hm_grp 3 4.963 0.15488 8.6137 0.005 **
Residual 141 27.081 0.84512
Total 144 32.044 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bonferroni correction of p-values for multiple comparisons (pairwise PERMANOVA tests)
ADJUSTED P-VALUES: 0.01 0.01Fig. 4 - nMDS all samples
In order to investigate the OTUs that may be the highest contributors to community structure, we collapsed the samples to into group based on the results above.
beetle samples were collapsed, as visual inspection of nMDS ordination, PERMANOVA and tests of homogeneity of dispersion suggest that sex (male/female) and region (body/gut) do not explain a large proportion of the observed community variation. We also chose to retain both samples (body and gut) from individual beetle samples when available, as this aspect does not appear to contribute to variation in community composition, either. nMDS does, however indicate a difference between samples from burnt and unburnt sites along axis 1, which is supported by PERMANOVA.
uncolonized phloem samples taken at pre-colonization, 1st instar, and 3rd instar sampling time points were collapsed into a single group because they are so few and appear to cluster together in ordination, and the most notable pattern revealed is a general shift along axis 1 as uncolonized – colonized 1st instar – colonized, 3rd instar, and burnt – unburnt along axis 2, a pattern supported by PERMANOVA. Fire history is also supported by PERMANOVA as an explanatory factor of variation in fungal communities.
Therefore, we analyzed the OTUs that were both present in 30% of samples and had an average relative read abundance of 1% or greater in any of the following eight groups:
beetles, burnt
beetles, unburnt
phloem, no larvae, burnt
phloem, no larvae, unburnt
phloem, 1st instar, burnt
phloem, 1st instar, unburnt
phloem, 3rd instar, burnt
phloem, 3rd instar, unburnt
It is important to note that the beetle and phloem samples were run on two separate sequencing plates, so this could contribute to the differences between them. However, nMDS indicates that the 3rd instar samples are most like the beetle samples, followed by 3rd instar, and then uncolonized samples, which makes sense biologically. For this reason, we believe that the signal here is reliable.
nMDS results
Call:
metaMDS(comm = otutab_min10_ra, distance = "bray", k = 2, trymax = 200)
global Multidimensional Scaling using monoMDS
Data: otutab_min10_ra
Distance: bray
Dimensions: 2
Stress: 0.2230584
Stress type 1, weak ties
No convergent solutions - best solution after 200 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'otutab_min10_ra'nMDS plot
plot saved as ‘tp_Fig4_nMDS_all_samples.pdf’
Figure 4. Non-metric multidimensional scaling (nMDS) ordination displaying fungal com-munity shifts associated with forest fire, sampling type, and substrate (phloem vs. beetle samples). Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 331 samples, 286 OTUs, and 11962sequence reads represented in the plot. Plot stress = 0.223.
This is the same plot, just simplified to show 1) only burnt vs unburnt, and 2) only beetle vs. phloem:
PERMANOVA
Phloem vs beetle samples
# PERMAVONA test of significant difference between fungal communities in phloem vs beetle samples (from all four sites)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_min10, site)
adonis2(
otutab_min10_ra ~ seq_plate,
data = meta_seq_min10,
by = 'terms',
method = 'bray',
permutations = perm
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: with(meta_seq_min10, site)
Permutation: free
Number of permutations: 199
adonis2(formula = otutab_min10_ra ~ seq_plate, data = meta_seq_min10, permutations = perm, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
seq_plate 1 16.327 0.14256 54.701 0.005 **
Residual 329 98.198 0.85744
Total 330 114.525 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fire history (burnt vs unburnt)
# PERMAVONA test of significant difference between fungal communities in burnt vs unburnt sites (from all fungal communities, samples from both phloem and beetles)
adonis2(
otutab_min10_ra ~ fire_hist,
data = meta_seq_min10,
by = 'terms',
method = 'bray'
) Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_min10_ra ~ fire_hist, data = meta_seq_min10, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 6.137 0.05359 18.63 0.001 ***
Residual 329 108.388 0.94641
Total 330 114.525 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Phloem vs beetle samples plus fire history (marginal variation)
# marginal PERMANOVA for the two terms (amount of variation explained by each term after other term(s) accounted for)
adonis2(
otutab_min10_ra ~ fire_hist + seq_plate,
data = meta_seq_min10,
by = 'margin',
method = 'bray'
) Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_min10_ra ~ fire_hist + seq_plate, data = meta_seq_min10, method = "bray", by = "margin")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 5.128 0.04478 18.072 0.001 ***
seq_plate 1 15.317 0.13375 53.982 0.001 ***
Residual 328 93.070 0.81266
Total 330 114.525 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pairwise PERMANOVA
Fire history for uncolonized, colonized first, and colonized third separately
pairwise PERMANOVA test of significance between burnt and unburnt phloem samples as grouped in heatmap below
# pairwise PERMANOVA, uncolonized burnt vs uncolonized unburnt
permanova_nolarvae <- adonis2(
otutab_phloem_nolarvae_ra ~ fire_hist,
data = meta_seq_phloem_nolarvae,
by = 'terms',
method = 'bray')
permanova_nolarvae Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_phloem_nolarvae_ra ~ fire_hist, data = meta_seq_phloem_nolarvae, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 0.1392 0.04317 0.767 0.675
Residual 17 3.0863 0.95683
Total 18 3.2255 1.00000# pairwise PERMANOVA, 1st instar burnt vs 1st instar unburnt
meta_seq_phloem_1 <- subset(meta_seq_phloem,
substrate_collapsed_4 == "phloem, 1st instar")
meta_seq_phloem_1 <- droplevels(meta_seq_phloem_1)
otutab_phloem_1 <- otutab_phloem[row.names(meta_seq_phloem_1),]
otutab_phloem_1 <- otutab_phloem_1[, colSums(otutab_phloem_1) > 0]
otutab_phloem_1_ra <- decostand(otutab_phloem_1, method = "total")
permanova_1st <- adonis2(
otutab_phloem_1_ra ~ fire_hist,
data = meta_seq_phloem_1,
by = 'terms',
method = 'bray')
permanova_1st Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_phloem_1_ra ~ fire_hist, data = meta_seq_phloem_1, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 1.4026 0.09677 6.6425 0.001 ***
Residual 62 13.0913 0.90323
Total 63 14.4939 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# pairwise PERMANOVA, 3rd instar burnt vs 3rd instar unburnt
meta_seq_phloem_3 <- subset(meta_seq_phloem,
substrate_collapsed_4 == "phloem, 3rd instar")
meta_seq_phloem_3 <- droplevels(meta_seq_phloem_3)
otutab_phloem_3 <- otutab_phloem[row.names(meta_seq_phloem_3),]
otutab_phloem_3 <- otutab_phloem_3[, colSums(otutab_phloem_3) > 0]
otutab_phloem_3_ra <- decostand(otutab_phloem_3, method = "total")
permanova_3rd <- adonis2(
otutab_phloem_3_ra ~ fire_hist,
data = meta_seq_phloem_3,
by = 'terms',
method = 'bray')
permanova_3rd Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = otutab_phloem_3_ra ~ fire_hist, data = meta_seq_phloem_3, method = "bray", by = "terms")
Df SumOfSqs R2 F Pr(>F)
fire_hist 1 1.8826 0.11861 10.631 0.001 ***
Residual 79 13.9898 0.88139
Total 80 15.8724 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bonferroni correction of p-values for multiple comparisons
ADJUSTED P-VALUES: 1 0.003 0.003Fig. 5 - heatmap of core community OTUs
These are the core OTUs:
[1] "scata3997_0" "scata3997_6" "scata3997_5" "scata3997_7" "scata3997_11" "scata3997_8" "scata3997_28" "scata3997_36" "scata3997_2" "scata3997_3"
[11] "scata3997_13" "scata3997_22" "scata3997_15" "scata3997_17" "scata3997_19"These OTUs account for 72.86 percent of the filtered OTU table (otutab_min10).
heatmap plot
plot saved as ‘tp_Fig5_heatmap.pdf’
Figure 5. Heatmap of the average relative sequence read abundance across samples for 15 OTUs identified as the core community in this study. To the right, OTU number and assigned taxonomy.
Number of samples included in each category:
| site | tree | substrate | n |
|---|---|---|---|
| KB | K2 | 1st instar, uncolonized | 1 |
| KB | K2 | 1st instar, colonized | 5 |
| KB | K2 | 3rd instar, uncolonized | 2 |
| KB | K2 | 3rd instar, colonized | 5 |
| KB | K3 | pre-colonization | 1 |
| KB | K3 | 1st instar, uncolonized | 1 |
| KB | K3 | 1st instar, colonized | 6 |
| KB | K3 | 3rd instar, uncolonized | 2 |
| KB | K3 | 3rd instar, colonized | 6 |
| KB | K5 | pre-colonization | 1 |
| KB | K5 | 1st instar, uncolonized | 1 |
| KB | K5 | 1st instar, colonized | 2 |
| KB | K5 | 3rd instar, colonized | 7 |
| KB | female, body | 21 | |
| KB | female, gut | 12 | |
| KB | male, body | 22 | |
| KB | male, gut | 3 | |
| SK | S2 | 1st instar, colonized | 4 |
| SK | S2 | 3rd instar, colonized | 8 |
| SK | S3 | 1st instar, colonized | 6 |
| SK | S3 | 3rd instar, uncolonized | 1 |
| SK | S3 | 3rd instar, colonized | 7 |
| SK | S5 | 1st instar, colonized | 7 |
| SK | S5 | 3rd instar, colonized | 7 |
| SK | female, body | 21 | |
| SK | female, gut | 11 | |
| SK | male, body | 12 | |
| SK | male, gut | 10 | |
| HSB | H1 | 1st instar, colonized | 8 |
| HSB | H1 | 3rd instar, uncolonized | 1 |
| HSB | H1 | 3rd instar, colonized | 6 |
| HSB | H2 | 1st instar, colonized | 7 |
| HSB | H2 | 3rd instar, colonized | 8 |
| HSB | H3 | 1st instar, uncolonized | 1 |
| HSB | H3 | 1st instar, colonized | 5 |
| HSB | H3 | 3rd instar, colonized | 6 |
| HSB | female, body | 7 | |
| HSB | female, gut | 3 | |
| HSB | male, body | 5 | |
| HSB | male, gut | 2 | |
| OJ | O1 | 1st instar, uncolonized | 1 |
| OJ | O1 | 1st instar, colonized | 1 |
| OJ | O1 | 3rd instar, uncolonized | 2 |
| OJ | O1 | 3rd instar, colonized | 7 |
| OJ | O2 | 1st instar, colonized | 5 |
| OJ | O2 | 3rd instar, uncolonized | 1 |
| OJ | O2 | 3rd instar, colonized | 6 |
| OJ | O4 | 1st instar, uncolonized | 1 |
| OJ | O4 | 1st instar, colonized | 8 |
| OJ | O4 | 3rd instar, uncolonized | 2 |
| OJ | O4 | 3rd instar, colonized | 8 |
| OJ | female, body | 12 | |
| OJ | female, gut | 9 | |
| OJ | male, body | 10 | |
| OJ | male, gut | 7 |
Analysis Part II: phloem chemistry
C, N, and P concentrations were measured for each phloem sample. Three samples—SK5TPWC1, SK3TPWC3, and K3TPWC3—were excluded from analysis because either N or C was measured to be zero. (N and C were measured simultaneously for each sample, so the sample had to be removed entirely from analysis.) For the remaining 249 samples, molar ratios (C:N, C:P, N:P) were calculated.
Quick look at the phloem chemistry data
In general, more tree-to-tree variation is visible than site-to-site in this quick overview of the data.
C:N
C:P
N:P
Burnt vs unburnt C:N, C:P and N:P pre-colonization
Here are the pre-col values again:
One-way ANOVA tests were used to test for difference in nutrient ratios in trees from burnt vs unburnt sites prior to bark beetle colonization, with site treated as a random effects factor.
Based on the results, we fail to reject \(H_0\) that there is no difference in C:P, C:N and N:P ratios of the pre-colonization samples between the burnt and the unburnt sites. These results should be interpreted with caution, however, as our sample size may result in low power of detection.
ANOVA, burnt vs unburnt C:N pre-colonization
# ANOVA, burnt vs unburnt C:N pre-colonization
aov(cn_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0"))
Call:
aov(formula = cn_ratio ~ fire_hist + Error(site), data = subset(meta_chem,
sampling_time == "t0"))
Grand Mean: 95.71526
Stratum 1: site
Terms:
fire_hist Residuals
Sum of Squares 133.5371 215.4130
Deg. of Freedom 1 2
Residual standard error: 10.37817
Estimated effects are balanced
Stratum 2: Within
Terms:
Residuals
Sum of Squares 892.5092
Deg. of Freedom 8
Residual standard error: 10.56237summary(aov(cn_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0")))
Error: site
Df Sum Sq Mean Sq F value Pr(>F)
fire_hist 1 133.5 133.5 1.24 0.381
Residuals 2 215.4 107.7
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 8 892.5 111.6ANOVA, burnt vs unburnt C:P pre-colonization
# ANOVA, burnt vs unburnt C:N pre-colonization
aov(cp_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0"))
Call:
aov(formula = cp_ratio ~ fire_hist + Error(site), data = subset(meta_chem,
sampling_time == "t0"))
Grand Mean: 1993.308
Stratum 1: site
Terms:
fire_hist Residuals
Sum of Squares 1703021.7 945897.1
Deg. of Freedom 1 2
Residual standard error: 687.7125
Estimated effects are balanced
Stratum 2: Within
Terms:
Residuals
Sum of Squares 2262637
Deg. of Freedom 8
Residual standard error: 531.8172summary(aov(cp_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0")))
Error: site
Df Sum Sq Mean Sq F value Pr(>F)
fire_hist 1 1703022 1703022 3.601 0.198
Residuals 2 945897 472949
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 8 2262637 282830ANOVA, burnt vs unburnt N:P pre-colonization
# ANOVA, burnt vs unburnt C:N pre-colonization
aov(np_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0"))
Call:
aov(formula = np_ratio ~ fire_hist + Error(site), data = subset(meta_chem,
sampling_time == "t0"))
Grand Mean: 20.81932
Stratum 1: site
Terms:
fire_hist Residuals
Sum of Squares 135.1451 118.5604
Deg. of Freedom 1 2
Residual standard error: 7.699363
Estimated effects are balanced
Stratum 2: Within
Terms:
Residuals
Sum of Squares 243.4461
Deg. of Freedom 8
Residual standard error: 5.516409summary(aov(np_ratio ~ fire_hist + Error(site),
data = subset(meta_chem, sampling_time == "t0")))
Error: site
Df Sum Sq Mean Sq F value Pr(>F)
fire_hist 1 135.2 135.15 2.28 0.27
Residuals 2 118.6 59.28
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 8 243.4 30.43Fig. 6 - nutrient ratio values for all samples
plot saved as ‘tp_Fig6_chem_by_fire_history_tree_means.pdf’
Figure 6. C:N, C:P and N:P ratios (molar) for all samples. Lines indicate mean ratio value for colonized (solid line) and uncolonized (dashed line) phloem samples for each tree. Faint line connects the pre-colonization value to both the colonized (closed symbols) and uncolonized (open symbols) value for each tree. See Fig. S2 to examine values for samples from each tree separately, and Fig. S3 for model-predicted values plotted.
Fig. S2 - nutrient ratio values for all samples organized by tree within burnt vs unburnt sites
plot saved as ‘tp_FigS2_chem_by_tree.pdf’
Figure S2. C:N, C:P, and N:P ratios for each tree separately. See Fig. 6.
Linear mixed-effects models, partially Bayesian implementation
Here, the relationships between each nutrient ratio (cn_ratio, cp_ratio, and np_ratio) and sampling time (sampling_time; 1st vs 3rd instar), fire history (fire_hist; burnt vs unburnt), and beetle infestation (bb_colonization; colonized vs uncolonized phloem) are examined. The nutrient ratio measured for each tree pre-colonization are treated as a nuisance covariate.
# center (but not scale) pre-col values and code fixed effects factors as sum-to-zero contrasts
meta_chem_t1t2_cent <- meta_chem_t1t2
contrasts(meta_chem_t1t2_cent$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent$bb_colonization) <- contr.sum
contrasts(meta_chem_t1t2_cent$fire_hist) <- contr.sum
meta_chem_t1t2_cent$precol_cn_c <- scale(meta_chem_t1t2_cent$precol_cn, scale = FALSE)
meta_chem_t1t2_cent$precol_cp_c <- scale(meta_chem_t1t2_cent$precol_cp, scale = FALSE)
meta_chem_t1t2_cent$precol_np_c <- scale(meta_chem_t1t2_cent$precol_np, scale = FALSE)# test the random effects structure using intercept-only model
# interaction term [(1|site:tree)]
cn_fit_err1 <- blme::blmer(cn_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)
cp_fit_err1 <- blme::blmer(cp_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)
np_fit_err1 <- blme::blmer(np_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)
# nested term [(1|site/tree)]
# equivalent in this case to [(1|site) + (1|tree)]
cn_fit_err2 <- blme::blmer(cn_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent)
cp_fit_err2 <- blme::blmer(cp_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent)
np_fit_err2 <- blme::blmer(np_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent) # this model fails to converge Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.600704 (tol = 0.002, component 1)# compare models
anova(cn_fit_err1, cn_fit_err2, refit = FALSE) Data: meta_chem_t1t2_cent
Models:
cn_fit_err1: cn_ratio ~ 1 + (1 | site:tree)
cn_fit_err2: cn_ratio ~ 1 + (1 | site/tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cn_fit_err1 3 2024.7 2035.1 -1009.4 2018.7
cn_fit_err2 4 2029.0 2042.8 -1010.5 2021.0 0 1 1anova(cp_fit_err1, cp_fit_err2, refit = FALSE) Data: meta_chem_t1t2_cent
Models:
cp_fit_err1: cp_ratio ~ 1 + (1 | site:tree)
cp_fit_err2: cp_ratio ~ 1 + (1 | site/tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cp_fit_err1 3 3467.1 3477.5 -1730.5 3461.1
cp_fit_err2 4 3468.7 3482.6 -1730.4 3460.7 0.3301 1 0.5656anova(np_fit_err1, np_fit_err2, refit = FALSE) Data: meta_chem_t1t2_cent
Models:
np_fit_err1: np_ratio ~ 1 + (1 | site:tree)
np_fit_err2: np_ratio ~ 1 + (1 | site/tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
np_fit_err1 3 1183.1 1193.5 -588.57 1177.1
np_fit_err2 4 1185.1 1199.0 -588.56 1177.1 0.011 1 0.9166fit the model of interest and test the assumption of equal slopes
In order to treat precol_xx as a nuisance covariate, we must first check that the assumption of equal slopes is met by comparing the models above to models that include an interaction between the pre-colonization value and each of the three explanatory variables (fire_hist, bb_colonization, and sampling_time).
# fit a model with variables of interest plus nuisance covariate pre-col value
cn_fit_nus <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
(1|site:tree), data = meta_chem_t1t2_cent)
cp_fit_nus <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
(1|site:tree), data = meta_chem_t1t2_cent)
np_fit_nus <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
(1|site:tree), data = meta_chem_t1t2_cent)
# fit a model that includes an interaction between the covariate and each of the three fixed effects variables of interest
# see if the model with the interactions is significantly better to test assumption of equal slopes
cn_fit_int <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time +
(1|site:tree), data = meta_chem_t1t2_cent)
cp_fit_int <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time +
(1|site:tree), data = meta_chem_t1t2_cent)
np_fit_int <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time +
(1|site:tree), data = meta_chem_t1t2_cent)
anova(cn_fit_nus, cn_fit_int) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
cn_fit_nus: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + (1 | site:tree)
cn_fit_int: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cn_fit_nus 11 1853.0 1891.1 -915.48 1831.0
cn_fit_int 14 1849.2 1897.8 -910.61 1821.2 9.7412 3 0.0209 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(cp_fit_nus, cp_fit_int) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
cp_fit_nus: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + (1 | site:tree)
cp_fit_int: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cp_fit_nus 11 3378.7 3416.8 -1678.3 3356.7
cp_fit_int 14 3378.5 3427.0 -1675.2 3350.5 6.2103 3 0.1018anova(np_fit_nus, np_fit_int) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
np_fit_nus: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + (1 | site:tree)
np_fit_int: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
np_fit_nus 11 1097.4 1135.6 -537.71 1075.4
np_fit_int 14 1085.5 1134.0 -528.73 1057.5 17.974 3 0.0004454 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model with the interactions is significantly better in the model comparisons for C:N and N:P, indicating that the ANCOVA assumption of equal slopes is violated for these. For the C:P model, both alternatives have nearly the same AIC, but with different number of parameters. Therefore, a mixed-effects model is retained with the interactions included in all three cases. However, it seems that it is mostly driven by the interaction between precol_xx and sampling_time, so models are fitted with only the precol_xx:sampling_time interaction included to see if the model with all three interactions is still significantly better.
# fit a model with only the precol_xx:sampling_time interaction and see if the model with all three interactions is still significantly better
cn_fit_int_red <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)
cp_fit_int_red <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)
np_fit_int_red <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)
anova(cn_fit_int, cn_fit_int_red) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
cn_fit_int_red: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time + (1 | site:tree)
cn_fit_int: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cn_fit_int_red 12 1845.5 1887.1 -910.76 1821.5
cn_fit_int 14 1849.2 1897.8 -910.61 1821.2 0.3075 2 0.8575anova(cp_fit_int, cp_fit_int_red) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
cp_fit_int_red: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time + (1 | site:tree)
cp_fit_int: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
cp_fit_int_red 12 3380.3 3421.9 -1678.2 3356.3
cp_fit_int 14 3378.5 3427.0 -1675.2 3350.5 5.8374 2 0.054 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(np_fit_int, np_fit_int_red) refitting model(s) with ML (instead of REML) Data: meta_chem_t1t2_cent
Models:
np_fit_int_red: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time + (1 | site:tree)
np_fit_int: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time + (1 | site:tree)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
np_fit_int_red 12 1083.9 1125.5 -529.93 1059.9
np_fit_int 14 1085.5 1134.0 -528.73 1057.5 2.4065 2 0.3002# assign the final models (model assumptions tested below)
cn_fit <- cn_fit_int_red
cp_fit <- cp_fit_int_red
np_fit <- np_fit_int_redFig. S3 - model-predicted values
plot saved as ‘tp_FigS3_chem_by_fire_history_predicted_values.pdf’
Figure S3. Predicted C:N, C:P and N:P ratios based on partially Bayesian linear mixed-effects models. Measured C:N, C:P and N:P ratio values for each sample are plotted in Fig. 6 and Fig. S2.
Comparison of model-predicted vs mean nutrient ratio values (subsets of Fig. 5 and Fig. S3)
Checking model assumptions
C:N
1. linearity of predictors
visual inspection, plot of residuals against response values:
2. homogeneity of variance
visual inspection, plot of residuals against fitted values:
statistical test, Levene’s test:
The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).
leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$tree) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 11 3.9057 0.00003454 ***
225
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$bb_colonization) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.7226 0.3962
235leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$sampling_time) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 6.7237 0.01011 *
235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$site) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 2.7622 0.04282 *
233
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. normality of errors
visual inspection, histogram of residuals:
visual inspection, plot of residuals for each explanatory factor:
visual inspection, plot of measured value against residuals:
statistical test, Shapiro-Wilk test:
As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.
shapiro.test(meta_w_preds_and_resids$residuals_cn)
Shapiro-Wilk normality test
data: meta_w_preds_and_resids$residuals_cn
W = 0.98414, p-value = 0.009555shapiro.test(scale(meta_w_preds_and_resids$residuals_cn))
Shapiro-Wilk normality test
data: scale(meta_w_preds_and_resids$residuals_cn)
W = 0.98414, p-value = 0.009555Also, check the normality of random effects:
qqnorm(ranef(cn_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(cn_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(cn_fit)$`site:tree`$`(Intercept)`)
Shapiro-Wilk normality test
data: ranef(cn_fit)$`site:tree`$`(Intercept)`
W = 0.90692, p-value = 0.1948
C:P
1. linearity of predictors
visual inspection, plot of residuals against response values:
2. homogeneity of variance
visual inspection, plot of residuals against fitted values:
statistical test, Levene’s test:
The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).
leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$tree) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 11 1.1438 0.3283
225leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$bb_colonization) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 20.811 0.000008167 ***
235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$sampling_time) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0 0.998
235leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$site) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 4.1986 0.006443 **
233
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. normality of errors
visual inspection, histogram of residuals:
visual inspection, plot of residuals for each explanatory factor:
visual inspection, plot of measured value against residuals:
statistical test, Shapiro-Wilk test:
As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.
shapiro.test(meta_w_preds_and_resids$residuals_cp)
Shapiro-Wilk normality test
data: meta_w_preds_and_resids$residuals_cp
W = 0.96513, p-value = 0.00001518shapiro.test(scale(meta_w_preds_and_resids$residuals_cp))
Shapiro-Wilk normality test
data: scale(meta_w_preds_and_resids$residuals_cp)
W = 0.96513, p-value = 0.00001518Also, check the normality of random effects:
qqnorm(ranef(cp_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(cp_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(cp_fit)$`site:tree`$`(Intercept)`)
Shapiro-Wilk normality test
data: ranef(cp_fit)$`site:tree`$`(Intercept)`
W = 0.97804, p-value = 0.9746
N:P
1. linearity of predictors
visual inspection, plot of residuals against response values:
2. homogeneity of variance
visual inspection, plot of residuals against fitted values:
statistical test, Levene’s test:
The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).
leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$tree) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 11 1.4346 0.1585
225leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$bb_colonization) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 4.0556 0.04516 *
235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$sampling_time) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.132 0.7167
235leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$site) Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 4.6263 0.003654 **
233
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. normality of errors
visual inspection, histogram of residuals:
visual inspection, plot of residuals for each explanatory factor:
visual inspection, plot of measured value against residuals:
statistical test, Shapiro-Wilk test:
As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.
shapiro.test(meta_w_preds_and_resids$residuals_np)
Shapiro-Wilk normality test
data: meta_w_preds_and_resids$residuals_np
W = 0.98206, p-value = 0.004281shapiro.test(scale(meta_w_preds_and_resids$residuals_np))
Shapiro-Wilk normality test
data: scale(meta_w_preds_and_resids$residuals_np)
W = 0.98206, p-value = 0.004281Also, check the normality of random effects:
qqnorm(ranef(np_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(np_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(np_fit)$`site:tree`$`(Intercept)`)
Shapiro-Wilk normality test
data: ranef(np_fit)$`site:tree`$`(Intercept)`
W = 0.95388, p-value = 0.6942
\(R^{2}\), ICC, Wald \(\chi^{2}\) test, 95% CI, and summary for each model
C:N
# R^2
performance::r2(cn_fit) # R2 for Mixed Models
Conditional R2: 0.731
Marginal R2: 0.469# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit) # Intraclass Correlation Coefficient
Adjusted ICC: 0.494
Conditional ICC: 0.262# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cn_ratio
Chisq Df Pr(>Chisq)
fire_hist 0.1576 1 0.691353
bb_colonization 51.5438 1 0.0000000000007002 ***
sampling_time 160.1106 1 < 0.00000000000000022 ***
precol_cn_c 5.2685 1 0.021714 *
fire_hist:bb_colonization 21.3155 1 0.0000038956150564 ***
fire_hist:sampling_time 44.9590 1 0.0000000000201202 ***
bb_colonization:sampling_time 22.4943 1 0.0000021076720233 ***
sampling_time:precol_cn_c 9.4124 1 0.002155 **
fire_hist:bb_colonization:sampling_time 0.7451 1 0.388047
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cn_fit) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 5.5040364 13.32601345
.sigma 9.6762915 11.64216078
(Intercept) 103.3788808 114.05217704
fire_hist1 -6.8616237 4.41613764
bb_colonization1 -7.7790426 -4.31815700
sampling_time1 -7.9413362 -4.47807442
precol_cn_c 0.2093577 1.30118611
fire_hist1:bb_colonization1 2.3967124 5.85759801
fire_hist1:sampling_time1 -6.2650854 -2.67321686
bb_colonization1:sampling_time1 -5.9932267 -2.52999415
sampling_time1:precol_cn_c -0.3662111 -0.08199132
fire_hist1:bb_colonization1:sampling_time1 -2.5007257 0.96250693# model summary
summary(cn_fit) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : 0.0377
Linear mixed model fit by REML ['blmerMod']
Formula: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent
REML criterion at convergence: 1809.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.3942 -0.5249 0.0247 0.5969 2.8124
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 111.9 10.58
Residual 114.8 10.71
Number of obs: 237, groups: site:tree, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 108.71630 3.18115 34.175
fire_hist1 -1.22186 3.36210 -0.363
bb_colonization1 -6.05156 0.88995 -6.800
sampling_time1 -6.21011 0.89058 -6.973
precol_cn_c 0.75530 0.32614 2.316
fire_hist1:bb_colonization1 4.12419 0.88995 4.634
fire_hist1:sampling_time1 -4.46995 0.92365 -4.839
bb_colonization1:sampling_time1 -4.26122 0.89058 -4.785
sampling_time1:precol_cn_c -0.22422 0.07308 -3.068
fire_hist1:bb_colonization1:sampling_time1 -0.76871 0.89058 -0.863
Correlation of Fixed Effects:
(Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
fire_hist1 0.005
bb_colnztn1 -0.174 -0.015
samplng_tm1 -0.006 -0.006 0.021
precol_cn_c 0.004 0.324 0.000 0.000
fr_hst1:b_1 -0.016 -0.165 0.056 0.021 0.000
fr_hst1:s_1 -0.006 -0.006 0.020 0.057 -0.001 0.020
bb_clnz1:_1 0.006 0.006 -0.021 -0.623 0.000 -0.021 -0.056
smplng_1:__ 0.000 -0.002 -0.002 0.004 -0.007 -0.002 0.265 -0.001
fr_h1:_1:_1 0.006 0.006 -0.021 -0.058 0.000 -0.021 -0.601 0.058 -0.001C:P
# R^2
performance::r2(cp_fit) # R2 for Mixed Models
Conditional R2: 0.774
Marginal R2: 0.352# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit) # Intraclass Correlation Coefficient
Adjusted ICC: 0.651
Conditional ICC: 0.422# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cp_ratio
Chisq Df Pr(>Chisq)
fire_hist 0.3713 1 0.5422934
bb_colonization 5.6866 1 0.0170945 *
sampling_time 78.4676 1 < 0.00000000000000022 ***
precol_cp_c 2.5643 1 0.1093029
fire_hist:bb_colonization 0.0402 1 0.8410526
fire_hist:sampling_time 3.2554 1 0.0711892 .
bb_colonization:sampling_time 33.0173 1 0.000000009134 ***
sampling_time:precol_cp_c 0.3703 1 0.5428501
fire_hist:bb_colonization:sampling_time 14.7234 1 0.0001245 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cp_fit) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 197.14654597 461.47439522
.sigma 242.66215053 291.96579135
(Intercept) 1704.17086516 2067.98868988
fire_hist1 -303.35346749 145.34857058
bb_colonization1 13.59968590 100.43771280
sampling_time1 -113.89826848 -27.03452531
precol_cp_c -0.01600956 0.68101675
fire_hist1:bb_colonization1 -46.18613658 40.65189032
fire_hist1:sampling_time1 -35.44100908 64.82133279
bb_colonization1:sampling_time1 -176.63560881 -89.77510786
sampling_time1:precol_cp_c -0.08646784 0.04552046
fire_hist1:bb_colonization1:sampling_time1 -129.24737792 -42.38687697# model summary
summary(cp_fit) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -0.9366
Linear mixed model fit by REML ['blmerMod']
Formula: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent
REML criterion at convergence: 3294.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.7138 -0.6336 -0.0479 0.5273 4.2803
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 134799 367.1
Residual 72198 268.7
Number of obs: 237, groups: site:tree, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 1886.05435 108.31653 17.412
fire_hist1 -79.03938 133.64113 -0.591
bb_colonization1 57.12225 22.33091 2.558
sampling_time1 -70.57560 22.33759 -3.159
precol_cp_c 0.33247 0.20776 1.600
fire_hist1:bb_colonization1 -2.66357 22.33091 -0.119
fire_hist1:sampling_time1 14.51387 25.78271 0.563
bb_colonization1:sampling_time1 -133.09680 22.33676 -5.959
sampling_time1:precol_cp_c -0.02065 0.03394 -0.609
fire_hist1:bb_colonization1:sampling_time1 -85.70857 22.33676 -3.837
Correlation of Fixed Effects:
(Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
fire_hist1 0.006
bb_colnztn1 -0.128 -0.010
samplng_tm1 -0.004 0.000 0.021
precol_cp_c 0.007 0.586 -0.001 0.005
fr_hst1:b_1 -0.012 -0.104 0.057 0.021 -0.001
fr_hst1:s_1 -0.001 0.003 0.003 0.055 0.005 0.003
bb_clnz1:_1 0.004 0.001 -0.021 -0.623 -0.005 -0.021 -0.052
smplng_1:__ 0.006 0.006 -0.031 0.010 0.002 -0.031 0.499 -0.004
fr_h1:_1:_1 0.004 0.001 -0.021 -0.058 -0.005 -0.021 -0.542 0.058 -0.004N:P
# R^2
performance::r2(np_fit) # R2 for Mixed Models
Conditional R2: 0.777
Marginal R2: 0.320# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit) # Intraclass Correlation Coefficient
Adjusted ICC: 0.672
Conditional ICC: 0.457# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: np_ratio
Chisq Df Pr(>Chisq)
fire_hist 1.6250 1 0.2023915
bb_colonization 72.6535 1 < 0.00000000000000022 ***
sampling_time 0.2238 1 0.6361339
precol_np_c 0.8075 1 0.3688732
fire_hist:bb_colonization 23.4658 1 0.000001272 ***
fire_hist:sampling_time 14.2598 1 0.0001592 ***
bb_colonization:sampling_time 7.5743 1 0.0059206 **
sampling_time:precol_np_c 15.6855 1 0.000074797 ***
fire_hist:bb_colonization:sampling_time 18.7486 1 0.000014913 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(np_fit) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 1.62591532 3.7927453
.sigma 1.90514682 2.2922253
(Intercept) 16.03280506 19.0184422
fire_hist1 -2.52923167 0.9616255
bb_colonization1 1.12649924 1.8082652
sampling_time1 0.07541901 0.7573734
precol_np_c -0.12398286 0.4149505
fire_hist1:bb_colonization1 -1.17498193 -0.4932160
fire_hist1:sampling_time1 0.70462091 1.4620188
bb_colonization1:sampling_time1 -0.86713001 -0.1851935
sampling_time1:precol_np_c 0.05063376 0.1480716
fire_hist1:bb_colonization1:sampling_time1 -1.10066365 -0.4187272# model summary
summary(np_fit) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -1.0744
Linear mixed model fit by REML ['blmerMod']
Formula: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent
REML criterion at convergence: 1076.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.9613 -0.6221 -0.0219 0.5938 3.2688
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 9.109 3.018
Residual 4.450 2.110
Number of obs: 237, groups: site:tree, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 17.52536 0.88872 19.720
fire_hist1 -0.78413 1.03940 -0.754
bb_colonization1 1.46847 0.17532 8.376
sampling_time1 0.41604 0.17537 2.372
precol_np_c 0.14547 0.16059 0.906
fire_hist1:bb_colonization1 -0.83301 0.17532 -4.751
fire_hist1:sampling_time1 1.08258 0.19477 5.558
bb_colonization1:sampling_time1 -0.52581 0.17537 -2.998
sampling_time1:precol_np_c 0.09924 0.02506 3.960
fire_hist1:bb_colonization1:sampling_time1 -0.75934 0.17537 -4.330
Correlation of Fixed Effects:
(Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
fire_hist1 0.005
bb_colnztn1 -0.123 -0.010
samplng_tm1 -0.004 -0.001 0.021
precol_np_c 0.006 0.519 -0.001 0.005
fr_hst1:b_1 -0.011 -0.105 0.057 0.021 -0.001
fr_hst1:s_1 -0.001 0.002 0.006 0.056 0.005 0.006
bb_clnz1:_1 0.004 0.001 -0.021 -0.623 -0.005 -0.021 -0.054
smplng_1:__ 0.006 0.006 -0.031 0.008 0.002 -0.031 0.435 -0.004
fr_h1:_1:_1 0.004 0.001 -0.021 -0.058 -0.005 -0.021 -0.563 0.058 -0.004Linear mixed-effects models, partially Bayesian implementation, burnt and unburnt subsets
Fire history does not have a significant main effect in any of the three models, but there are significant interactions that include fire history, in different combinations for each model. To begin teasing apart the relationships between fire history and the remaining two explanatory variables, bark beetle colonization status and sampling time, we can analyze the data for samples from burnt and unburnt sites separately. Fitted models follow the same formula as the models above, only with fire_hist terms excluded.
fit models
# subset by fire history, center precol_XX variables, and apply sum-to-zero contrasts
meta_chem_t1t2_cent_burnt <- subset(meta_chem_t1t2, fire_hist == "burnt") %>% droplevels()
contrasts(meta_chem_t1t2_cent_burnt$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent_burnt$bb_colonization) <- contr.sum
meta_chem_t1t2_cent_burnt$precol_cn_c <- scale(meta_chem_t1t2_cent_burnt$precol_cn, scale = FALSE)
meta_chem_t1t2_cent_burnt$precol_cp_c <- scale(meta_chem_t1t2_cent_burnt$precol_cp, scale = FALSE)
meta_chem_t1t2_cent_burnt$precol_np_c <- scale(meta_chem_t1t2_cent_burnt$precol_np, scale = FALSE)
meta_chem_t1t2_cent_unburnt <- subset(meta_chem_t1t2, fire_hist == "unburnt") %>% droplevels()
contrasts(meta_chem_t1t2_cent_unburnt$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent_unburnt$bb_colonization) <- contr.sum
meta_chem_t1t2_cent_unburnt$precol_cn_c <- scale(meta_chem_t1t2_cent_unburnt$precol_cn, scale = FALSE)
meta_chem_t1t2_cent_unburnt$precol_cp_c <- scale(meta_chem_t1t2_cent_unburnt$precol_cp, scale = FALSE)
meta_chem_t1t2_cent_unburnt$precol_np_c <- scale(meta_chem_t1t2_cent_unburnt$precol_np, scale = FALSE)# fit models following structure from full model (model including samples from both burnt and unburnt sites)
cn_fit_burnt <- blme::blmer(cn_ratio ~ bb_colonization * sampling_time + precol_cn_c +
precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)
cp_fit_burnt <- blme::blmer(cp_ratio ~ bb_colonization * sampling_time + precol_cp_c +
precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)
np_fit_burnt <- blme::blmer(np_ratio ~ bb_colonization * sampling_time + precol_np_c +
precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)
cn_fit_unburnt <- blme::blmer(cn_ratio ~ bb_colonization * sampling_time + precol_cn_c +
precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)
cp_fit_unburnt <- blme::blmer(cp_ratio ~ bb_colonization * sampling_time + precol_cp_c +
precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)
np_fit_unburnt <- blme::blmer(np_ratio ~ bb_colonization * sampling_time + precol_np_c +
precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)Model results
C:N burnt
# R^2
performance::r2(cn_fit_burnt) # R2 for Mixed Models
Conditional R2: 0.762
Marginal R2: 0.209# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit_burnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.699
Conditional ICC: 0.553# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit_burnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cn_ratio
Chisq Df Pr(>Chisq)
bb_colonization 59.7424 1 0.00000000000001081 ***
sampling_time 18.8727 1 0.00001397384680382 ***
precol_cn_c 0.3335 1 0.563587
bb_colonization:sampling_time 7.0376 1 0.007982 **
sampling_time:precol_cn_c 2.7833 1 0.095251 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cn_fit_burnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 6.59502112 22.69988827
.sigma 10.07671203 13.06995771
(Intercept) 101.56126974 123.29901915
bb_colonization1 -12.74964396 -7.60186605
sampling_time1 -5.05382053 0.09395739
precol_cn_c -0.88052444 1.99789667
bb_colonization1:sampling_time1 -6.06639004 -0.91861212
sampling_time1:precol_cn_c -0.04028855 0.50908560# model summary
summary(cn_fit_burnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -1.2626
Linear mixed model fit by REML ['blmerMod']
Formula: cn_ratio ~ bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_burnt
REML criterion at convergence: 932.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.64507 -0.55870 -0.08217 0.63479 2.49597
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 308.9 17.57
Residual 133.1 11.54
Number of obs: 120, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 112.4301 7.2946 15.413
bb_colonization1 -10.1758 1.3165 -7.729
sampling_time1 -2.4799 1.3165 -1.884
precol_cn_c 0.5587 0.9674 0.578
bb_colonization1:sampling_time1 -3.4925 1.3165 -2.653
sampling_time1:precol_cn_c 0.2344 0.1405 1.668
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.108
samplng_tm1 0.000 0.000
precol_cn_c 0.000 0.000 0.000
bb_clnz1:_1 0.000 0.000 -0.600 0.000
smplng_1:__ 0.000 0.000 0.000 0.000 0.000C:N unburnt
# R^2
performance::r2(cn_fit_unburnt) # R2 for Mixed Models
Conditional R2: 0.814
Marginal R2: 0.742# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit_unburnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.277
Conditional ICC: 0.072# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit_unburnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cn_ratio
Chisq Df Pr(>Chisq)
bb_colonization 4.003 1 0.04542 *
sampling_time 265.206 1 < 0.00000000000000022 ***
precol_cn_c 16.026 1 0.00006246412 ***
bb_colonization:sampling_time 22.664 1 0.00000192930 ***
sampling_time:precol_cn_c 36.421 1 0.00000000159 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cn_fit_unburnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.9878653 6.8873900
.sigma 7.6181318 9.9154532
(Intercept) 101.3352645 108.4777333
bb_colonization1 -3.9318497 0.1943767
sampling_time1 -11.9814505 -7.8535201
precol_cn_c 0.5513502 1.1491094
bb_colonization1:sampling_time1 -7.0953268 -2.9673907
sampling_time1:precol_cn_c -0.5742674 -0.2929042# model summary
summary(cn_fit_unburnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : 1.4356
Linear mixed model fit by REML ['blmerMod']
Formula: cn_ratio ~ bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_unburnt
REML criterion at convergence: 841.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.1720 -0.5802 -0.0368 0.5781 2.1912
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 29.29 5.412
Residual 76.26 8.733
Number of obs: 117, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 104.91608 2.44795 42.859
bb_colonization1 -1.90597 1.05431 -1.808
sampling_time1 -9.92338 1.05561 -9.401
precol_cn_c 0.85035 0.20777 4.093
bb_colonization1:sampling_time1 -5.02546 1.05562 -4.761
sampling_time1:precol_cn_c -0.43399 0.07191 -6.035
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.277
samplng_tm1 -0.017 0.040
precol_cn_c -0.001 -0.001 0.001
bb_clnz1:_1 0.017 -0.040 -0.644 -0.001
smplng_1:__ 0.001 -0.003 0.001 -0.015 -0.002C:P burnt
# R^2
performance::r2(cp_fit_burnt) # R2 for Mixed Models
Conditional R2: 0.724
Marginal R2: 0.310# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit_burnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.601
Conditional ICC: 0.415# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit_burnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cp_ratio
Chisq Df Pr(>Chisq)
bb_colonization 3.0299 1 0.08174 .
sampling_time 19.4656 1 0.00001024 ***
precol_cp_c 3.5148 1 0.06082 .
bb_colonization:sampling_time 1.9036 1 0.16768
sampling_time:precol_cp_c 0.1096 1 0.74056
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cp_fit_burnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 133.09046863 476.3141237
.sigma 262.91147842 341.0082471
(Intercept) 1859.40513592 2318.8910835
bb_colonization1 -7.36951045 126.9411554
sampling_time1 -159.95116186 -25.6404960
precol_cp_c 0.11671905 0.9987833
bb_colonization1:sampling_time1 -114.54356382 19.7671020
sampling_time1:precol_cp_c -0.08664678 0.1219758# model summary
summary(cp_fit_burnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -0.613
Linear mixed model fit by REML ['blmerMod']
Formula: cp_ratio ~ bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_burnt
REML criterion at convergence: 1691.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.3517 -0.6346 -0.0307 0.5378 3.7755
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 136338 369.2
Residual 90599 301.0
Number of obs: 120, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 2089.14801 154.60499 13.513
bb_colonization1 59.78582 34.34644 1.741
sampling_time1 -92.79583 34.34644 -2.702
precol_cp_c 0.55775 0.29750 1.875
bb_colonization1:sampling_time1 -47.38823 34.34644 -1.380
sampling_time1:precol_cp_c 0.01766 0.05335 0.331
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.133
samplng_tm1 0.000 0.000
precol_cp_c 0.000 0.000 0.000
bb_clnz1:_1 0.000 0.000 -0.600 0.000
smplng_1:__ 0.000 0.000 0.000 0.000 0.000C:P unburnt
# R^2
performance::r2(cp_fit_unburnt) # R2 for Mixed Models
Conditional R2: 0.843
Marginal R2: 0.209# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit_unburnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.802
Conditional ICC: 0.634# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit_unburnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cp_ratio
Chisq Df Pr(>Chisq)
bb_colonization 2.9804 1 0.08428 .
sampling_time 80.2984 1 < 0.00000000000000022 ***
precol_cp_c 0.0930 1 0.76045
bb_colonization:sampling_time 62.6113 1 0.000000000000002518 ***
sampling_time:precol_cp_c 2.1756 1 0.14021
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(cp_fit_unburnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 177.2625408 593.65591054
.sigma 198.8528551 258.82711679
(Intercept) 1395.5402688 1961.49638884
bb_colonization1 2.1180279 110.01652103
sampling_time1 -102.4289405 5.41979579
precol_cp_c -0.4326445 0.65172278
bb_colonization1:sampling_time1 -272.2899292 -164.43880385
sampling_time1:precol_cp_c -0.1405676 0.02004044# model summary
summary(cp_fit_unburnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -2.0988
Linear mixed model fit by REML ['blmerMod']
Formula: cp_ratio ~ bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_unburnt
REML criterion at convergence: 1590.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.83998 -0.65983 -0.09913 0.62621 2.39668
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 210862 459.2
Residual 52040 228.1
Number of obs: 117, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 1678.47939 189.48902 8.858
bb_colonization1 56.28936 27.59790 2.040
sampling_time1 -48.58482 27.58574 -1.761
precol_cp_c 0.10952 0.36345 0.301
bb_colonization1:sampling_time1 -218.28334 27.58634 -7.913
sampling_time1:precol_cp_c -0.06059 0.04108 -1.475
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.094
samplng_tm1 -0.006 0.040
precol_cp_c -0.004 -0.001 0.007
bb_clnz1:_1 0.006 -0.040 -0.644 -0.007
smplng_1:__ 0.009 -0.061 0.004 0.002 -0.008N:P burnt
# R^2
performance::r2(np_fit_burnt) # R2 for Mixed Models
Conditional R2: 0.733
Marginal R2: 0.348# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit_burnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.591
Conditional ICC: 0.385# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit_burnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: np_ratio
Chisq Df Pr(>Chisq)
bb_colonization 78.4014 1 < 0.00000000000000022 ***
sampling_time 0.8924 1 0.3448
precol_np_c 1.7991 1 0.1798
bb_colonization:sampling_time 0.8073 1 0.3689
sampling_time:precol_np_c 21.2492 1 0.000004033 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(np_fit_burnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.98100045 3.5286402
.sigma 1.98965206 2.5806700
(Intercept) 17.08966961 20.4967022
bb_colonization1 1.79326545 2.8096969
sampling_time1 -0.84477429 0.1716571
precol_np_c -0.02755255 0.5535212
bb_colonization1:sampling_time1 -0.27468208 0.7417494
sampling_time1:precol_np_c 0.09525869 0.2355934# model summary
summary(np_fit_burnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -0.5496
Linear mixed model fit by REML ['blmerMod']
Formula: np_ratio ~ bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_burnt
REML criterion at convergence: 559.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.0286 -0.6498 -0.1108 0.6674 2.9046
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 7.485 2.736
Residual 5.189 2.278
Number of obs: 120, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 18.79319 1.14673 16.389
bb_colonization1 2.30148 0.25992 8.854
sampling_time1 -0.33656 0.25992 -1.295
precol_np_c 0.26298 0.19607 1.341
bb_colonization1:sampling_time1 0.23353 0.25992 0.898
sampling_time1:precol_np_c 0.16543 0.03589 4.610
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.136
samplng_tm1 0.000 0.000
precol_np_c 0.000 0.000 0.000
bb_clnz1:_1 0.000 0.000 -0.600 0.000
smplng_1:__ 0.000 0.000 0.000 0.000 0.000N:P unburnt
# R^2
performance::r2(np_fit_unburnt) # R2 for Mixed Models
Conditional R2: 0.854
Marginal R2: 0.059# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit_unburnt) # Intraclass Correlation Coefficient
Adjusted ICC: 0.844
Conditional ICC: 0.795# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit_unburnt) Analysis of Deviance Table (Type II Wald chisquare tests)
Response: np_ratio
Chisq Df Pr(>Chisq)
bb_colonization 8.1073 1 0.004409 **
sampling_time 3.9313 1 0.047395 *
precol_np_c 0.0002 1 0.988471
bb_colonization:sampling_time 34.0067 1 0.000000005492 ***
sampling_time:precol_np_c 0.1284 1 0.720080
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 95% confident interval
confint(np_fit_unburnt) Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 1.64842715 5.46636294
.sigma 1.58115782 2.05802842
(Intercept) 13.60947431 18.80887503
bb_colonization1 0.24428682 1.10241383
sampling_time1 0.72730751 1.58486209
precol_np_c -0.50505491 0.49601447
bb_colonization1:sampling_time1 -1.70776804 -0.85019469
sampling_time1:precol_np_c -0.05216145 0.07605762# model summary
summary(np_fit_unburnt) Cov prior : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev : -2.5378
Linear mixed model fit by REML ['blmerMod']
Formula: np_ratio ~ bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time + (1 | site:tree)
Data: meta_chem_t1t2_cent_unburnt
REML criterion at convergence: 499.9
Scaled residuals:
Min 1Q Median 3Q Max
-4.0005 -0.5508 0.0111 0.5687 3.3721
Random effects:
Groups Name Variance Std.Dev.
site:tree (Intercept) 17.87 4.227
Residual 3.29 1.814
Number of obs: 117, groups: site:tree, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 16.208832 1.739468 9.318
bb_colonization1 0.675333 0.219492 3.077
sampling_time1 1.156293 0.219352 5.271
precol_np_c -0.004533 0.335174 -0.014
bb_colonization1:sampling_time1 -1.279185 0.219357 -5.832
sampling_time1:precol_np_c 0.011753 0.032796 0.358
Correlation of Fixed Effects:
(Intr) bb_cl1 smpl_1 prcl__ b_1:_1
bb_colnztn1 -0.081
samplng_tm1 -0.005 0.040
precol_np_c -0.004 -0.001 0.006
bb_clnz1:_1 0.005 -0.040 -0.644 -0.006
smplng_1:__ 0.008 -0.064 0.004 0.003 -0.008Session Info
devtools::session_info() ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
setting value
version R version 3.6.2 (2019-12-12)
os macOS High Sierra 10.13.6
system x86_64, darwin15.6.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Stockholm
date 2021-12-02
─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
package * version date lib source
abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0)
assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0)
backports 1.2.1 2020-12-09 [1] CRAN (R 3.6.2)
blme * 1.0-5 2021-01-05 [1] CRAN (R 3.6.2)
boot 1.3-28 2021-05-03 [1] CRAN (R 3.6.2)
broom 0.7.7 2021-06-13 [1] CRAN (R 3.6.2)
bslib 0.2.5.1 2021-05-18 [1] CRAN (R 3.6.2)
cachem 1.0.5 2021-05-15 [1] CRAN (R 3.6.2)
callr 3.7.0 2021-04-20 [1] CRAN (R 3.6.2)
car * 3.0-10 2020-09-29 [1] CRAN (R 3.6.2)
carData * 3.0-4 2020-05-22 [1] CRAN (R 3.6.2)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.0)
cli 2.5.0 2021-04-26 [1] CRAN (R 3.6.2)
cluster 2.1.2 2021-04-17 [1] CRAN (R 3.6.2)
colorspace 2.0-1 2021-05-04 [1] CRAN (R 3.6.2)
cowplot 1.1.1 2020-12-30 [1] CRAN (R 3.6.2)
crayon 1.4.1 2021-02-08 [1] CRAN (R 3.6.2)
curl 4.3.1 2021-04-30 [1] CRAN (R 3.6.2)
data.table 1.14.0 2021-02-21 [1] CRAN (R 3.6.2)
DBI 1.1.1 2021-01-15 [1] CRAN (R 3.6.2)
desc 1.3.0 2021-03-05 [1] CRAN (R 3.6.2)
devtools 2.4.2 2021-06-07 [1] CRAN (R 3.6.2)
digest 0.6.27 2020-10-24 [1] CRAN (R 3.6.2)
dplyr * 1.0.6 2021-05-05 [1] CRAN (R 3.6.2)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 3.6.2)
evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0)
fansi 0.5.0 2021-05-25 [1] CRAN (R 3.6.2)
farver 2.1.0 2021-02-28 [1] CRAN (R 3.6.2)
fastmap 1.1.0 2021-01-25 [1] CRAN (R 3.6.2)
forcats 0.5.1 2021-01-27 [1] CRAN (R 3.6.2)
foreign 0.8-76 2020-03-03 [1] CRAN (R 3.6.0)
fs 1.5.0 2020-07-31 [1] CRAN (R 3.6.2)
generics 0.1.0 2020-10-31 [1] CRAN (R 3.6.2)
ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 3.6.2)
ggpubr 0.4.0 2020-06-27 [1] CRAN (R 3.6.2)
ggsignif 0.6.2 2021-06-14 [1] CRAN (R 3.6.2)
glue 1.4.2 2020-08-27 [1] CRAN (R 3.6.2)
gridExtra 2.3 2017-09-09 [1] CRAN (R 3.6.0)
gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0)
haven 2.4.1 2021-04-23 [1] CRAN (R 3.6.2)
highr 0.9 2021-04-16 [1] CRAN (R 3.6.2)
hms 1.1.0 2021-05-17 [1] CRAN (R 3.6.2)
htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 3.6.2)
insight 0.14.1 2021-05-28 [1] CRAN (R 3.6.2)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 3.6.2)
jsonlite 1.7.2 2020-12-09 [1] CRAN (R 3.6.2)
knitr 1.33 2021-04-24 [1] CRAN (R 3.6.2)
labeling 0.4.2 2020-10-20 [1] CRAN (R 3.6.2)
lattice * 0.20-44 2021-05-02 [1] CRAN (R 3.6.2)
lifecycle 1.0.0 2021-02-15 [1] CRAN (R 3.6.2)
lme4 * 1.1-27 2021-05-15 [1] CRAN (R 3.6.2)
magrittr 2.0.1 2020-11-17 [1] CRAN (R 3.6.2)
maps 3.3.0 2018-04-03 [1] CRAN (R 3.6.0)
MASS 7.3-54 2021-05-03 [1] CRAN (R 3.6.2)
Matrix * 1.3-4 2021-06-01 [1] CRAN (R 3.6.2)
memoise 2.0.0 2021-01-26 [1] CRAN (R 3.6.2)
mgcv 1.8-36 2021-06-01 [1] CRAN (R 3.6.2)
minqa 1.2.4 2014-10-09 [1] CRAN (R 3.6.0)
modelr * 0.1.8 2020-05-19 [1] CRAN (R 3.6.2)
munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0)
nlme 3.1-152 2021-02-04 [1] CRAN (R 3.6.2)
nloptr 1.2.2.2 2020-07-02 [1] CRAN (R 3.6.2)
openxlsx 4.2.3 2020-10-27 [1] CRAN (R 3.6.2)
performance 0.7.2 2021-05-17 [1] CRAN (R 3.6.2)
permute * 0.9-5 2019-03-12 [1] CRAN (R 3.6.0)
pillar 1.6.1 2021-05-16 [1] CRAN (R 3.6.2)
pkgbuild 1.2.0 2020-12-15 [1] CRAN (R 3.6.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0)
pkgload 1.2.1 2021-04-06 [1] CRAN (R 3.6.2)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.0)
processx 3.5.2 2021-04-30 [1] CRAN (R 3.6.2)
ps 1.6.0 2021-02-28 [1] CRAN (R 3.6.2)
purrr 0.3.4 2020-04-17 [1] CRAN (R 3.6.2)
R6 2.5.0 2020-10-28 [1] CRAN (R 3.6.2)
Rcpp 1.0.6 2021-01-15 [1] CRAN (R 3.6.2)
readxl 1.3.1 2019-03-13 [1] CRAN (R 3.6.0)
remotes 2.4.0 2021-06-02 [1] CRAN (R 3.6.2)
rio 0.5.26 2021-03-01 [1] CRAN (R 3.6.2)
rlang 0.4.11 2021-04-30 [1] CRAN (R 3.6.2)
rmarkdown 2.8 2021-05-07 [1] CRAN (R 3.6.2)
rprojroot 2.0.2 2020-11-15 [1] CRAN (R 3.6.2)
rstatix 0.7.0 2021-02-13 [1] CRAN (R 3.6.2)
sass 0.4.0 2021-05-12 [1] CRAN (R 3.6.2)
scales 1.1.1 2020-05-11 [1] CRAN (R 3.6.2)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0)
stringi 1.6.2 2021-05-17 [1] CRAN (R 3.6.2)
stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0)
testthat 3.0.2 2021-02-14 [1] CRAN (R 3.6.2)
tibble 3.1.2 2021-05-16 [1] CRAN (R 3.6.2)
tidyr * 1.1.3 2021-03-03 [1] CRAN (R 3.6.2)
tidyselect 1.1.1 2021-04-30 [1] CRAN (R 3.6.2)
usethis 2.0.1 2021-02-10 [1] CRAN (R 3.6.2)
utf8 1.2.1 2021-03-12 [1] CRAN (R 3.6.2)
vctrs 0.3.8 2021-04-29 [1] CRAN (R 3.6.2)
vegan * 2.5-7 2020-11-28 [1] CRAN (R 3.6.2)
viridis * 0.6.1 2021-05-11 [1] CRAN (R 3.6.2)
viridisLite * 0.4.0 2021-04-13 [1] CRAN (R 3.6.2)
withr 2.4.2 2021-04-18 [1] CRAN (R 3.6.2)
xfun 0.23 2021-05-15 [1] CRAN (R 3.6.2)
yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.0)
zip 2.2.0 2021-05-31 [1] CRAN (R 3.6.2)
[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library