Derived Root Depths - GCREW CO2 by Competition 1999
James R Holmquist
3/15/2021
# co2_community_roots
co2_community_root_cores <- read_csv("../1999 CO2xCommunity Root Biomass Data Assembled.csv",
na = c("", "NA", -99))## Parsed with column specification:
## cols(
## community = col_character(),
## plot = col_double(),
## treatment = col_character(),
## core = col_double(),
## depth_center = col_double(),
## target_width = col_double(),
## volume = col_double(),
## washer = col_character(),
## plucker = col_character(),
## category = col_character(),
## mass = col_double(),
## seeds = col_double()
## )Top Level Curation
First, I'm making sure that there are no redundant entries. This was an issue with an earlier version. I used R (R Core Team 2019) and tidyverse (Wickham 2017) packages for curation and graphing.
# First curate the CO2 x Community Data
co2_by_community_n <- co2_community_root_cores %>%
# Group by community, plot, treamtement, core, depth increment, and category
group_by(community, plot, treatment, category, core, depth_center) %>%
# Get the mean biomass
summarise(n = n()) %>%
arrange(n)
head(co2_by_community_n)## # A tibble: 6 x 7
## # Groups: community, plot, treatment, category, core [1]
## community plot treatment category core depth_center n
## <chr> <dbl> <chr> <chr> <dbl> <dbl> <int>
## 1 SC 1 E DRT 1 2.5 1
## 2 SC 1 E DRT 1 5 1
## 3 SC 1 E DRT 1 7.5 1
## 4 SC 1 E DRT 1 10 1
## 5 SC 1 E DRT 1 12.5 1
## 6 SC 1 E DRT 1 15 1tail(co2_by_community_n)## # A tibble: 6 x 7
## # Groups: community, plot, treatment, category, core [1]
## community plot treatment category core depth_center n
## <chr> <dbl> <chr> <chr> <dbl> <dbl> <int>
## 1 SP 14 E WRH 2 15 1
## 2 SP 14 E WRH 2 21.2 1
## 3 SP 14 E WRH 2 32.5 1
## 4 SP 14 E WRH 2 52.5 1
## 5 SP 14 E WRH 2 72.5 1
## 6 SP 14 E WRH 2 92.5 1Then I did some quick renaming, rearranging of variables and filtering.
co2_by_community <- co2_community_root_cores %>%
# create top and bottom depth increments
mutate(top = depth_center - (target_width)/2,
bottom = depth_center + (target_width)/2) %>%
# Filter out litter and seeds
filter(! category %in% c("LITT", "SCseeds")) %>%
ungroup() %>%
# Rename abbreviated samples to make the more readable
mutate(category = recode(category,
"LRT"="light-colored roots",
"RRT"="red roots",
"DRT"="dark roots",
"WRH"="white rhizome",
"RRH"="red rhizome",
"RHC"="rhizomotous culms",
"NRHC"="non-rhizomotous culms"),
community = recode(community, "SC"= "Schoenoplectus americanus",
"SP"="Spartina patens"),
treatment = recode(treatment, "E"="Elevated",
"A" = "Ambient")
)
biomass_com_core_cat <- ggplot(data = co2_by_community,
aes(x = bottom,
y = mass,
color = category)) +
geom_point(alpha=0.66) +
geom_line(aes(group = paste(community, plot, treatment, category, core, sep="-")),
lty = 2,
alpha=0.66) +
facet_wrap(community~treatment) +
xlab("Depth (cm)") +
ylab("Biomass (g)") +
scale_x_reverse() +
coord_flip() +
theme(legend.position = "right",
legend.title = element_blank()) +
guides(col = guide_legend(ncol = 1))
(biomass_com_core_cat)
Scale Up Measured Depth Increments to Representative Depth Increments
The measurements for each core are not all continuous. Near the top of the core measurements are continuous and 2.5 cm wide, but deeper, they are discontinous, 5 cm wide, representing wider depth increments. For biomass measurements, I scaled weights by the ratio of their measured width to their representative width. For each core, I calculated the cumulative sum of biomass for each depth increment, and normalized it by the total sum of biomass. I then rescaled this normalized biomass to equal 0.999, or 99.9% of the total biomass. I am assuming that the shape of the root zone is exponential, and that the core captures nearly 100% (%99.9) of the total live roots of the plants.
co2_by_community_scaled <- co2_by_community %>%
group_by(community, plot, core, treatment, depth_center, target_width, top, bottom, volume) %>%
summarise(total_mass = sum(mass)) %>%
ungroup() %>%
group_by(community, plot, core, treatment) %>%
mutate(increment_n = 1:n()) %>%
mutate(representative_top = ifelse(increment_n == 1,
0,
(lag(bottom)-top)/2+top),
representative_bottom = ifelse(increment_n == max(increment_n),
max(bottom),
(bottom-lead(top))/2+lead(top)),
representative_volume = (representative_bottom-representative_top)*(5.1/2)^2*pi) %>%
ungroup() %>%
mutate(scale_factor = representative_volume/volume,
scaled_mass = total_mass * scale_factor) %>%
group_by(community, plot, core, treatment) %>%
mutate(cumSumBmass = cumsum(scaled_mass) ,
cumSumBmassNormalized = cumSumBmass / max(cumSumBmass)*0.999
)
cumSumBmassPlot <- ggplot(data = co2_by_community_scaled, aes(x = representative_bottom, y = cumSumBmassNormalized)) +
geom_line(aes(group = paste(community, plot, treatment, core, sep="-")), color="grey") +
geom_segment(aes(xend = representative_top, yend=cumSumBmassNormalized), color="darkred") +
facet_wrap(community~treatment) +
xlab("Depth (cm)") +
ylab("Cumulative Biomass Sum (proportion)") +
scale_x_reverse() +
coord_flip()
(cumSumBmassPlot)
Calculate Total Biomass By Core
To calculate total biomass by core I summarised the scaled biomass for each core and divided by the area of the 5.1 cm (=0.051 m) diameter core (\(A = \pi(0.051/2)^2\)).
CO2_community_bg_summary_by_core <- co2_by_community_scaled %>%
group_by(community, plot, treatment, core) %>%
summarise(BGB_total = sum(scaled_mass)/((0.051/2)^2*pi))
bgbTotalPlot <- ggplot(CO2_community_bg_summary_by_core, aes(x=BGB_total, color = community, lty = treatment)) +
geom_density() +
xlab(expression(paste("Total Below Ground Biomass (g m"^"-2",")", sep=""))) +
theme(legend.title = element_blank())
(bgbTotalPlot)
Model 95% max root depth
To model the root zone, for each core, I treated the normalized, cumulative belowground biomass (\(C\)) as a function of depth (\(d\)) using the cumulative density function (cdf) for an exponential distribution.
\[ C = 1 - e^{-\lambda d}\]
I solved for \(\lambda\) using the non-linear solver (nls) in R.
I then extrapolated the maximum rooting depth of 95% of the root mass \(D\), using the quantile for the exponential distribution, using \(\lambda\) and a probability (\(p\)) of 0.95. The max root depth of 95% of the biomass is how exponential root distribution is parameterized in some versions of the marsh equilibrium model (MEM) (Schile et al. 2014).
\[ D = log(1-p)/\lambda \]
# For each unique community, plot, core
for (i in 1:nrow(CO2_community_bg_summary_by_core)) {
depth_series_i <- filter(co2_by_community_scaled,
community == CO2_community_bg_summary_by_core$community[i] &
plot == CO2_community_bg_summary_by_core$plot[i] &
core == CO2_community_bg_summary_by_core$core[i]) %>%
filter(complete.cases(.))
# Fit the CDF for an exponential distribution
# {1 - exp(-lambda * x)
f1 <- formula("cumSumBmassNormalized~(1-exp(-lambda * representative_bottom))")
model_i <- nls(f1, data = depth_series_i, start = list(lambda = 0.05))
lambda_i <- coef(model_i)[1]
# Then get the quantile for p = 0.95
# - (log(1-p)/lambda)
max_rooting_depth_i <- - (log(1-0.95)/lambda_i)
# Save and add to file
CO2_community_bg_summary_by_core$max_rooting_depth[i] <- max_rooting_depth_i
}
MaxRootingDepthPlot <- ggplot(CO2_community_bg_summary_by_core,
aes(x=max_rooting_depth, color = community, lty = treatment)) +
geom_density() +
xlab("95% Max Rooting Depth (cm)") +
theme(legend.title = element_blank())
(MaxRootingDepthPlot) 
Run Stats 1
I used the the lme4 package (Bates et al. 2015) to fit hierarchical linear models, and the lmerTest package (Kuznetsova, Brockhoff, and Christensen 2017) to test the significance of random effects. I assumed normality for response variables based on my visual inspection of their density plots. I fit linear models with the experimental treatement as the independent variable and total belowground biomass or max rooting depth as the dependent variable. I treated plant community as a random effect.
Experimental Effects on BG total, Heirarchical for Species
bg_total_test <- lmer(BGB_total ~ treatment + (1|community), data = CO2_community_bg_summary_by_core)
summary(bg_total_test)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BGB_total ~ treatment + (1 | community)
## Data: CO2_community_bg_summary_by_core
##
## REML criterion at convergence: 640
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.72275 -0.57512 -0.01722 0.50231 2.74737
##
## Random effects:
## Groups Name Variance Std.Dev.
## community (Intercept) 1228554 1108.4
## Residual 945347 972.3
## Number of obs: 40, groups: community, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2725.514 813.354 1.075 3.351 0.1705
## treatmentElevated 767.126 307.465 37.000 2.495 0.0172 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## trtmntElvtd -0.189ranova(bg_total_test)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## BGB_total ~ treatment + (1 | community)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 -319.99 647.98
## (1 | community) 3 -328.24 662.49 16.509 1 4.843e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Total below ground biomass was significantly (p=0.0172) and positively affected by the experimental treatement. Spartina patens had significantly less below ground biomass than Scirpus olneyi (p<0.0001).
root_depth_test <- lmer(max_rooting_depth ~ treatment + (1|community), data = CO2_community_bg_summary_by_core)## boundary (singular) fit: see ?isSingularsummary(root_depth_test)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: max_rooting_depth ~ treatment + (1 | community)
## Data: CO2_community_bg_summary_by_core
##
## REML criterion at convergence: 321.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8110 -0.6985 -0.1787 0.5552 2.9876
##
## Random effects:
## Groups Name Variance Std.Dev.
## community (Intercept) 0.0 0.00
## Residual 238.5 15.44
## Number of obs: 40, groups: community, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 50.403 3.453 38.000 14.60 <2e-16 ***
## treatmentElevated -2.782 4.883 38.000 -0.57 0.572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## trtmntElvtd -0.707
## convergence code: 0
## boundary (singular) fit: see ?isSingularranova(root_depth_test)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## max_rooting_depth ~ treatment + (1 | community)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 -160.93 329.85
## (1 | community) 3 -160.93 327.85 0 1 1While in the plots, it appeared that maximum root depth was shallower for the experimental treatment, this difference was not significant (p=0.572), and maximum rooting depth was indistinguishable between species (p=1).
Works Cited
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. doi:10.18637/jss.v067.i01.
Kuznetsova, Alexandra, Per B. Brockhoff, and Rune H. B. Christensen. 2017. “lmerTest Package: Tests in Linear Mixed Effects Models.” Journal of Statistical Software 82 (13): 1–26. doi:10.18637/jss.v082.i13.
R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Schile, Lisa M, John C Callaway, James T Morris, Diana Stralberg, V Thomas Parker, and Maggi Kelly. 2014. “Modeling Tidal Marsh Distribution with Sea-Level Rise: Evaluating the Role of Vegetation, Sediment, and Upland Habitat in Marsh Resiliency.” PloS One 9 (2). Public Library of Science: e88760.
Wickham, Hadley. 2017. Tidyverse: Easily Install and Load the ’Tidyverse’. https://CRAN.R-project.org/package=tidyverse.