https://github.com/tom-n-walker/uphill-plants-soil-carbon
Tip revision: 05e1672664e587a3808d326d6f78d29ca5e44c01 authored by Tom Walker on 07 September 2021, 14:31:53 UTC
Added all analyses, streamlined pipelines and improved statistical models. Removed data export (all done via drake cache).
Added all analyses, streamlined pipelines and improved statistical models. Removed data export (all done via drake cache).
Tip revision: 05e1672
field_explain_soil_carbon_loss.R
################################################################################
#### Project: Lowland plant migrations alpine soil C loss
#### Title: Explaining field soil C loss
#### Author: Tom Walker (thomas.walker@usys.ethz.ch)
#### Date: 26 May 2021
#### ---------------------------------------------------------------------------
#### PROLOGUE ------------------------------------------------------------------
## Options ----
# remove objects from global environment
rm(list = ls())
# set random seed
set.seed(3)
# R session options (no factors, bias against scientific #s)
options(
stringsAsFactors = F,
scipen = 6
)
## Libraries ----
# standard library set
library(nlme)
library(emmeans)
library(tidyverse)
library(vegan)
#### DATA ----------------------------------------------------------------------
## Load from drake plan ----
# main targets
fieldData <- drake::readd(field_data_subset)
fieldTraits <- drake::readd(trait_data)
ghTraits <- drake::readd(gh_plants)
# separate plant datasets
westFocalsLong <- fieldData$focals$all_long[[2]]
centFocalsLong <- fieldData$focals$all_long[[1]]
westFocalsWide <- fieldData$focals$all_wide[[2]]
centFocalsWide <- fieldData$focals$all_wide[[1]]
responses <- fieldData$plantsRR
# remove block for wide datasets (ease of modelling)
westFocalsWide <- westFocalsWide %>%
select(-block)
centFocalsWide <- centFocalsWide %>%
select(-block)
# remove Dactylis, Viola, Plantago and Trifolium from central alps (very rare)
onlyThere <- filter(centFocalsLong, cover > 0)
table(onlyThere$species)
centFocalsLong <- centFocalsLong %>%
filter(species != "Dactylis.glomerata") %>%
filter(species != "Plantago.lanceolata") %>%
filter(species != "Trifolium.montanum") %>%
filter(species != "Viola.hirta")
toRemove <- colnames(centFocalsWide) %in% c("Dactylis.glomerata", "Plantago.lanceolata", "Trifolium.montanum", "Viola.hirta")
centFocalsWide <- centFocalsWide[, !toRemove]
#### BROAD TRAIT DIFFERENCES ---------------------------------------------------
## Field traits ----
# format traits
ftOnly <- fieldTraits %>%
select(leaf_area:SLA) %>%
mutate(leaf_area = log10(leaf_area)) %>%
mutate(plant_height = log10(plant_height)) %>%
mutate(seed_mass = log10(seed_mass))
# PCA
ftPCA <- prcomp(ftOnly, center = T, scale = T)
# PERMANOVA
adonis2(ftOnly ~ is_focal, fieldTraits, method = "gower")
## Glasshouse traits ----
# format traits
ghtOnly <- ghTraits %>%
select(rAGB.mg:gsmax)
# PCA
ghtPCA <- prcomp(ghtOnly, center = T, scale = T)
# PERMANOVA
adonis2(ghtOnly ~ Origin, ghTraits)
## Biplots ----
par(mfrow = c(1,1))
biplot(ftPCA)
biplot(ghtPCA)
#### ANALYSE IDENTITY EFFECTS - RANDOM INTERCEPT -------------------------------
## West Alps ----
# build model
m1 <- lme(
soilCloss ~ cover,
random = ~ 1 | species,
data = westFocalsLong,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r1 <- residuals(m1, type = "normalized")
par(mfrow = c(1, 3))
plot(r1 ~ fitted(m1))
plot(r1 ~ westFocalsLong$cover)
hist(r1)
# test main effects
m1a <- update(m1, ~.- cover)
anova(m1, m1a)
## Central Alps ----
# build model
m2 <- lme(
soilCloss ~ cover,
random = ~ 1 | species,
data = centFocalsLong,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r2 <- residuals(m2, type = "normalized")
par(mfrow = c(1, 3))
plot(r2 ~ fitted(m2))
plot(r2 ~ centFocalsLong$cover)
hist(r2)
# test main effects
m2a <- update(m2, ~.- cover)
anova(m2, m2a)
#### ANALYSE IDENTITY EFFECTS - SPECIES SEPARATE -------------------------------
## West Alps ----
# build full and null models
m3x <- lm(soilCloss ~ 1., westFocalsWide)
m3 <- lm(soilCloss ~ ., westFocalsWide)
# diagnose model fit
par(mfrow = c(1, 4))
plot(m3)
# overall model fit
anova(m3x, m3)
# individual main effects
anova(m3)
## Central Alps ----
# build full and null models
m4x <- lm(soilCloss ~ 1., centFocalsWide)
m4 <- lm(soilCloss ~ ., centFocalsWide)
# diagnose model fit
par(mfrow = c(1, 4))
plot(m4)
# overall model fit
anova(m4x, m4)
# individual main effects
anova(m4)
#### ANALYSE TOTAL LOWLAND COVER -----------------------------------------------
# build model
m5 <- lm(
Csoil ~ focal_bio + bkgnd_bio + vege_bio + site,
data = responses
)
# diagnose
plot(m5)
# main effects
anova(m5)
#### ANALYSE BACKGROUND COMMUNITY COMPOSITION ----------------------------------
## Analyse treatment effects on raw NMDS 1 scores ----
# build model
m6 <- lme(
bgPCC1 ~ treatment * site,
random = ~ 1 | site_block,
data = fieldData$plantsFull,
method = "ML"
)
# diagnose model
r6 <- residuals(m6, type = "normalized")
par(mfrow = c(1, 3))
plot(r6 ~ fitted(m6))
boxplot(r6 ~ fieldData$plantsFull$treatment)
hist(r6)
# test main effects
m6a <- update(m6, ~.- treatment:site)
anova(m6, m6a)
anova(m6a, update(m6a, ~.- treatment))
anova(m6a, update(m6a, ~.- site))
## Analyse treatment effects on raw NMDS 2 scores ----
# build model
m7 <- lme(
bgPCC2 ~ treatment * site,
random = ~ 1 | site_block,
data = fieldData$plantsFull,
method = "ML"
)
# diagnose model
r7 <- residuals(m7, type = "normalized")
par(mfrow = c(1, 3))
plot(r7 ~ fitted(m7))
boxplot(r7 ~ fieldData$plantsFull$treatment)
hist(r7)
# test main effects
m7a <- update(m7, ~.- treatment:site)
anova(m7, m7a)
anova(m7a, update(m7a, ~.- treatment))
anova(m7a, update(m7a, ~.- site))
## Analyse NMDS change effects on soil C loss ----
# remove one very large value for central alps (no effect on result)
responses[responses$bgPCC2 > 30, "bgPCC2"] <- NA
# build model
m8 <- lm(
Csoil ~ bgPCC1 + bgPCC2 + site,
data = responses
)
# diagnose
par(mfrow = c(1, 4))
plot(m8)
# main effects
anova(m8)