Nothing
inst/doc/JEMBE.R
In NitrogenUptake2016: Data and Source Code From: Nitrogen Uptake and Allocation Estimates for Spartina Alterniflora and Distichlis Spicata
## ----setup, include=FALSE, echo=FALSE------------------------------------
if(!require(knitr)){
install.packages("knitr", repos='http://cran.us.r-project.org')
}
if(!require(rmarkdown)){
install.packages("rmarkdown", repos='http://cran.us.r-project.org')
}
if(!require(plyr)){
install.packages("plyr", repos='http://cran.us.r-project.org')
}
if(!require(reshape2)){
install.packages("reshape2", repos='http://cran.us.r-project.org')
}
if(!require(ggplot2)){
install.packages("ggplot2", repos='http://cran.us.r-project.org')
}
if(!require(scales)){
install.packages("scales", repos='http://cran.us.r-project.org')
}
if(!require(MuMIn)){
install.packages("MuMIn", repos='http://cran.us.r-project.org')
}
if(!require(rsq)){
install.packages("rsq", repos='http://cran.us.r-project.org')
}
# then load the package:
library(knitr)
library(rmarkdown)
library(plyr)
library(reshape2)
library(ggplot2)
library(scales)
library(MuMIn)
library(rsq)
library(NitrogenUptake2016)
knitr::opts_chunk$set(echo = TRUE, comment=NA)
## ---- echo = FALSE, include=FALSE----------------------------------------
# if the NitrogenUptake2016 package isn't installed, use devtools to do so:
# devtools::install_github("troyhill/NitrogenUptake2016", build_vignettes = TRUE)
# set some constants
todaysDate <- substr(as.character(Sys.time()), 1, 10)
core.area <- pot.m2 <- 0.00801185 # mesocosm surface area (m2)
top.vol <- core.area * 0.05 * 1e6 # cm3 in core: top 5 cm only
pointSize <- 2 # for ggplot graphics
pd <- pd2 <- position_dodge(1.2)
pd3 <- position_dodge(0.8)
grayColor <- 0.55
fig2Col <- "gray55"
## ----allometry, include=FALSE, echo=FALSE--------------------------------
# Allometry ----------------------------------------------------------------
##### Establish relationships following Lu et al. 2016
CSP <- plyr::dlply(allometry, c("spp"), bCM)
CSP.coef <- plyr::ldply(CSP, stats::coef)
CSP.coef$lam <- plyr::ddply(allometry, c("spp"), function(df) bCM(df, lam.only = TRUE))[, "V1"]
CSP.coef$BIC <- plyr::ldply(CSP, stats::BIC)[, "V1"]
CSP.coef$rsq <- plyr::ldply(CSP, MuMIn::r.squaredGLMM)[, "R2m"]
#CSP.coef
## ----applying allometry to estimate plant masses, include=FALSE, echo=FALSE----
### estimate plant masses from allometry
for (i in 1:nrow(stemHeights)) {
if (stemHeights$species[i] %in% "DS") {
stemHeights$mass[i] <- (stemHeights$height_cm[i] * CSP.coef[1, 3] + CSP.coef[1, "(Intercept)"])^(1/CSP.coef$lam[1]) # DISP
} else if (stemHeights$species[i] %in% "SA") {
stemHeights$mass[i] <- (stemHeights$height_cm[i] * CSP.coef[2, 3] + CSP.coef[2, "(Intercept)"])^(1/CSP.coef$lam[2]) # SPAL
}
}
## ----Allometry validation: comparison between predicted and harvested masses, include=FALSE, echo=FALSE----
###
### compare predicted (allometry) and observed masses at harvest
###
# get live + dead biomass for each pot at each date
dat.ld <- plyr::ddply(stemHeights, plyr::.(core_num, species, date), plyr::summarise,
live = sum(mass[dead_live %in% "L"], na.rm = TRUE),
dead = sum(mass[dead_live %in% "D"], na.rm = TRUE)
)
# add cohort number
dat.ld$day <- as.POSIXct(as.character(dat.ld$date), format = "%y%m%d", origin = "1960-01-01")
dat.ld$cohort <- NA
for (i in 1:length(unique(dat.ld$day))) {
uniqueIDS <- unique(dat.ld$core_num[dat.ld$day == sort(unique(dat.ld$day))[i]])
if (length(uniqueIDS) == 3) {
dat.ld$cohort[dat.ld$core_num %in% uniqueIDS] <- ifelse(sum(is.finite(dat.ld$cohort)) == 0, 1, max(dat.ld$cohort, na.rm = TRUE) + 1)
}
}
dat.ld$cohort[is.na(dat.ld$cohort)] <- 0
dat.ld$cohort[dat.ld$cohort == 5] <- 4
dat.ld$year <- 2016
dat.ld$dead <- 0
dat.ld$pot2 <- paste0(dat.ld$core_num, "-", dat.ld$species)
### allometry estimates
for (i in 1:length(unique(dat.ld$pot2))) {
targPot <- unique(dat.ld$pot2)[i]
maxDate <- max(dat.ld$day[dat.ld$pot2 %in% targPot])
tempDat <- dat.ld[(dat.ld$pot2 %in% targPot) & (dat.ld$day == maxDate), ]
if (i == 1) {
dat.ld.sub <- tempDat
}
if (i > 1) {
dat.ld.sub <- rbind(dat.ld.sub, tempDat)
}
}
dat.ld.sub$id <- paste0(dat.ld.sub$species, "-", dat.ld.sub$core_num)
### observed values
obs <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("leaf 1", "leaf 2", "leaf 3", "leaf 4", # "dead leaf",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
"stems"), ], plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T)
)
obs$core_num <- as.character(as.integer(substr(obs$new.core.id, 3, 4)))
obs$core_num[as.integer(obs$core_num) > 12] <- paste0(obs$species[as.integer(obs$core_num) > 12], as.integer(obs$core_num[as.integer(obs$core_num) > 12]) - 12, "_T0")
obs$id <- paste0(obs$species, "-", obs$core_num)
obs <- plyr::join_all(list(obs, dat.ld.sub[, c(4,9:10)]), by = "id")
names(obs)[which(names(obs) %in% "live")] <- "allom.est"
obs$diff <- (obs$allom.est - obs$g) # magnitude accuracy of allometry prediction (g)
obs$diff.pct <- (obs$allom.est - obs$g) / obs$g # percent accuracy of allometry prediction
obs$species2 <- ifelse(obs$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
summary(obs$diff)
## ----Figure 1 (Data In Brief), fig.width = 6, fig.height = 4, include=FALSE, echo=FALSE----
# Figure 1 from Hill et al. (Data In Brief) - predicted vs obs biomass ---------------------------------------
ggplot(obs, aes(x = allom.est / pot.m2, y = g / pot.m2, colour = species2, shape = species2)) +
geom_point(size = pointSize) + theme_classic() %+replace% theme(legend.title = element_blank()) +
labs(x = expression("Predicted biomass (allometry; g"%.%m^-2~")"), y = expression("Measured biomass (g "%.%m^-2~")")) +
ylim(0, 650) + xlim(0, 650) +
scale_colour_grey(start = 0.5, end = 0.1, name = "", breaks = c(unique(obs$species2)[1], unique(obs$species2)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(obs$species2)[1], unique(obs$species2)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
theme(legend.position = c(0.3, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA)) +
geom_abline(slope = 1, intercept = 0, linetype = 2)
# ggsave(paste0("C:/RDATA/greenhouse/output/GRO/FigureS1_", todaysDate, ".eps"), width = 90, height= 90, units = "mm", dpi = 400)
## ----Relative growth rates, include=FALSE, echo=FALSE--------------------
# Relative growth rates ---------------------------------------------------
stemHeights$RGR <- stemHeights$lateStartDate <- NA
for(i in 1:length(unique(stemHeights$id[stemHeights$dead_live %in% "L"]))) { # 830 plants
targPlant <- unique(stemHeights$id[stemHeights$dead_live %in% "L"])[i]
subDat <- stemHeights[stemHeights$id %in% targPlant, ]
if (nrow(subDat) > 1) {
subDat$mass2 <- c(NA, subDat$mass[1:(nrow(subDat)-1)])
subDat$day2 <- c(subDat$day[1], subDat$day[1:(nrow(subDat)-1)])
subDat$RGR <- (log(subDat$mass) - log(subDat$mass2)) / as.numeric(difftime(subDat$day, subDat$day2, units = "days"))
### add RGRs to dat.live object
stemHeights$RGR[(which(rownames(stemHeights) %in% rownames(subDat)))] <- subDat$RGR
if (!is.na(subDat$height[1]) & (subDat$height[1] < 6)) {
stemHeights$lateStartDate[which(rownames(stemHeights) %in% rownames(subDat))] <- 1 # indicates if high-growth period is captured
}
}
}
stemHeights$RGR[!is.finite(stemHeights$RGR)] <- NA
# mean growth rate for each plant, then each mesocosm
mean.rgr <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(id, new.core.id, species), plyr::summarise, meanRGR = mean(RGR, na.rm = TRUE))
mean.rgr.pot <- plyr::ddply(mean.rgr, plyr::.(new.core.id, species), plyr::summarise, RGR = mean(meanRGR, na.rm = TRUE))
mean.rgr.pot <- mean.rgr.pot[1:24, ]
## ----stem density and primary production, include=FALSE, echo=FALSE------
# Stem density and biomass ------------------------------------------------
### summarize by pot and session
dat.live <- stemHeights[stemHeights$dead_live %in% "L", ]
ddHgt2 <- plyr::ddply(dat.live, plyr::.(date, day, core_num, species), plyr::numcolwise(sum, na.rm = TRUE))
ddHgt2$day <- as.POSIXct(ddHgt2$day, origin = "1960-01-01")
ddHgt2$cohort <- NA
for (i in 1:length(unique(ddHgt2$date))) {
uniqueIDS <- unique(ddHgt2$core_num[ddHgt2$date == sort(unique(ddHgt2$date))[i]])
if (length(uniqueIDS) == 3) {
ddHgt2$cohort[ddHgt2$core_num %in% uniqueIDS] <- i
}
}
ddHgt2$cohort <- (ddHgt2$cohort - 1) / 2
ddHgt2$cohort[is.na(ddHgt2$cohort)] <- 0
ddHgt4 <- plyr::ddply(ddHgt2, plyr::.(day, cohort, species), plyr::numcolwise(mean, na.rm = TRUE))
ddHgt4.se <- plyr::ddply(ddHgt2, plyr::.(day, cohort, species), plyr::numcolwise(se))
names(ddHgt4.se) <- paste0(names(ddHgt4.se), ".se")
ddHgt4 <- cbind(ddHgt4, ddHgt4.se[, c(5, 7, 9)])
ddHgt4$session <- as.POSIXct(ddHgt4$day, origin = "1960-01-01")
ddHgt4$species2 <- ifelse(ddHgt4$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
####
#### Calculate Milner-Hughes NAPP
####
test <- nappCalc2(dataset = dat.ld[dat.ld$cohort > 0, ], liveCol = "live", siteCol = "pot2", timeCol = "day", summarize = "TRUE", annualReset = "FALSE")
napp <- test$summary
napp$pot2 <- as.character(napp$site)
napp$core_num <- as.numeric(sapply(strsplit(napp$pot2, "-"), "[[", 1))
napp$species <- sapply(strsplit(napp$pot2, "-"), "[[", 2)
napp$cohort <- dat.ld$cohort[match(napp$pot2, dat.ld$pot2)]
napp$new.core.id <- ifelse(napp$core_num > 9, paste0(napp$species, napp$core_num), paste0(napp$species, "0", napp$core_num))
### build object to merge with master: average RGR, stem characteristics
# mean growth rate for each plant, then each mesocosm
mean.rgr <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(id, new.core.id, species), plyr::summarise, meanRGR = mean(RGR, na.rm = TRUE))
mean.rgr.pot <- plyr::ddply(mean.rgr, plyr::.(new.core.id, species), plyr::summarise, rgr = mean(meanRGR, na.rm = TRUE))
mean.rgr.pot <- mean.rgr.pot[1:24, ]
mean.dens <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(new.core.id, date), plyr::summarise, stem.density = sum(!is.na(height_cm)))
mean.dens.pot <- plyr::ddply(mean.dens, plyr::.(new.core.id), plyr::summarise, dens = mean(stem.density / pot.m2, na.rm = TRUE)) # mean of all stem density msrmts for a pot
mean.dens.pot <- mean.dens.pot[1:24, ]
rgr.stems <- plyr::join_all(list(mean.rgr.pot, mean.dens.pot), by = "new.core.id")
## ----Denitrification enzyme activity and ambient N2O flux, include = FALSE, echo = FALSE, message = FALSE----
# Denitrification enzyme activity -----------------------------------------
dea$species <- substr(dea$pot, 1, 2)
dea$species2 <- ifelse(dea$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
dea$DEA.m2 <- dea$DEA * dea$bd_gcm3 * (top.vol / core.area) # * 5 * 1e4 # [flux per gram] * [g/cm3] * [5 cm depth] * [cm2/m2]
dea$IV.m2 <- dea$IV * dea$bd_gcm3 * (top.vol / core.area) # * 5 * 1e4
dea$gap.g <- dea$DEA - dea$IV
dea$gap.m2 <- dea$DEA.m2 - dea$IV.m2
dd.dea <- plyr::ddply(dea, plyr::.(species, species2), plyr::summarise,
flux.g = mean(DEA),
flux.m2 = mean(DEA.m2),
flux.g.se = se(DEA),
flux.m2.se = se(DEA.m2)
)
dd.dea$type <- "DEA"
dd.iv <- plyr::ddply(dea, plyr::.(species, species2), plyr::summarise,
flux.g = mean(IV),
flux.m2 = mean(IV.m2),
flux.g.se = se(IV),
flux.m2.se = se(IV.m2)
)
dd.iv$type <- "GHG"
dd.dea <- rbind(dd.dea, dd.iv)
### stats on DEA/IV N2O flux
dea.m <- reshape2::melt(dea, id.vars = c("species", "species2", "pot"), measure.vars = c(2, 3, 8:11))
dea.m$ID <- paste0(dea.m$species, "-", dea.m$variable)
# dea.g.aov <- stats::aov(value ~ factor(variable) * factor(species), data = dea.m[(dea.m$variable %in% c("IV.m2", "DEA.m2")), ])
# car::Anova(dea.g.aov, type = 2, singular = TRUE) # difference between methods, but not species
## ----15N spike, include=FALSE, echo=FALSE--------------------------------
# how much 15N was added to each pot?
KNO3 <- 2.24047 # grams of KNO3 added
soln <- 2027.5 # final volume of spike solution (mls)
KNO3_mw <- 102.1032 # g/mol
N_mw <- 15 # gN/mol KNO3
mls <- 12 # mls added per mesocosm
spike <- KNO3 * (N_mw / KNO3_mw) / soln * mls # total grams of 15N added to each mesocosm
# calculate atom pct from per mil values
CN_mass_data$n_ap <- ap(CN_mass_data$d15n)
## ----15N background, include=FALSE, echo=FALSE---------------------------
# get background 15N AP
bkd <- plyr::ddply(CN_mass_data[as.character(CN_mass_data$time) %in% "t0", ], plyr::.(species, pool_label), plyr::summarise,
n15 = mean(n_ap, na.rm = T),
se.n15 = se(n_ap)
)
## ----15N corrections, include=FALSE, echo=FALSE--------------------------
# subtract background to get 15N excess (if result is negative, 15Nxs = 0)
CN_mass_data$n15xs <- as.numeric(NA)
for(i in 1:nrow(CN_mass_data)) {
if(!is.na(CN_mass_data$species[i])) {
spp <- CN_mass_data$species[i]
pool <- CN_mass_data$pool_label[i]
a <- bkd[(bkd$species %in% spp) & (bkd$pool_label %in% pool), "n15"]
if((is.numeric(a) & (length(a) > 0))) { # excludes all pools we don't have bkd for
difference <- CN_mass_data$n_ap[i] - a # NOT decimal fraction
ifelse (!is.na(difference) & (difference > 0), CN_mass_data$n15xs[i] <- difference, CN_mass_data$n15xs[i] <- 0)
} else {
difference <- CN_mass_data$n_ap[i] - 0.370317701 # assume bkgd of 11 d15N
ifelse ((difference > 0) & !is.na(difference), CN_mass_data$n15xs[i] <- difference, CN_mass_data$n15xs[i] <- 0)
}
}
}
### Additional correction to belowground live biomass, subtracting 15Nxs of dead MOM to account for sorption
CN_mass_data$n15xs_MOM <- CN_mass_data$n15xs # n15xs still useful for total recoveries
for(i in 1:nrow(CN_mass_data)) {
if (CN_mass_data$time[i] %in% paste0("t", 1:4)) {
if(!is.na(CN_mass_data$depth_top[i])) {
spp <- CN_mass_data$species[i]
coreID <- CN_mass_data$new.core.id[i]
pool <- CN_mass_data$pool_label[i]
depth_bottom <- CN_mass_data$depth_bottom[i]
a <- CN_mass_data[(CN_mass_data$new.core.id %in% coreID) & (CN_mass_data$pool_label %in% paste0("dead biomass ", depth_bottom, "cm")) &
(CN_mass_data$depth_bottom == depth_bottom), "n15xs"]
if((is.numeric(a) & (length(a) > 0))) { # excludes all pools we don't have bkd for
difference <- CN_mass_data$n15xs[i] - a # NOT decimal fraction
ifelse (!is.na(difference) & (difference > 0), CN_mass_data$n15xs_MOM[i] <- difference, CN_mass_data$n15xs_MOM[i] <- 0)
} else {
difference <- CN_mass_data$n_ap[i] - 0.370317701 # assume bkgd of 11 d15N
ifelse ((difference > 0) & !is.na(difference), CN_mass_data$n15xs_MOM[i] <- difference, CN_mass_data$n15xs_MOM[i] <- 0)
}
}
}
}
CN_mass_data$n15_g_pg_recov <- CN_mass_data$n15xs / 100 * CN_mass_data$n_pct # for total recovery
CN_mass_data$n15_g_recov <- CN_mass_data$n15_g_pg_recov * CN_mass_data$g_core # for total recovery
CN_mass_data$n15_g_pg <- CN_mass_data$n15xs_MOM / 100 * CN_mass_data$n_pct
CN_mass_data$n15_g <- CN_mass_data$n15_g_pg * CN_mass_data$g_core
CN_mass_data$n_core <- CN_mass_data$n_pct * CN_mass_data$g_core
## ----15N recovery by pool, include=FALSE, echo=FALSE---------------------
# get 15Nxs mass in each pool in each core
# total recovery using belowground pools (rather than bulk sediment)
mat2 <- plyr::ddply(CN_mass_data[!CN_mass_data$sample.type %in% c("bulk sediment"), ], plyr::.(time, species, new.core.id), plyr::summarise,
n15_2 = sum(n15_g_recov, na.rm = T)
)
mat2$recovery2 <- mat2$n15_2 / (spike) # "recovery" uses bulk belowground data, "recovery2" uses root pools
mat2$days <- as.numeric(substr(mat2$time, 2, 2)) * 7
### relps between napp and 15n recoveries (excludes dead biomass)
bgd <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ],
plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n_core = sum(n_core, na.rm = T)
)
bgd$g_pg <- bgd$n15 / bgd$g
# species-level belowground inventory over time
abv <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("decomp layer", "leaf 1", "leaf 2", "leaf 3", "leaf 4",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
# "dead leaf", "standing dead", "microbe mat",
"standing dead", "stems"), ], plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n_core = sum(n_core, na.rm = T)
)
abv$g_pg <- abv$n15 / abv$g
names(abv)[4:7] <- paste0(names(abv)[4:7], ".ag")
names(bgd)[4:7] <- paste0(names(bgd)[4:7], ".bg")
comb <- plyr::join_all(list(napp, abv, bgd, mat2), by = "new.core.id")
## ----15N inventories: aboveground, include=FALSE, echo=FALSE-------------
###
### aboveground data
###
### mass over time in aboveground compartments. mean +- se by species (sum across depth intervals)
ag2 <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type2 %in% c("stems", "leaf"), ],
plyr::.(time, species, new.core.id, sample.type2), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n15_pg = n15 / (g * pot.m2),
c13 = mean(d13c, na.rm = TRUE),
cn = sum(c_pct * g_core, na.rm = TRUE) / n
)
ag2.sp <- plyr::ddply(ag2, plyr::.(time, species, sample.type2), plyr::numcolwise(mean))
ag2.se <- plyr::ddply(ag2, plyr::.(time, species, sample.type2), plyr::numcolwise(se))
names(ag2.se)[4:10] <- paste0(names(ag2.sp)[4:10], ".se")
ag2.sp <- cbind(ag2.sp, ag2.se[, 4:10])
ag2.sp$t <- as.integer(substr(ag2.sp$time, 2, 2))
ag2.sp$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(ag2.sp$time))]
ag2.sp$species2 <- ifelse(ag2.sp$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
###
### belowground data
###
### mass over time in belowground compartments. mean +- se by species (sum across depth intervals)
bg <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ],
plyr::.(time, species, new.core.id, sample.type), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n15_pg = n15 / (g * pot.m2)
)
bg.sp <- plyr::ddply(bg, plyr::.(time, species, sample.type), plyr::numcolwise(mean))
bg.se <- plyr::ddply(bg, plyr::.(time, species, sample.type), plyr::numcolwise(se))
names(bg.se)[4:8] <- paste0(names(bg.se)[4:8], ".se")
bg.sp <- cbind(bg.sp, bg.se[, 4:8])
bg.sp$t <- as.integer(substr(bg.sp$time, 2, 2))
bg.sp$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(bg.sp$time))]
bg.sp$species2 <- ifelse(bg.sp$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
# Total 15N inventories ---------------------------------------------------------
# aboveground
ag <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems", "leaf", paste0("leaf ", 1:10)), ],
plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T) # g per pot
)
agd1 <- plyr::ddply(ag, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
# species-level belowground inventory over time
bgd <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ], #, "dead biomass"
plyr::.(time, species, new.core.id), plyr::summarise,
n15 = sum(n15_g, na.rm = T)
)
bgd1 <- plyr::ddply(bgd, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
bgd1$type <- "Belowground"
agd1$type <- "Aboveground"
mgd <- rbind(agd1, bgd1)
tot <- plyr::ddply(mat2, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15_2, na.rm = T),
n15.se = se(n15_2),
recovery = mean(n15_2 / spike, na.rm = TRUE),
recovery.se = se(n15_2 / spike)
)
tot$type <- "Total"
# total recovery combined above + belowground (includes dead aboveground but only live belowground)
tot.live <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes",
"leaf 1", "leaf 2", "leaf 3", "leaf 4",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
"stems"), ], #, "dead biomass"
plyr::.(time, species, new.core.id), plyr::summarise,
n15 = sum(n15_g, na.rm = T)
)
tot.live1 <- plyr::ddply(tot.live, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
mgd <- rbind(mgd, tot)
# replace unique times with actual sampling dates using indexing
mgd$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(mgd$time))]
## ----N uptake rates, include=FALSE, echo=FALSE---------------------------
# N uptake rates ----------------------------------------------------------
# Aboveground: estimate from biomass and weighted average biomass N concentration
# get weighted average of aboveground N concentrations for each core
n_mean <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems", "leaf", paste0("leaf ", 1:10)), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
n_bg_mean <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("fine roots", "coarse roots", "rhizomes"), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass_bg = sum(g_core / pot.m2, na.rm = TRUE),
tot_n_bg = sum(n_core / pot.m2, na.rm = TRUE),
tot_c_bg = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa_bg = tot_n_bg / tot_mass_bg,
c_wa_bg = tot_c_bg / tot_mass_bg
)
bg_del <- plyr::ddply(bg, plyr::.(new.core.id), plyr::summarise, # I think all these values are per m2
g_bg = sum(g, na.rm = TRUE),
g_froots = sum(g[sample.type %in% "fine roots"], na.rm = TRUE),
g_roots = sum(g[sample.type %in% c("fine roots", "coarse roots")], na.rm = TRUE),
n_bg = sum(n, na.rm = TRUE),
n_pct_bg = mean(n_pct * g, na.rm = TRUE) / g_bg,
n_pct_bg2 = n_bg / g_bg,
n15_bg = sum(n15, na.rm = TRUE)
)
# all relevant values are per m2
master <- plyr::join_all(list(napp[, c("napp.MH", "new.core.id")], # primary production
rgr.stems,
n_mean[, -c(1)], n_bg_mean[, -c(1:2)], # N inventories above and belowground
ag[, 3:7], bg_del # 15N recoveries and unweighted average pct N
), by = "new.core.id")
master$napp.MH <- master$napp.MH / pot.m2 # g/m2
master$n15_core <- master$n15 + master$n15_bg # total live biomass inventory
master$recovery <- master$n15_core / spike
master$session <- unique(ddHgt4$session)[c(3, 5, 7, 9)][as.numeric(as.factor(master$time))]
# use master$timeDiff for 15N rate calculations
# use master$timeDiff + 2 for NAPP rate calculations
master$timeDiff <- as.numeric(difftime(master$session, "2016-06-24", units = "days"))
master$n_uptake_15n <- master$n15 * 1e3 / master$timeDiff / pot.m2 # mg N/day/m2 accumulating in aboveground tissue, using 15N
master$n_uptake_15n_bg <- master$n15_bg * 1e3 / master$timeDiff /pot.m2 # mg N/day/m2
master$prodn_rate <- master$napp.MH / (master$timeDiff + 2) # g/m2/day
master$n_uptake_biomass <- master$prodn_rate * 1e3 * master$n_wa # mg N/m2/day accumulating in aboveground tissue, using primary production rate and weighted average pct N
master$cn_ag <- master$c_wa / master$n_wa # C:N ratio in aboveground tissue
master$cn_bg <- master$c_wa_bg / master$n_wa_bg
master$cluster <- ifelse(master$time %in% c("t1", "t2"), "1-2 weeks", "3-4 weeks")
master$n15_pgBG <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_bg # 15N uptake per gram belowground biomass: mg 15N/m2/day/g TOTAL BG biomass
master$n15_pgFR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_froots # 15N uptake per gram coarse roots: mg 15N/m2/day/g fine roots
master$n15_pgCR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / (master$g_roots - master$g_froots) # 15N uptake per gram coarse roots: mg 15N/m2/day/g coarse roots
master$n15_pgR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_roots # 15N uptake per gram fine+coarse root biomass: mg 15N/m2/day/g root biomass
master$n15_pgFR2 <- master$n15_core * 1e3 / master$g_froots
master$tot_15n_uptake <- master$n_uptake_15n + master$n_uptake_15n_bg # mg N/day/m2
master$days <- as.numeric(substr(master$time, 2, 2)) * 7
### apply relationship to estimate belowground production
summary(lm3_4 <- stats::lm(n_uptake_biomass ~ I(n_uptake_15n ) , data = master[master$time %in% c("t3", "t4"), ]))
lm3_4
master$bg_n_est <- as.numeric(master$n_uptake_15n_bg) * coefficients(lm3_4)[2] + coefficients(lm3_4)[1] # total N uptake belowground: mg N/m2/day
master$bg_n_est[master$time %in% c("t1", "t2")] <- NA
master$bgp_est <- master$bg_n_est / master$n_pct_bg2 / 1e3 # [estimated total N uptake (mg N/day)] / [mg N / mg biomass] / [1000 mg/g] = [g belowground biomass / day (/m2)]
master$bgp_biomass_est <- master$bgp_est * master$timeDiff # [g belowground biomass / day] * [timeDiff (days)] = total belowground production during period [g/m2]
### differences in total nitrogen interception between species
master$tot_N_int <- master$bg_n_est + master$n_uptake_biomass + ifelse(master$species %in% "SA", nmolHr_mgDay(dd.dea$flux.m2[2]), nmolHr_mgDay(dd.dea$flux.m2[1]))
master$tot_N_int[master$time %in% c("t1", "t2")] <- NA
dd.master <- ddply(master, .(time, species), numcolwise(mean))
dd.master.se <- ddply(master, .(time, species), numcolwise(se))
dd.master$n_uptake_biomass.se <- dd.master.se$n_uptake_biomass
dd.master$prodn_rate.se <- dd.master.se$prodn_rate
dd.master$session <- unique(ddHgt4$session)[c(3, 5, 7, 9)][as.numeric(as.factor(dd.master$time))]
## ----Results section A, include=FALSE, echo=FALSE------------------------
# Results: a. Stem allometry, aboveground biomass and production ----------
table(allometry$spp[!is.na(allometry$sample)])
plyr::ldply(CSP, MuMIn::r.squaredGLMM)
# change over time in production?
summary(lm.sa.prod <- lm(napp.MH ~ session, data = master[master$species %in% "SA",]))
summary(lm.ds.prod <- lm(napp.MH ~ session, data = master[master$species %in% "DS",]))
plyr::ddply(master, plyr::.(species), plyr::summarise,
production = mean(prodn_rate, na.rm = TRUE),
production.se = se(prodn_rate)
)
t.test(prodn_rate ~ species, data = master)
for (i in 1:length(unique(master$session))) {
print(unique(master$session)[i])
print(t.test(prodn_rate ~ species, data = master[(master$session %in% c(unique(master$session)[i])), ]))
}
## ----Results section B, include = FALSE, echo = FALSE--------------------
# Results: b. Leaf and stem biomass, N pools, and uptake ------------------
# leaf, stem biomass differences
leaf <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("leaf", paste0("leaf ", 1:10)), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass)
leaf$t <- as.numeric(substr(leaf$time, 2, 2)) * 7
stem <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems"), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
stem$t <- as.numeric(substr(stem$time, 2, 2)) * 7
leaf2 <- leaf
names(leaf2)[4:8] <- paste0(names(leaf)[4:8], ".leaf")
leaf2 <- cbind(leaf2, stem[4:8])
leaf2$n_wa_ag <- (leaf2$tot_n.leaf + leaf2$tot_n) / (leaf2$tot_mass.leaf + leaf2$tot_mass)
plyr::ddply(leaf[, -c(2, 3)], plyr::.(species), plyr::numcolwise(mean))
plyr::ddply(leaf[, -c(2, 3)], plyr::.(species), plyr::numcolwise(se))
# leaf biomass
summary(lm.sa.leaf <- lm(tot_mass ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leaf <- lm(tot_mass ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(tot_mass ~ species, data = leaf)
# leaf N content
summary(lm.sa.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(n_wa ~ species, data = leaf)
# stems
plyr::ddply(stem[, -c(2, 3)], plyr::.(species), plyr::numcolwise(mean))
plyr::ddply(stem[, -c(2, 3)], plyr::.(species), plyr::numcolwise(se))
# stem biomass
summary(lm.sa.stem <- lm(tot_mass ~ t, data = stem[stem$species %in% "SA",]))
summary(lm.ds.stem <- lm(tot_mass ~ t, data = stem[stem$species %in% "DS",]))
t.test(tot_mass ~ species, data = stem)
# stem N content
summary(lm.sa.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(n_wa ~ species, data = leaf)
# mean N content for all
plyr::ddply(stem[, -c(2, 3)], plyr::.(), plyr::numcolwise(mean))
plyr::ddply(stem[, -c(2, 3)], plyr::.(), plyr::numcolwise(se))
# stem N stock
t.test(tot_n ~ species, data = stem)
# aboveground weighted average N content
summary(lm.sa.N.ag <- lm(n_wa_ag ~ t, data = leaf2[leaf2$species %in% "SA",]))
summary(lm.ds.N.ag <- lm(n_wa_ag ~ t, data = leaf2[leaf2$species %in% "DS",]))
t.test(n_wa_ag ~ species, data = leaf2)
plyr::ddply(leaf2, plyr::.(species, time), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
plyr::ddply(leaf2, plyr::.(species), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
plyr::ddply(leaf2, plyr::.(), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
### total N uptake (NAPP * %N)
summary(lm.sa.Nupt <- lm(n_uptake_biomass ~ session, data = master[master$species %in% "SA",]))
summary(lm.ds.Nupt <- lm(n_uptake_biomass ~ session, data = master[master$species %in% "DS",]))
t.test(n_uptake_biomass ~ species, data = master) # no difference bt species as a whole
for (i in 1:length(unique(master$session))) {
print(unique(master$session)[i])
print(t.test(n_uptake_biomass ~ species, data = master[(master$session %in% c(unique(master$session)[i])), ]))
}
plyr::ddply(master, plyr::.(species), plyr::summarise,
n_uptake = mean(n_uptake_biomass),
n_uptake.se = se(n_uptake_biomass))
## ----Figure 2 (Data In Brief), fig.width = 6, fig.height = 4, message = FALSE, include = FALSE, echo = FALSE----
# Figure 2 in Data In Brief manuscript - aboveground pools ---------------------------------------
ggplot(ag2.sp, aes(y = g, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() %+replace% theme(legend.title = element_blank()) +
facet_grid(sample.type2 ~ .) + geom_errorbar(aes(ymin = g - g.se, ymax = g + g.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab(expression("Biomass (g "%.%m^-2~")")) + xlab("") + ylim(0, 325) +
scale_x_date(breaks = as.Date(unique(ag2.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(2.2, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA))
ggplot(ag2.sp, aes(y = n_pct, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() %+replace% theme(legend.title = element_blank()) +
facet_grid(sample.type2 ~ .) + geom_errorbar(aes(ymin = n_pct - n_pct.se, ymax = n_pct + n_pct.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab("Tissue N") + xlab("") + scale_y_continuous(labels = scales::percent) +
scale_x_date(breaks = as.Date(unique(ag2.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(2.2, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA))
## ----Results: C: Belowground biomass and N pools, include = FALSE, echo = FALSE----
# Results: c. Belowground biomass and N pools -----------------------------
# belowground biomass and tissue N
rts <- ddply(CN_mass_data[CN_mass_data$sample.type %in% c("fine roots", "coarse roots", "rhizomes"), ],
.(species, time, new.core.id, sample.type), summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
rts$t <- as.numeric(substr(rts$time, 2, 2)) * 7
ddply(rts[, -c(2, 3)], .(species, sample.type), numcolwise(mean))
ddply(rts[, -c(2, 3)], .(species, sample.type), numcolwise(se))
# fine roots
organ <- "fine roots"
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
# coarse roots
organ <- "coarse roots"
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
for (i in 1:length(unique(rts$t))) {
print(unique(rts$t)[i])
print(t.test(n_wa ~ species, data = rts[(rts$t %in% c(unique(rts$t)[i])) & (rts$sample.type %in% organ), ]))
}
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
# rhizomes
organ <- "rhizomes"
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
## ----Figure 3 (Data In Brief), fig.width = 6, fig.height = 4, include = FALSE, echo = FALSE----
# Figure 3 (Data In Brief) - belowground biomass pools --------------------------------------
ggplot(bg.sp, aes(y = g, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() +
facet_grid(sample.type ~ .) + geom_errorbar(aes(ymin = g - g.se, ymax = g + g.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab(expression("Biomass (g "%.%m^-2~")")) + xlab("") +
scale_x_date(breaks = as.Date(unique(bg.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(2.2, 1),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
ggplot(bg.sp, aes(y = n_pct, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() +
facet_grid(sample.type ~ .) + geom_errorbar(aes(ymin = n_pct - n_pct.se, ymax = n_pct + n_pct.se), width = 0, position = pd2) +
scale_colour_grey(start = 0.5, end = 0.1, name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(D.~spicata)), expression(italic(S.~alterniflora)))) +
scale_x_date(breaks = as.Date(unique(bg.sp$session)), labels = scales::date_format("%b-%d")) +
ylab("Tissue N content") + xlab("") + scale_y_continuous(labels = scales::percent) +
theme(legend.text.align = 0, legend.position = c(2.2, 1),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Recoveries and sorption, include = FALSE, echo = FALSE--------------
# comparison of total 15N recovery by species
t.test(recovery2 ~ factor(species), data = mat2)
# variation over time in total 15N recovery
Anova(aov.tot <- aov(recovery2 ~ factor(species) * days, data = mat2), type = 2)
summary(lm(recovery2 ~ days, data = mat2[mat2$species %in% "DS",]))
summary(lm(recovery2 ~ days, data = mat2[mat2$species %in% "SA",]))
# comparison of belowground 15N recoveries between species
t.test(n15_bg ~ factor(species), data = master)
Anova(aov.bg <- aov(n15_bg ~ factor(species) * days, data = master), type = 2)
# comparison of aboveground 15N recoveries between species
t.test(n15 ~ factor(species), data = master)
Anova(aov.ag.n15 <- aov(n15 ~ factor(species) * days, data = master), type = 2)
summary(lm(n15 ~ days, data = master[master$species %in% "DS",]))
summary(lm(n15 ~ days, data = master[master$species %in% "SA",]))
t.test(n15 ~ factor(species), data = master)
# Q: how much sorption is there? (look just at belowground data from t1-4)
summary(CN_mass_data$n15xs[!is.na(CN_mass_data$depth_bottom) & (CN_mass_data$time %in% paste0("t", 1:4))])
summary(CN_mass_data$n15xs_MOM[!is.na(CN_mass_data$depth_bottom) & (CN_mass_data$time %in% paste0("t", 1:4))])
# Q: does sorption vary by species? No.
summary(CN_mass_data$n15xs[CN_mass_data$sample.type %in% "dead biomass"])
summary(stats::aov(n15xs ~ species, data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ]))
summary(stats::aov(n15xs ~ species, data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0") & (CN_mass_data$depth_bottom < 11), ])) # doesn't change if focused only on top 10 cm
### does sorption change over time? No
summary(stats::lm(n15xs ~ species + as.numeric(substr(time, 2, 2)), data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ])) # no
summary(stats::lm(n15xs ~ as.numeric(substr(time, 2, 2)), data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ]))
### differences in 15N-nitrogen interception between species (weeks 3 & 4)
# t.test(n_uptake_15n ~ factor(species), data = master[!is.na(master$bgp_biomass_est), ])
# t.test(n_uptake_15n_bg ~ factor(species), data = master[!is.na(master$bgp_biomass_est), ])
ddply(master[!is.na(master$bgp_biomass_est), ], .(species), summarise,
bg.15n.uptake = mean(n_uptake_15n_bg),
bg.15n.uptake.se = se(n_uptake_15n_bg),
ag.15n.uptake = mean(n_uptake_15n),
ag.15n.uptake.se = se(n_uptake_15n))
# 15N uptake per gram fine-roots
master$spp2 <- factor(master$species)
contrasts(master$spp2) <- contr.poly(2)
Anova(aov(n15_pgFR ~ timeDiff * spp2, data = master), type = "III")
summary(lm.sa <- lm(n15_pgFR ~ timeDiff, data = master[master$species %in% "SA", ]))
summary(lm.ds <- lm(n15_pgFR ~ timeDiff, data = master[master$species %in% "DS", ]))
# Calculate intersection of lines of best fit
rightSide <- lm.ds$coefficients[2] - lm.sa$coefficients[2]
leftSide <- lm.sa$coefficients[1] - lm.ds$coefficients[1]
leftSide / rightSide # days until no difference detectable
## ----Drivers of N uptake, include = FALSE, echo = FALSE------------------
summary(lm.mv.sa <- lm(tot_15n_uptake ~ napp.MH + rgr + g_roots + timeDiff, data = master[master$species %in% "SA", ]))
summary(lm.mv.ds <- lm(tot_15n_uptake ~ napp.MH + rgr + g_roots + timeDiff, data = master[master$species %in% "DS", ]))
rsq::rsq.partial(lm.mv.sa)
rsq::rsq.partial(lm.mv.ds)
# belowground production
summary(lm3_4)
### differences in total nitrogen interception between species
t.test(tot_N_int ~ species, data = master)
## ----Table 1, include = FALSE, echo = FALSE------------------------------
### assemble Table 1: 15N recovery and uptake by compartment
# mean(CN_mass_data[(CN_mass_data$sample.type2 %in% "leaf") & (!CN_mass_data$time %in% "t0"), "n15xs_MOM"], na.rm = T)
CN_mass_data$day.no <- as.numeric(substr(CN_mass_data$time, 2, 2)) * 7
tbl.prep <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("coarse roots", "fine roots", "rhizomes",
"stems", "leaf")), ],
.(day.no, species, new.core.id, sample.type2), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep$uptake_rate <- tbl.prep$n15.mg / tbl.prep$day.no
tbl.prep$sp_uptake_rate <- tbl.prep$n15.perg / tbl.prep$day.no
tbl.prep.for.stats <- ddply(tbl.prep, .(species, sample.type2), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.all <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("coarse roots", "fine roots", "rhizomes",
"stems", "leaf")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
tbl.prep.all$uptake_rate <- tbl.prep.all$n15.mg / tbl.prep.all$day.no
tbl.prep.all$sp_uptake_rate <- tbl.prep.all$n15.perg / tbl.prep.all$day.no
tbl.prep.all2 <- ddply(tbl.prep.all, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.all2$sample.type2 <- "Total"
tbl.prep.ag <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("stems", "leaf")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep.ag$uptake_rate <- tbl.prep.ag$n15.mg / tbl.prep.ag$day.no
tbl.prep.ag$sp_uptake_rate <- tbl.prep.ag$n15.perg / tbl.prep.ag$day.no
tbl.prep.ag2 <- ddply(tbl.prep.ag, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.ag2$sample.type2 <- "Total aboveground"
tbl.prep.bg <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("fine roots", "coarse roots", "rhizomes")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep.bg$uptake_rate <- tbl.prep.bg$n15.mg / tbl.prep.bg$day.no
tbl.prep.bg$sp_uptake_rate <- tbl.prep.bg$n15.perg / tbl.prep.bg$day.no
tbl.prep.bg2 <- ddply(tbl.prep.bg, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.bg2$sample.type2 <- "Total belowground"
tbl.final <- plyr::rbind.fill(list(tbl.prep.for.stats, tbl.prep.ag2, tbl.prep.bg2, tbl.prep.all2))
# write.csv(tbl.final, "C:/RDATA/greenhouse/output/GRO/table_15n_180524.csv")
for (i in 1:length(unique(tbl.prep$sample.type2))) {
print(unique(tbl.prep$sample.type2)[i])
print(paste("n15.mg P-value: ", t.test(n15.mg ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
print(paste("del n-15 P-value: ", t.test(n15.xs ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
print(paste("uptake_rate P-value: ", t.test(uptake_rate ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
}
setDat <- tbl.prep.all # tbl.prep.all, tbl.prep.ag, tbl.prep.bg
paste("del n-15 P-value: ", t.test(n15.xs ~ species, data = setDat)$p.value)
paste("n15.mg P-value: ", t.test(n15.mg ~ species, data = setDat)$p.value)
paste("uptake_rate P-value: ", t.test(uptake_rate ~ species, data = setDat)$p.value)
# average recovery across species
ddply(tbl.prep.ag, .(), summarise,
n15 = mean(n15.mg) / (spike * 1e3 / pot.m2),
n15.se = se(n15.mg) / (spike * 1e3 / pot.m2)
)
ddply(tbl.prep.bg, .(), summarise,
n15 = mean(n15.mg) / (spike * 1e3 / pot.m2),
n15.se = se(n15.mg) / (spike * 1e3 / pot.m2)
)
# difference between above and belowground recovery by species
ag$type <- "aboveground"
bgd$type <- "belowground"
bgag <- rbind(ag[(!ag$time %in% "t0"), c(1:3, 7:8)], bgd[(!bgd$time %in% "t0"), ])
Anova(aov2 <- aov(n15 ~ type * species, data = bgag), type = 2)
TukeyHSD(aov2)
## ----Estimating belowground production, include = FALSE, echo = FALSE----
### differences in total nitrogen interception between species
summary(aov.sp1 <- stats::aov(tot_N_int ~ species, data = master))
car::Anova(aov.sp2 <- stats::aov(tot_N_int ~ species, data = master))
stats::TukeyHSD(aov.sp2)
stats::t.test(tot_N_int ~ species, data = master)
plyr::ddply(master, plyr::.(species), plyr::summarise,
Nint = mean(tot_N_int, na.rm = TRUE),
Nint.se = se(tot_N_int))
bgp.out <- plyr::ddply(master[master$time %in% c("t3", "t4"), ], plyr::.(species), plyr::summarise,
bg_n_uptake = mean(bg_n_est), # mg/m2/d
bg_n_uptake.se = se(bg_n_est),
ag_n_uptake = mean(n_uptake_biomass),
ag_n_uptake.se = se(n_uptake_biomass),
bg_production = mean(bgp_est), # g/m2/day
bg_production.se = se(bgp_est),
ag_production = mean(prodn_rate), # g/m2/day
ag_production.se = se(prodn_rate),
new_bg_biomass = mean(bgp_biomass_est), # g/m2
new_bg_biomass.se = se(bgp_biomass_est),
root_shoot_ratio = mean(bgp_est / prodn_rate),
root_shoot_ratio.se = se(bgp_est / prodn_rate)
)
# write.csv(bgp.out, file = "C:/RDATA/greenhouse/output/GRO/bgp_estimates.csv", row.names = FALSE)
# stats
stats::t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DS higher than SA
stats::t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DISP is more productive than SPAL
# stats for table 2
t.test(bg_n_est ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DS higher than SA
t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(bgp_est ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DISP is more productive than SPAL
t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(DEA.m2 ~ species, data = dea)
t.test(IV.m2 ~ species, data = dea)
## ----Results section E, include = FALSE, echo = FALSE, message = FALSE----
# Results e. Denitrification enzyme activity and total N interception ----
t.test(nmolHr_mgDay(DEA.m2) ~ species, data = dea)
t.test(nmolHr_mgDay(IV.m2) ~ species, data = dea)
dea.x <- ddply(dea, .(species), numcolwise(mean))
dea.se <- ddply(dea, .(species), numcolwise(se))
# dea.x
# dea.se
# dea.se[, 6:7] / dea.x[, 6:7]
tableCols <- c("species", "bg_n_est", "n_uptake_biomass", "bgp_est", "prodn_rate")
ddply(master[master$time %in% c("t3", "t4"), tableCols], .(species), numcolwise(mean))
# ddply(master[master$time %in% c("t3", "t4"), tableCols], .(species), numcolwise(se))
# nmolHr_mgDay(dea.x[, c("DEA.m2", "IV.m2")])
# nmolHr_mgDay(dea.se[, c("DEA.m2", "IV.m2")])
## ----Figure 1 - Allometry, fig.width = 6, fig.height = 4, message = FALSE, echo = FALSE----
# png(filename = paste0("C:/RDATA/greenhouse/output/GRO/Figure1_", todaysDate, ".png"), width = 90, height = 90, units = "mm", res = 1000)
# tiff(file = paste0("C:/RDATA/greenhouse/output/GRO/Figure1_", todaysDate, ".png"), width = 90, height = 90, units = "mm", res = 800)
graphics::par(fig = c(0, 1, 0, 1), mar = c(4, 4, 0.5, 0.5))
graphics::plot(sample ~ height_cm, data = allometry[(allometry$spp %in% "DISP"), ],
cex = pointSize / 2, pch = 19, col = fig2Col,
ylab = "Total mass (g)", xlab = "Height (cm)", xlim = c(0, 80), ylim = c(0, 0.65), xaxt = "n",
las = 1, tcl = 0.25, tck = 0.01, bty = "n", yaxs = "i", xaxs = "i")
graphics::abline(h = 0)
graphics::abline(v = 0)
graphics::axis(side = 1, tcl = 0.25, tck = 0.01, at = axTicks(1), labels = graphics::axTicks(1))
graphics::points(sample ~ height_cm, data = allometry[(allometry$spp %in% "SPAL"), ], cex = pointSize / 2, pch = 17)
x <- allometry[(allometry$spp %in% "DISP"), "height_cm"]
x.spal <- allometry[(allometry$spp %in% "SPAL"), "height_cm"]
y.pred2 <- (x * CSP.coef[1, 3] + CSP.coef[1, "(Intercept)"])^(1/CSP.coef$lam[1]) # DISP
y.pred3 <- (x.spal * CSP.coef[2, 3] + CSP.coef[2, "(Intercept)"])^(1/CSP.coef$lam[2]) # SPAL
graphics::lines(x = x[order(y.pred2)], y = y.pred2[order(y.pred2)], lty = 1, lwd = 2, col = fig2Col)
graphics::lines(x = x.spal[order(y.pred3)], y = y.pred3[order(y.pred3)], lty = 1, lwd = 2)
graphics::text(x = 20, y = 0.52, cex = 0.95, expression(italic("Spartina")))
graphics::text(x = 65, y = 0.3, cex = 0.95, expression(italic("Distichlis")))
# dev.off()
## ----Figure 2 (JEMBE), fig.width = 6, fig.height = 4, message = FALSE, echo=FALSE----
# Figure 2 - Aboveground biomass over time ----------------------------------
ggplot2::ggplot(ddHgt4[ddHgt4$cohort > 0, ], ggplot2::aes(y = mass / pot.m2, x = as.Date(session), colour = species, shape = species)) +
ggplot2::geom_point(size = pointSize, position = pd) + ggplot2::theme_classic() +
ggplot2::geom_errorbar(ggplot2::aes(ymin = (mass - mass.se) / pot.m2, ymax = (mass + mass.se) / pot.m2), width = 0, position = pd) +
ggplot2::scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ddHgt4$species[ddHgt4$cohort > 0])[1], unique(ddHgt4$species[ddHgt4$cohort > 0])[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ggplot2::scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ddHgt4$species[ddHgt4$cohort > 0])[1], unique(ddHgt4$species[ddHgt4$cohort > 0])[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ggplot2::scale_x_date(breaks = as.Date(unique(ddHgt4$session))[c(1, 2, 4, 6, 8)], labels = scales::date_format("%b-%d")) +
ggplot2::ylim(0, 600) + ggplot2::facet_grid(. ~ cohort, labeller = label_parsed) +
ggplot2::labs(y = expression("Biomass (g "%.%m^-2~")"), x = "") +
ggplot2::theme(legend.position = c(0.125, 0.87), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Figure 3 (JEMBE), fig.width = 6, fig.height = 4, message = FALSE, echo=FALSE----
# Figure 3A (JEMBE) - NAPP ------------------------------------------
# set dodge width for points & error bars
pd3 <- position_dodge(0.8)
Fig3A <- ggplot(dd.master, aes(y = prodn_rate, x = as.Date(session), shape = species, colour = species)) + geom_point(size = pointSize) + theme_classic() +
# facet_wrap(~ species) +
geom_errorbar(aes(ymin = prodn_rate - prodn_rate.se, ymax = prodn_rate + prodn_rate.se), width = 0) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylim(0, 10.2) +
ylab(expression("Primary production (g"%.%m^{-2}%.%"d"^{-1}~")")) + xlab("") + ggplot2::scale_x_date(breaks = as.Date(unique(dd.master$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(5,5),
legend.background = element_rect(fill = NA, colour = NA)) +
annotate("text", x = as.Date(unique(dd.master$session)[3]), y = 10, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[4]), y = 8, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[1]), y = 8, label = "A")
#Fig3A
# Figure 3B - Aboveground N uptake ----------------------------------------
Fig3B <- ggplot(dd.master, aes(y = n_uptake_biomass, x = as.Date(session), shape = species, colour = species)) + geom_point(size = pointSize, position = pd3) + theme_classic() +
# facet_wrap(~ species) +
geom_errorbar(aes(ymin = n_uptake_biomass - n_uptake_biomass.se, ymax = n_uptake_biomass + n_uptake_biomass.se), width = 0, position = pd3) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylim(0, 150) +
ylab(expression("N uptake (mg N"%.%d^{-1}%.%"m"^{-2}*")")) + xlab("") +
scale_x_date(breaks = as.Date(unique(dd.master$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(5,5),
legend.background = element_rect(fill = NA, colour = NA)) +
annotate("text", x = as.Date(unique(dd.master$session)[3]), y = 147, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[4]), y = 125, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[1]), y = 150, label = "B")
Fig3B
## ----Figures 4-5 (JEMBE), fig.width = 6, fig.height = 4, echo=FALSE------
# Figure 4 - 15N recoveries ------------------------------------------------
ggplot(mgd[!mgd$time %in% "t0", ], aes(x = as.Date(session), y = recovery, colour = species, shape = species)) +
geom_point(size = pointSize, position = pd2) + ylab(expression(" "^15~"N recovery")) + xlab("") +
geom_errorbar(aes(ymin = recovery - recovery.se, ymax = recovery + recovery.se), width = 0, position = pd2) + theme_classic() +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = unique(mgd$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = unique(mgd$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
facet_grid(. ~ type) + scale_y_continuous(labels = scales::percent, lim = c(0, 0.85)) +
scale_x_date(breaks = as.Date(unique(mgd[!mgd$time %in% "t0", "session"])), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(0.15, 0.9), legend.text.align = 0, axis.title.x=element_blank(),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Fig5, fig.width = 6, fig.height = 4, echo=FALSE---------------------
# Figure 5 - 15N uptake per gram fine roots ---------------------------------
rt2 <- plyr::ddply(master, plyr::.(species, timeDiff), plyr::summarise,
xbar = mean(n15_pgFR, na.rm = TRUE),
se = se(n15_pgFR),
session = mean(session))
ggplot(rt2, aes(x = timeDiff, y = xbar, colour = species, shape = species)) + geom_point(size = pointSize) + theme_classic() +
ylim(0, 0.15) +
geom_errorbar(aes(ymin = xbar - se, ymax = xbar + se), width = 0) +
xlab("Experiment duration (days)") +
ylab(expression("N uptake (mg "^{15}~N%.%d^{-1}%.%"gdw"^{-1}~")")) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = unique(rt2$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(rt2$species)), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_x_continuous(breaks = unique(rt2$timeDiff), lim = c(7, 28)) +
theme(legend.text.align = 0, legend.position = c(0.8, 0.8),
legend.background = element_rect(fill = NA, colour = NA)) +
geom_smooth(method = "lm", se = FALSE)
## ----Figure6, fig.width = 6, fig.height = 4, echo=FALSE------------------
# Figure 6 - N uptake estimated by 2 methods --------------------------------
ggplot(master[master$cluster %in% "3-4 weeks", ], aes(x = n_uptake_15n, y = n_uptake_biomass, colour = species, shape = species)) + geom_point(size = pointSize) + theme_classic() +
ylim(0, 150) + # facet_wrap(~ cluster, scales = "free_x") +
ylab(expression("N uptake (mg N "%.%d^{-1}%.%"m"^{-2}*")")) +
xlab(expression(""^15*"N uptake (mg "^15*"N"%.%d^{-1}%.%"m"^{-2}*")")) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(master$species)[1], unique(master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(master$species)[1], unique(master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
theme(legend.text.align = 0, legend.position = c(0.2, 0.9),
legend.background = element_rect(fill = NA, colour = NA)) +
geom_smooth(data = subset(master, cluster %in% "3-4 weeks"), aes(group = 1), method = "lm", se = FALSE, colour = "black") +
geom_text(data = data.frame(cluster = "3-4 weeks", n_uptake_biomass = 50, n_uptake_15n = 3, species = "SA"), label =
"paste(italic(R) ^ 2, \" = 0.45\")", parse = TRUE, colour = "black")
## ----Allometry model table, results = "asis", message = FALSE, echo = FALSE----
knitr::kable(CSP.coef)
## ----Table2, results = "asis", echo = FALSE------------------------------
kable(tbl.final[, 1:6])
## ----Table 3, results = "asis", echo = FALSE-----------------------------
Species <- c("Belowground uptake (mg N)", "Aboveground uptake (mg N)",
"NBPP (g dw)", "NAPP (g dw)",
"DEA (mg N)", "Unamended N2O+N2 flux (mg N)")
Distichlis <- c("62.2 (20.5)", "115.5*** (7.7)", "7.2 (5.9)", "7.5*** (0.7)", "99.4 (32.4)", "39.4 (23.0)")
Spartina <- c("39.6 (7.1)", "49.2*** (9.8)", "5.9 (1.3)", "2.9*** (0.5)", "65.4 (23.6)", "8.9 (1.1)"
)
t3 <- rbind(Species, Distichlis)
t3.2 <- rbind(t3, Spartina)
knitr::kable(t3.2)
Try the NitrogenUptake2016 package in your browser
Any scripts or data that you put into this service are public.
NitrogenUptake2016 documentation built on May 1, 2019, 11:31 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----setup, include=FALSE, echo=FALSE------------------------------------
if(!require(knitr)){
install.packages("knitr", repos='http://cran.us.r-project.org')
}
if(!require(rmarkdown)){
install.packages("rmarkdown", repos='http://cran.us.r-project.org')
}
if(!require(plyr)){
install.packages("plyr", repos='http://cran.us.r-project.org')
}
if(!require(reshape2)){
install.packages("reshape2", repos='http://cran.us.r-project.org')
}
if(!require(ggplot2)){
install.packages("ggplot2", repos='http://cran.us.r-project.org')
}
if(!require(scales)){
install.packages("scales", repos='http://cran.us.r-project.org')
}
if(!require(MuMIn)){
install.packages("MuMIn", repos='http://cran.us.r-project.org')
}
if(!require(rsq)){
install.packages("rsq", repos='http://cran.us.r-project.org')
}
# then load the package:
library(knitr)
library(rmarkdown)
library(plyr)
library(reshape2)
library(ggplot2)
library(scales)
library(MuMIn)
library(rsq)
library(NitrogenUptake2016)
knitr::opts_chunk$set(echo = TRUE, comment=NA)
## ---- echo = FALSE, include=FALSE----------------------------------------
# if the NitrogenUptake2016 package isn't installed, use devtools to do so:
# devtools::install_github("troyhill/NitrogenUptake2016", build_vignettes = TRUE)
# set some constants
todaysDate <- substr(as.character(Sys.time()), 1, 10)
core.area <- pot.m2 <- 0.00801185 # mesocosm surface area (m2)
top.vol <- core.area * 0.05 * 1e6 # cm3 in core: top 5 cm only
pointSize <- 2 # for ggplot graphics
pd <- pd2 <- position_dodge(1.2)
pd3 <- position_dodge(0.8)
grayColor <- 0.55
fig2Col <- "gray55"
## ----allometry, include=FALSE, echo=FALSE--------------------------------
# Allometry ----------------------------------------------------------------
##### Establish relationships following Lu et al. 2016
CSP <- plyr::dlply(allometry, c("spp"), bCM)
CSP.coef <- plyr::ldply(CSP, stats::coef)
CSP.coef$lam <- plyr::ddply(allometry, c("spp"), function(df) bCM(df, lam.only = TRUE))[, "V1"]
CSP.coef$BIC <- plyr::ldply(CSP, stats::BIC)[, "V1"]
CSP.coef$rsq <- plyr::ldply(CSP, MuMIn::r.squaredGLMM)[, "R2m"]
#CSP.coef
## ----applying allometry to estimate plant masses, include=FALSE, echo=FALSE----
### estimate plant masses from allometry
for (i in 1:nrow(stemHeights)) {
if (stemHeights$species[i] %in% "DS") {
stemHeights$mass[i] <- (stemHeights$height_cm[i] * CSP.coef[1, 3] + CSP.coef[1, "(Intercept)"])^(1/CSP.coef$lam[1]) # DISP
} else if (stemHeights$species[i] %in% "SA") {
stemHeights$mass[i] <- (stemHeights$height_cm[i] * CSP.coef[2, 3] + CSP.coef[2, "(Intercept)"])^(1/CSP.coef$lam[2]) # SPAL
}
}
## ----Allometry validation: comparison between predicted and harvested masses, include=FALSE, echo=FALSE----
###
### compare predicted (allometry) and observed masses at harvest
###
# get live + dead biomass for each pot at each date
dat.ld <- plyr::ddply(stemHeights, plyr::.(core_num, species, date), plyr::summarise,
live = sum(mass[dead_live %in% "L"], na.rm = TRUE),
dead = sum(mass[dead_live %in% "D"], na.rm = TRUE)
)
# add cohort number
dat.ld$day <- as.POSIXct(as.character(dat.ld$date), format = "%y%m%d", origin = "1960-01-01")
dat.ld$cohort <- NA
for (i in 1:length(unique(dat.ld$day))) {
uniqueIDS <- unique(dat.ld$core_num[dat.ld$day == sort(unique(dat.ld$day))[i]])
if (length(uniqueIDS) == 3) {
dat.ld$cohort[dat.ld$core_num %in% uniqueIDS] <- ifelse(sum(is.finite(dat.ld$cohort)) == 0, 1, max(dat.ld$cohort, na.rm = TRUE) + 1)
}
}
dat.ld$cohort[is.na(dat.ld$cohort)] <- 0
dat.ld$cohort[dat.ld$cohort == 5] <- 4
dat.ld$year <- 2016
dat.ld$dead <- 0
dat.ld$pot2 <- paste0(dat.ld$core_num, "-", dat.ld$species)
### allometry estimates
for (i in 1:length(unique(dat.ld$pot2))) {
targPot <- unique(dat.ld$pot2)[i]
maxDate <- max(dat.ld$day[dat.ld$pot2 %in% targPot])
tempDat <- dat.ld[(dat.ld$pot2 %in% targPot) & (dat.ld$day == maxDate), ]
if (i == 1) {
dat.ld.sub <- tempDat
}
if (i > 1) {
dat.ld.sub <- rbind(dat.ld.sub, tempDat)
}
}
dat.ld.sub$id <- paste0(dat.ld.sub$species, "-", dat.ld.sub$core_num)
### observed values
obs <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("leaf 1", "leaf 2", "leaf 3", "leaf 4", # "dead leaf",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
"stems"), ], plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T)
)
obs$core_num <- as.character(as.integer(substr(obs$new.core.id, 3, 4)))
obs$core_num[as.integer(obs$core_num) > 12] <- paste0(obs$species[as.integer(obs$core_num) > 12], as.integer(obs$core_num[as.integer(obs$core_num) > 12]) - 12, "_T0")
obs$id <- paste0(obs$species, "-", obs$core_num)
obs <- plyr::join_all(list(obs, dat.ld.sub[, c(4,9:10)]), by = "id")
names(obs)[which(names(obs) %in% "live")] <- "allom.est"
obs$diff <- (obs$allom.est - obs$g) # magnitude accuracy of allometry prediction (g)
obs$diff.pct <- (obs$allom.est - obs$g) / obs$g # percent accuracy of allometry prediction
obs$species2 <- ifelse(obs$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
summary(obs$diff)
## ----Figure 1 (Data In Brief), fig.width = 6, fig.height = 4, include=FALSE, echo=FALSE----
# Figure 1 from Hill et al. (Data In Brief) - predicted vs obs biomass ---------------------------------------
ggplot(obs, aes(x = allom.est / pot.m2, y = g / pot.m2, colour = species2, shape = species2)) +
geom_point(size = pointSize) + theme_classic() %+replace% theme(legend.title = element_blank()) +
labs(x = expression("Predicted biomass (allometry; g"%.%m^-2~")"), y = expression("Measured biomass (g "%.%m^-2~")")) +
ylim(0, 650) + xlim(0, 650) +
scale_colour_grey(start = 0.5, end = 0.1, name = "", breaks = c(unique(obs$species2)[1], unique(obs$species2)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(obs$species2)[1], unique(obs$species2)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
theme(legend.position = c(0.3, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA)) +
geom_abline(slope = 1, intercept = 0, linetype = 2)
# ggsave(paste0("C:/RDATA/greenhouse/output/GRO/FigureS1_", todaysDate, ".eps"), width = 90, height= 90, units = "mm", dpi = 400)
## ----Relative growth rates, include=FALSE, echo=FALSE--------------------
# Relative growth rates ---------------------------------------------------
stemHeights$RGR <- stemHeights$lateStartDate <- NA
for(i in 1:length(unique(stemHeights$id[stemHeights$dead_live %in% "L"]))) { # 830 plants
targPlant <- unique(stemHeights$id[stemHeights$dead_live %in% "L"])[i]
subDat <- stemHeights[stemHeights$id %in% targPlant, ]
if (nrow(subDat) > 1) {
subDat$mass2 <- c(NA, subDat$mass[1:(nrow(subDat)-1)])
subDat$day2 <- c(subDat$day[1], subDat$day[1:(nrow(subDat)-1)])
subDat$RGR <- (log(subDat$mass) - log(subDat$mass2)) / as.numeric(difftime(subDat$day, subDat$day2, units = "days"))
### add RGRs to dat.live object
stemHeights$RGR[(which(rownames(stemHeights) %in% rownames(subDat)))] <- subDat$RGR
if (!is.na(subDat$height[1]) & (subDat$height[1] < 6)) {
stemHeights$lateStartDate[which(rownames(stemHeights) %in% rownames(subDat))] <- 1 # indicates if high-growth period is captured
}
}
}
stemHeights$RGR[!is.finite(stemHeights$RGR)] <- NA
# mean growth rate for each plant, then each mesocosm
mean.rgr <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(id, new.core.id, species), plyr::summarise, meanRGR = mean(RGR, na.rm = TRUE))
mean.rgr.pot <- plyr::ddply(mean.rgr, plyr::.(new.core.id, species), plyr::summarise, RGR = mean(meanRGR, na.rm = TRUE))
mean.rgr.pot <- mean.rgr.pot[1:24, ]
## ----stem density and primary production, include=FALSE, echo=FALSE------
# Stem density and biomass ------------------------------------------------
### summarize by pot and session
dat.live <- stemHeights[stemHeights$dead_live %in% "L", ]
ddHgt2 <- plyr::ddply(dat.live, plyr::.(date, day, core_num, species), plyr::numcolwise(sum, na.rm = TRUE))
ddHgt2$day <- as.POSIXct(ddHgt2$day, origin = "1960-01-01")
ddHgt2$cohort <- NA
for (i in 1:length(unique(ddHgt2$date))) {
uniqueIDS <- unique(ddHgt2$core_num[ddHgt2$date == sort(unique(ddHgt2$date))[i]])
if (length(uniqueIDS) == 3) {
ddHgt2$cohort[ddHgt2$core_num %in% uniqueIDS] <- i
}
}
ddHgt2$cohort <- (ddHgt2$cohort - 1) / 2
ddHgt2$cohort[is.na(ddHgt2$cohort)] <- 0
ddHgt4 <- plyr::ddply(ddHgt2, plyr::.(day, cohort, species), plyr::numcolwise(mean, na.rm = TRUE))
ddHgt4.se <- plyr::ddply(ddHgt2, plyr::.(day, cohort, species), plyr::numcolwise(se))
names(ddHgt4.se) <- paste0(names(ddHgt4.se), ".se")
ddHgt4 <- cbind(ddHgt4, ddHgt4.se[, c(5, 7, 9)])
ddHgt4$session <- as.POSIXct(ddHgt4$day, origin = "1960-01-01")
ddHgt4$species2 <- ifelse(ddHgt4$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
####
#### Calculate Milner-Hughes NAPP
####
test <- nappCalc2(dataset = dat.ld[dat.ld$cohort > 0, ], liveCol = "live", siteCol = "pot2", timeCol = "day", summarize = "TRUE", annualReset = "FALSE")
napp <- test$summary
napp$pot2 <- as.character(napp$site)
napp$core_num <- as.numeric(sapply(strsplit(napp$pot2, "-"), "[[", 1))
napp$species <- sapply(strsplit(napp$pot2, "-"), "[[", 2)
napp$cohort <- dat.ld$cohort[match(napp$pot2, dat.ld$pot2)]
napp$new.core.id <- ifelse(napp$core_num > 9, paste0(napp$species, napp$core_num), paste0(napp$species, "0", napp$core_num))
### build object to merge with master: average RGR, stem characteristics
# mean growth rate for each plant, then each mesocosm
mean.rgr <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(id, new.core.id, species), plyr::summarise, meanRGR = mean(RGR, na.rm = TRUE))
mean.rgr.pot <- plyr::ddply(mean.rgr, plyr::.(new.core.id, species), plyr::summarise, rgr = mean(meanRGR, na.rm = TRUE))
mean.rgr.pot <- mean.rgr.pot[1:24, ]
mean.dens <- plyr::ddply(stemHeights[stemHeights$dead_live %in% "L", ], plyr::.(new.core.id, date), plyr::summarise, stem.density = sum(!is.na(height_cm)))
mean.dens.pot <- plyr::ddply(mean.dens, plyr::.(new.core.id), plyr::summarise, dens = mean(stem.density / pot.m2, na.rm = TRUE)) # mean of all stem density msrmts for a pot
mean.dens.pot <- mean.dens.pot[1:24, ]
rgr.stems <- plyr::join_all(list(mean.rgr.pot, mean.dens.pot), by = "new.core.id")
## ----Denitrification enzyme activity and ambient N2O flux, include = FALSE, echo = FALSE, message = FALSE----
# Denitrification enzyme activity -----------------------------------------
dea$species <- substr(dea$pot, 1, 2)
dea$species2 <- ifelse(dea$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
dea$DEA.m2 <- dea$DEA * dea$bd_gcm3 * (top.vol / core.area) # * 5 * 1e4 # [flux per gram] * [g/cm3] * [5 cm depth] * [cm2/m2]
dea$IV.m2 <- dea$IV * dea$bd_gcm3 * (top.vol / core.area) # * 5 * 1e4
dea$gap.g <- dea$DEA - dea$IV
dea$gap.m2 <- dea$DEA.m2 - dea$IV.m2
dd.dea <- plyr::ddply(dea, plyr::.(species, species2), plyr::summarise,
flux.g = mean(DEA),
flux.m2 = mean(DEA.m2),
flux.g.se = se(DEA),
flux.m2.se = se(DEA.m2)
)
dd.dea$type <- "DEA"
dd.iv <- plyr::ddply(dea, plyr::.(species, species2), plyr::summarise,
flux.g = mean(IV),
flux.m2 = mean(IV.m2),
flux.g.se = se(IV),
flux.m2.se = se(IV.m2)
)
dd.iv$type <- "GHG"
dd.dea <- rbind(dd.dea, dd.iv)
### stats on DEA/IV N2O flux
dea.m <- reshape2::melt(dea, id.vars = c("species", "species2", "pot"), measure.vars = c(2, 3, 8:11))
dea.m$ID <- paste0(dea.m$species, "-", dea.m$variable)
# dea.g.aov <- stats::aov(value ~ factor(variable) * factor(species), data = dea.m[(dea.m$variable %in% c("IV.m2", "DEA.m2")), ])
# car::Anova(dea.g.aov, type = 2, singular = TRUE) # difference between methods, but not species
## ----15N spike, include=FALSE, echo=FALSE--------------------------------
# how much 15N was added to each pot?
KNO3 <- 2.24047 # grams of KNO3 added
soln <- 2027.5 # final volume of spike solution (mls)
KNO3_mw <- 102.1032 # g/mol
N_mw <- 15 # gN/mol KNO3
mls <- 12 # mls added per mesocosm
spike <- KNO3 * (N_mw / KNO3_mw) / soln * mls # total grams of 15N added to each mesocosm
# calculate atom pct from per mil values
CN_mass_data$n_ap <- ap(CN_mass_data$d15n)
## ----15N background, include=FALSE, echo=FALSE---------------------------
# get background 15N AP
bkd <- plyr::ddply(CN_mass_data[as.character(CN_mass_data$time) %in% "t0", ], plyr::.(species, pool_label), plyr::summarise,
n15 = mean(n_ap, na.rm = T),
se.n15 = se(n_ap)
)
## ----15N corrections, include=FALSE, echo=FALSE--------------------------
# subtract background to get 15N excess (if result is negative, 15Nxs = 0)
CN_mass_data$n15xs <- as.numeric(NA)
for(i in 1:nrow(CN_mass_data)) {
if(!is.na(CN_mass_data$species[i])) {
spp <- CN_mass_data$species[i]
pool <- CN_mass_data$pool_label[i]
a <- bkd[(bkd$species %in% spp) & (bkd$pool_label %in% pool), "n15"]
if((is.numeric(a) & (length(a) > 0))) { # excludes all pools we don't have bkd for
difference <- CN_mass_data$n_ap[i] - a # NOT decimal fraction
ifelse (!is.na(difference) & (difference > 0), CN_mass_data$n15xs[i] <- difference, CN_mass_data$n15xs[i] <- 0)
} else {
difference <- CN_mass_data$n_ap[i] - 0.370317701 # assume bkgd of 11 d15N
ifelse ((difference > 0) & !is.na(difference), CN_mass_data$n15xs[i] <- difference, CN_mass_data$n15xs[i] <- 0)
}
}
}
### Additional correction to belowground live biomass, subtracting 15Nxs of dead MOM to account for sorption
CN_mass_data$n15xs_MOM <- CN_mass_data$n15xs # n15xs still useful for total recoveries
for(i in 1:nrow(CN_mass_data)) {
if (CN_mass_data$time[i] %in% paste0("t", 1:4)) {
if(!is.na(CN_mass_data$depth_top[i])) {
spp <- CN_mass_data$species[i]
coreID <- CN_mass_data$new.core.id[i]
pool <- CN_mass_data$pool_label[i]
depth_bottom <- CN_mass_data$depth_bottom[i]
a <- CN_mass_data[(CN_mass_data$new.core.id %in% coreID) & (CN_mass_data$pool_label %in% paste0("dead biomass ", depth_bottom, "cm")) &
(CN_mass_data$depth_bottom == depth_bottom), "n15xs"]
if((is.numeric(a) & (length(a) > 0))) { # excludes all pools we don't have bkd for
difference <- CN_mass_data$n15xs[i] - a # NOT decimal fraction
ifelse (!is.na(difference) & (difference > 0), CN_mass_data$n15xs_MOM[i] <- difference, CN_mass_data$n15xs_MOM[i] <- 0)
} else {
difference <- CN_mass_data$n_ap[i] - 0.370317701 # assume bkgd of 11 d15N
ifelse ((difference > 0) & !is.na(difference), CN_mass_data$n15xs_MOM[i] <- difference, CN_mass_data$n15xs_MOM[i] <- 0)
}
}
}
}
CN_mass_data$n15_g_pg_recov <- CN_mass_data$n15xs / 100 * CN_mass_data$n_pct # for total recovery
CN_mass_data$n15_g_recov <- CN_mass_data$n15_g_pg_recov * CN_mass_data$g_core # for total recovery
CN_mass_data$n15_g_pg <- CN_mass_data$n15xs_MOM / 100 * CN_mass_data$n_pct
CN_mass_data$n15_g <- CN_mass_data$n15_g_pg * CN_mass_data$g_core
CN_mass_data$n_core <- CN_mass_data$n_pct * CN_mass_data$g_core
## ----15N recovery by pool, include=FALSE, echo=FALSE---------------------
# get 15Nxs mass in each pool in each core
# total recovery using belowground pools (rather than bulk sediment)
mat2 <- plyr::ddply(CN_mass_data[!CN_mass_data$sample.type %in% c("bulk sediment"), ], plyr::.(time, species, new.core.id), plyr::summarise,
n15_2 = sum(n15_g_recov, na.rm = T)
)
mat2$recovery2 <- mat2$n15_2 / (spike) # "recovery" uses bulk belowground data, "recovery2" uses root pools
mat2$days <- as.numeric(substr(mat2$time, 2, 2)) * 7
### relps between napp and 15n recoveries (excludes dead biomass)
bgd <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ],
plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n_core = sum(n_core, na.rm = T)
)
bgd$g_pg <- bgd$n15 / bgd$g
# species-level belowground inventory over time
abv <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("decomp layer", "leaf 1", "leaf 2", "leaf 3", "leaf 4",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
# "dead leaf", "standing dead", "microbe mat",
"standing dead", "stems"), ], plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n_core = sum(n_core, na.rm = T)
)
abv$g_pg <- abv$n15 / abv$g
names(abv)[4:7] <- paste0(names(abv)[4:7], ".ag")
names(bgd)[4:7] <- paste0(names(bgd)[4:7], ".bg")
comb <- plyr::join_all(list(napp, abv, bgd, mat2), by = "new.core.id")
## ----15N inventories: aboveground, include=FALSE, echo=FALSE-------------
###
### aboveground data
###
### mass over time in aboveground compartments. mean +- se by species (sum across depth intervals)
ag2 <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type2 %in% c("stems", "leaf"), ],
plyr::.(time, species, new.core.id, sample.type2), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n15_pg = n15 / (g * pot.m2),
c13 = mean(d13c, na.rm = TRUE),
cn = sum(c_pct * g_core, na.rm = TRUE) / n
)
ag2.sp <- plyr::ddply(ag2, plyr::.(time, species, sample.type2), plyr::numcolwise(mean))
ag2.se <- plyr::ddply(ag2, plyr::.(time, species, sample.type2), plyr::numcolwise(se))
names(ag2.se)[4:10] <- paste0(names(ag2.sp)[4:10], ".se")
ag2.sp <- cbind(ag2.sp, ag2.se[, 4:10])
ag2.sp$t <- as.integer(substr(ag2.sp$time, 2, 2))
ag2.sp$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(ag2.sp$time))]
ag2.sp$species2 <- ifelse(ag2.sp$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
###
### belowground data
###
### mass over time in belowground compartments. mean +- se by species (sum across depth intervals)
bg <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ],
plyr::.(time, species, new.core.id, sample.type), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T),
n15_pg = n15 / (g * pot.m2)
)
bg.sp <- plyr::ddply(bg, plyr::.(time, species, sample.type), plyr::numcolwise(mean))
bg.se <- plyr::ddply(bg, plyr::.(time, species, sample.type), plyr::numcolwise(se))
names(bg.se)[4:8] <- paste0(names(bg.se)[4:8], ".se")
bg.sp <- cbind(bg.sp, bg.se[, 4:8])
bg.sp$t <- as.integer(substr(bg.sp$time, 2, 2))
bg.sp$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(bg.sp$time))]
bg.sp$species2 <- ifelse(bg.sp$species %in% "SA", "italic(S.~alterniflora)", "italic(D.~spicata)")
# Total 15N inventories ---------------------------------------------------------
# aboveground
ag <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems", "leaf", paste0("leaf ", 1:10)), ],
plyr::.(time, species, new.core.id), plyr::summarise,
g = sum(g_core / pot.m2, na.rm = T),
n_pct = mean(n_pct, na.rm = TRUE),
n = sum(n_core / pot.m2, na.rm = T),
n15 = sum(n15_g, na.rm = T) # g per pot
)
agd1 <- plyr::ddply(ag, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
# species-level belowground inventory over time
bgd <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes"), ], #, "dead biomass"
plyr::.(time, species, new.core.id), plyr::summarise,
n15 = sum(n15_g, na.rm = T)
)
bgd1 <- plyr::ddply(bgd, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
bgd1$type <- "Belowground"
agd1$type <- "Aboveground"
mgd <- rbind(agd1, bgd1)
tot <- plyr::ddply(mat2, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15_2, na.rm = T),
n15.se = se(n15_2),
recovery = mean(n15_2 / spike, na.rm = TRUE),
recovery.se = se(n15_2 / spike)
)
tot$type <- "Total"
# total recovery combined above + belowground (includes dead aboveground but only live belowground)
tot.live <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("coarse roots", "fine roots", "rhizomes",
"leaf 1", "leaf 2", "leaf 3", "leaf 4",
"leaf 5", "leaf 6", "leaf 7", "leaf 8", "leaf 9", "leaf 10", "leaf 11",
"stems"), ], #, "dead biomass"
plyr::.(time, species, new.core.id), plyr::summarise,
n15 = sum(n15_g, na.rm = T)
)
tot.live1 <- plyr::ddply(tot.live, plyr::.(time, species), plyr::summarise,
n15.mean = mean(n15, na.rm = T),
n15.se = se(n15),
recovery = mean(n15 / spike, na.rm = TRUE),
recovery.se = se(n15 / spike)
)
mgd <- rbind(mgd, tot)
# replace unique times with actual sampling dates using indexing
mgd$session <- unique(ddHgt4$session)[c(1, 3, 5, 7, 9)][as.numeric(as.factor(mgd$time))]
## ----N uptake rates, include=FALSE, echo=FALSE---------------------------
# N uptake rates ----------------------------------------------------------
# Aboveground: estimate from biomass and weighted average biomass N concentration
# get weighted average of aboveground N concentrations for each core
n_mean <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems", "leaf", paste0("leaf ", 1:10)), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
n_bg_mean <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("fine roots", "coarse roots", "rhizomes"), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass_bg = sum(g_core / pot.m2, na.rm = TRUE),
tot_n_bg = sum(n_core / pot.m2, na.rm = TRUE),
tot_c_bg = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa_bg = tot_n_bg / tot_mass_bg,
c_wa_bg = tot_c_bg / tot_mass_bg
)
bg_del <- plyr::ddply(bg, plyr::.(new.core.id), plyr::summarise, # I think all these values are per m2
g_bg = sum(g, na.rm = TRUE),
g_froots = sum(g[sample.type %in% "fine roots"], na.rm = TRUE),
g_roots = sum(g[sample.type %in% c("fine roots", "coarse roots")], na.rm = TRUE),
n_bg = sum(n, na.rm = TRUE),
n_pct_bg = mean(n_pct * g, na.rm = TRUE) / g_bg,
n_pct_bg2 = n_bg / g_bg,
n15_bg = sum(n15, na.rm = TRUE)
)
# all relevant values are per m2
master <- plyr::join_all(list(napp[, c("napp.MH", "new.core.id")], # primary production
rgr.stems,
n_mean[, -c(1)], n_bg_mean[, -c(1:2)], # N inventories above and belowground
ag[, 3:7], bg_del # 15N recoveries and unweighted average pct N
), by = "new.core.id")
master$napp.MH <- master$napp.MH / pot.m2 # g/m2
master$n15_core <- master$n15 + master$n15_bg # total live biomass inventory
master$recovery <- master$n15_core / spike
master$session <- unique(ddHgt4$session)[c(3, 5, 7, 9)][as.numeric(as.factor(master$time))]
# use master$timeDiff for 15N rate calculations
# use master$timeDiff + 2 for NAPP rate calculations
master$timeDiff <- as.numeric(difftime(master$session, "2016-06-24", units = "days"))
master$n_uptake_15n <- master$n15 * 1e3 / master$timeDiff / pot.m2 # mg N/day/m2 accumulating in aboveground tissue, using 15N
master$n_uptake_15n_bg <- master$n15_bg * 1e3 / master$timeDiff /pot.m2 # mg N/day/m2
master$prodn_rate <- master$napp.MH / (master$timeDiff + 2) # g/m2/day
master$n_uptake_biomass <- master$prodn_rate * 1e3 * master$n_wa # mg N/m2/day accumulating in aboveground tissue, using primary production rate and weighted average pct N
master$cn_ag <- master$c_wa / master$n_wa # C:N ratio in aboveground tissue
master$cn_bg <- master$c_wa_bg / master$n_wa_bg
master$cluster <- ifelse(master$time %in% c("t1", "t2"), "1-2 weeks", "3-4 weeks")
master$n15_pgBG <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_bg # 15N uptake per gram belowground biomass: mg 15N/m2/day/g TOTAL BG biomass
master$n15_pgFR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_froots # 15N uptake per gram coarse roots: mg 15N/m2/day/g fine roots
master$n15_pgCR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / (master$g_roots - master$g_froots) # 15N uptake per gram coarse roots: mg 15N/m2/day/g coarse roots
master$n15_pgR <- (master$n_uptake_15n + master$n_uptake_15n_bg) / master$g_roots # 15N uptake per gram fine+coarse root biomass: mg 15N/m2/day/g root biomass
master$n15_pgFR2 <- master$n15_core * 1e3 / master$g_froots
master$tot_15n_uptake <- master$n_uptake_15n + master$n_uptake_15n_bg # mg N/day/m2
master$days <- as.numeric(substr(master$time, 2, 2)) * 7
### apply relationship to estimate belowground production
summary(lm3_4 <- stats::lm(n_uptake_biomass ~ I(n_uptake_15n ) , data = master[master$time %in% c("t3", "t4"), ]))
lm3_4
master$bg_n_est <- as.numeric(master$n_uptake_15n_bg) * coefficients(lm3_4)[2] + coefficients(lm3_4)[1] # total N uptake belowground: mg N/m2/day
master$bg_n_est[master$time %in% c("t1", "t2")] <- NA
master$bgp_est <- master$bg_n_est / master$n_pct_bg2 / 1e3 # [estimated total N uptake (mg N/day)] / [mg N / mg biomass] / [1000 mg/g] = [g belowground biomass / day (/m2)]
master$bgp_biomass_est <- master$bgp_est * master$timeDiff # [g belowground biomass / day] * [timeDiff (days)] = total belowground production during period [g/m2]
### differences in total nitrogen interception between species
master$tot_N_int <- master$bg_n_est + master$n_uptake_biomass + ifelse(master$species %in% "SA", nmolHr_mgDay(dd.dea$flux.m2[2]), nmolHr_mgDay(dd.dea$flux.m2[1]))
master$tot_N_int[master$time %in% c("t1", "t2")] <- NA
dd.master <- ddply(master, .(time, species), numcolwise(mean))
dd.master.se <- ddply(master, .(time, species), numcolwise(se))
dd.master$n_uptake_biomass.se <- dd.master.se$n_uptake_biomass
dd.master$prodn_rate.se <- dd.master.se$prodn_rate
dd.master$session <- unique(ddHgt4$session)[c(3, 5, 7, 9)][as.numeric(as.factor(dd.master$time))]
## ----Results section A, include=FALSE, echo=FALSE------------------------
# Results: a. Stem allometry, aboveground biomass and production ----------
table(allometry$spp[!is.na(allometry$sample)])
plyr::ldply(CSP, MuMIn::r.squaredGLMM)
# change over time in production?
summary(lm.sa.prod <- lm(napp.MH ~ session, data = master[master$species %in% "SA",]))
summary(lm.ds.prod <- lm(napp.MH ~ session, data = master[master$species %in% "DS",]))
plyr::ddply(master, plyr::.(species), plyr::summarise,
production = mean(prodn_rate, na.rm = TRUE),
production.se = se(prodn_rate)
)
t.test(prodn_rate ~ species, data = master)
for (i in 1:length(unique(master$session))) {
print(unique(master$session)[i])
print(t.test(prodn_rate ~ species, data = master[(master$session %in% c(unique(master$session)[i])), ]))
}
## ----Results section B, include = FALSE, echo = FALSE--------------------
# Results: b. Leaf and stem biomass, N pools, and uptake ------------------
# leaf, stem biomass differences
leaf <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("leaf", paste0("leaf ", 1:10)), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass)
leaf$t <- as.numeric(substr(leaf$time, 2, 2)) * 7
stem <- plyr::ddply(CN_mass_data[CN_mass_data$sample.type %in% c("belowground stems", "stems"), ],
plyr::.(species, time, new.core.id), plyr::summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
stem$t <- as.numeric(substr(stem$time, 2, 2)) * 7
leaf2 <- leaf
names(leaf2)[4:8] <- paste0(names(leaf)[4:8], ".leaf")
leaf2 <- cbind(leaf2, stem[4:8])
leaf2$n_wa_ag <- (leaf2$tot_n.leaf + leaf2$tot_n) / (leaf2$tot_mass.leaf + leaf2$tot_mass)
plyr::ddply(leaf[, -c(2, 3)], plyr::.(species), plyr::numcolwise(mean))
plyr::ddply(leaf[, -c(2, 3)], plyr::.(species), plyr::numcolwise(se))
# leaf biomass
summary(lm.sa.leaf <- lm(tot_mass ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leaf <- lm(tot_mass ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(tot_mass ~ species, data = leaf)
# leaf N content
summary(lm.sa.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(n_wa ~ species, data = leaf)
# stems
plyr::ddply(stem[, -c(2, 3)], plyr::.(species), plyr::numcolwise(mean))
plyr::ddply(stem[, -c(2, 3)], plyr::.(species), plyr::numcolwise(se))
# stem biomass
summary(lm.sa.stem <- lm(tot_mass ~ t, data = stem[stem$species %in% "SA",]))
summary(lm.ds.stem <- lm(tot_mass ~ t, data = stem[stem$species %in% "DS",]))
t.test(tot_mass ~ species, data = stem)
# stem N content
summary(lm.sa.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "SA",]))
summary(lm.ds.leafN <- lm(n_wa ~ t, data = leaf[leaf$species %in% "DS",]))
t.test(n_wa ~ species, data = leaf)
# mean N content for all
plyr::ddply(stem[, -c(2, 3)], plyr::.(), plyr::numcolwise(mean))
plyr::ddply(stem[, -c(2, 3)], plyr::.(), plyr::numcolwise(se))
# stem N stock
t.test(tot_n ~ species, data = stem)
# aboveground weighted average N content
summary(lm.sa.N.ag <- lm(n_wa_ag ~ t, data = leaf2[leaf2$species %in% "SA",]))
summary(lm.ds.N.ag <- lm(n_wa_ag ~ t, data = leaf2[leaf2$species %in% "DS",]))
t.test(n_wa_ag ~ species, data = leaf2)
plyr::ddply(leaf2, plyr::.(species, time), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
plyr::ddply(leaf2, plyr::.(species), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
plyr::ddply(leaf2, plyr::.(), plyr::summarise,
n_pct = mean(n_wa_ag),
n_pct.se = se(n_wa_ag))
### total N uptake (NAPP * %N)
summary(lm.sa.Nupt <- lm(n_uptake_biomass ~ session, data = master[master$species %in% "SA",]))
summary(lm.ds.Nupt <- lm(n_uptake_biomass ~ session, data = master[master$species %in% "DS",]))
t.test(n_uptake_biomass ~ species, data = master) # no difference bt species as a whole
for (i in 1:length(unique(master$session))) {
print(unique(master$session)[i])
print(t.test(n_uptake_biomass ~ species, data = master[(master$session %in% c(unique(master$session)[i])), ]))
}
plyr::ddply(master, plyr::.(species), plyr::summarise,
n_uptake = mean(n_uptake_biomass),
n_uptake.se = se(n_uptake_biomass))
## ----Figure 2 (Data In Brief), fig.width = 6, fig.height = 4, message = FALSE, include = FALSE, echo = FALSE----
# Figure 2 in Data In Brief manuscript - aboveground pools ---------------------------------------
ggplot(ag2.sp, aes(y = g, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() %+replace% theme(legend.title = element_blank()) +
facet_grid(sample.type2 ~ .) + geom_errorbar(aes(ymin = g - g.se, ymax = g + g.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab(expression("Biomass (g "%.%m^-2~")")) + xlab("") + ylim(0, 325) +
scale_x_date(breaks = as.Date(unique(ag2.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(2.2, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA))
ggplot(ag2.sp, aes(y = n_pct, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() %+replace% theme(legend.title = element_blank()) +
facet_grid(sample.type2 ~ .) + geom_errorbar(aes(ymin = n_pct - n_pct.se, ymax = n_pct + n_pct.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ag2.sp$species2)[1], unique(ag2.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab("Tissue N") + xlab("") + scale_y_continuous(labels = scales::percent) +
scale_x_date(breaks = as.Date(unique(ag2.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(2.2, 0.9), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA))
## ----Results: C: Belowground biomass and N pools, include = FALSE, echo = FALSE----
# Results: c. Belowground biomass and N pools -----------------------------
# belowground biomass and tissue N
rts <- ddply(CN_mass_data[CN_mass_data$sample.type %in% c("fine roots", "coarse roots", "rhizomes"), ],
.(species, time, new.core.id, sample.type), summarise,
tot_mass = sum(g_core / pot.m2, na.rm = TRUE),
tot_n = sum(n_core / pot.m2, na.rm = TRUE),
tot_c = sum(c_pct * g_core / pot.m2, na.rm = TRUE),
n_wa = tot_n / tot_mass,
c_wa = tot_c / tot_mass
)
rts$t <- as.numeric(substr(rts$time, 2, 2)) * 7
ddply(rts[, -c(2, 3)], .(species, sample.type), numcolwise(mean))
ddply(rts[, -c(2, 3)], .(species, sample.type), numcolwise(se))
# fine roots
organ <- "fine roots"
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
# coarse roots
organ <- "coarse roots"
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
for (i in 1:length(unique(rts$t))) {
print(unique(rts$t)[i])
print(t.test(n_wa ~ species, data = rts[(rts$t %in% c(unique(rts$t)[i])) & (rts$sample.type %in% organ), ]))
}
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
# rhizomes
organ <- "rhizomes"
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "SA"), ]))
summary(lm(tot_mass ~ t, data = rts[(rts$sample.type %in% organ) & (rts$species %in% "DS"), ]))
t.test(tot_mass ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(n_wa ~ species, data = rts[(rts$sample.type %in% organ), ])
t.test(tot_n ~ species, data = rts[(rts$sample.type %in% organ), ])
## ----Figure 3 (Data In Brief), fig.width = 6, fig.height = 4, include = FALSE, echo = FALSE----
# Figure 3 (Data In Brief) - belowground biomass pools --------------------------------------
ggplot(bg.sp, aes(y = g, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() +
facet_grid(sample.type ~ .) + geom_errorbar(aes(ymin = g - g.se, ymax = g + g.se), width = 0, position = pd2) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylab(expression("Biomass (g "%.%m^-2~")")) + xlab("") +
scale_x_date(breaks = as.Date(unique(bg.sp$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(2.2, 1),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
ggplot(bg.sp, aes(y = n_pct, x = as.Date(session), colour = species2, shape = species2)) + geom_point(size = pointSize, position = pd2) + theme_classic() +
facet_grid(sample.type ~ .) + geom_errorbar(aes(ymin = n_pct - n_pct.se, ymax = n_pct + n_pct.se), width = 0, position = pd2) +
scale_colour_grey(start = 0.5, end = 0.1, name = "", breaks = c(unique(bg.sp$species2)[1], unique(bg.sp$species)[2]), labels = c(expression(italic(D.~spicata)), expression(italic(S.~alterniflora)))) +
scale_x_date(breaks = as.Date(unique(bg.sp$session)), labels = scales::date_format("%b-%d")) +
ylab("Tissue N content") + xlab("") + scale_y_continuous(labels = scales::percent) +
theme(legend.text.align = 0, legend.position = c(2.2, 1),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Recoveries and sorption, include = FALSE, echo = FALSE--------------
# comparison of total 15N recovery by species
t.test(recovery2 ~ factor(species), data = mat2)
# variation over time in total 15N recovery
Anova(aov.tot <- aov(recovery2 ~ factor(species) * days, data = mat2), type = 2)
summary(lm(recovery2 ~ days, data = mat2[mat2$species %in% "DS",]))
summary(lm(recovery2 ~ days, data = mat2[mat2$species %in% "SA",]))
# comparison of belowground 15N recoveries between species
t.test(n15_bg ~ factor(species), data = master)
Anova(aov.bg <- aov(n15_bg ~ factor(species) * days, data = master), type = 2)
# comparison of aboveground 15N recoveries between species
t.test(n15 ~ factor(species), data = master)
Anova(aov.ag.n15 <- aov(n15 ~ factor(species) * days, data = master), type = 2)
summary(lm(n15 ~ days, data = master[master$species %in% "DS",]))
summary(lm(n15 ~ days, data = master[master$species %in% "SA",]))
t.test(n15 ~ factor(species), data = master)
# Q: how much sorption is there? (look just at belowground data from t1-4)
summary(CN_mass_data$n15xs[!is.na(CN_mass_data$depth_bottom) & (CN_mass_data$time %in% paste0("t", 1:4))])
summary(CN_mass_data$n15xs_MOM[!is.na(CN_mass_data$depth_bottom) & (CN_mass_data$time %in% paste0("t", 1:4))])
# Q: does sorption vary by species? No.
summary(CN_mass_data$n15xs[CN_mass_data$sample.type %in% "dead biomass"])
summary(stats::aov(n15xs ~ species, data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ]))
summary(stats::aov(n15xs ~ species, data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0") & (CN_mass_data$depth_bottom < 11), ])) # doesn't change if focused only on top 10 cm
### does sorption change over time? No
summary(stats::lm(n15xs ~ species + as.numeric(substr(time, 2, 2)), data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ])) # no
summary(stats::lm(n15xs ~ as.numeric(substr(time, 2, 2)), data = CN_mass_data[(CN_mass_data$sample.type %in% "dead biomass") & (!CN_mass_data$time %in% "t0"), ]))
### differences in 15N-nitrogen interception between species (weeks 3 & 4)
# t.test(n_uptake_15n ~ factor(species), data = master[!is.na(master$bgp_biomass_est), ])
# t.test(n_uptake_15n_bg ~ factor(species), data = master[!is.na(master$bgp_biomass_est), ])
ddply(master[!is.na(master$bgp_biomass_est), ], .(species), summarise,
bg.15n.uptake = mean(n_uptake_15n_bg),
bg.15n.uptake.se = se(n_uptake_15n_bg),
ag.15n.uptake = mean(n_uptake_15n),
ag.15n.uptake.se = se(n_uptake_15n))
# 15N uptake per gram fine-roots
master$spp2 <- factor(master$species)
contrasts(master$spp2) <- contr.poly(2)
Anova(aov(n15_pgFR ~ timeDiff * spp2, data = master), type = "III")
summary(lm.sa <- lm(n15_pgFR ~ timeDiff, data = master[master$species %in% "SA", ]))
summary(lm.ds <- lm(n15_pgFR ~ timeDiff, data = master[master$species %in% "DS", ]))
# Calculate intersection of lines of best fit
rightSide <- lm.ds$coefficients[2] - lm.sa$coefficients[2]
leftSide <- lm.sa$coefficients[1] - lm.ds$coefficients[1]
leftSide / rightSide # days until no difference detectable
## ----Drivers of N uptake, include = FALSE, echo = FALSE------------------
summary(lm.mv.sa <- lm(tot_15n_uptake ~ napp.MH + rgr + g_roots + timeDiff, data = master[master$species %in% "SA", ]))
summary(lm.mv.ds <- lm(tot_15n_uptake ~ napp.MH + rgr + g_roots + timeDiff, data = master[master$species %in% "DS", ]))
rsq::rsq.partial(lm.mv.sa)
rsq::rsq.partial(lm.mv.ds)
# belowground production
summary(lm3_4)
### differences in total nitrogen interception between species
t.test(tot_N_int ~ species, data = master)
## ----Table 1, include = FALSE, echo = FALSE------------------------------
### assemble Table 1: 15N recovery and uptake by compartment
# mean(CN_mass_data[(CN_mass_data$sample.type2 %in% "leaf") & (!CN_mass_data$time %in% "t0"), "n15xs_MOM"], na.rm = T)
CN_mass_data$day.no <- as.numeric(substr(CN_mass_data$time, 2, 2)) * 7
tbl.prep <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("coarse roots", "fine roots", "rhizomes",
"stems", "leaf")), ],
.(day.no, species, new.core.id, sample.type2), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep$uptake_rate <- tbl.prep$n15.mg / tbl.prep$day.no
tbl.prep$sp_uptake_rate <- tbl.prep$n15.perg / tbl.prep$day.no
tbl.prep.for.stats <- ddply(tbl.prep, .(species, sample.type2), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.all <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("coarse roots", "fine roots", "rhizomes",
"stems", "leaf")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
tbl.prep.all$uptake_rate <- tbl.prep.all$n15.mg / tbl.prep.all$day.no
tbl.prep.all$sp_uptake_rate <- tbl.prep.all$n15.perg / tbl.prep.all$day.no
tbl.prep.all2 <- ddply(tbl.prep.all, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.all2$sample.type2 <- "Total"
tbl.prep.ag <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("stems", "leaf")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep.ag$uptake_rate <- tbl.prep.ag$n15.mg / tbl.prep.ag$day.no
tbl.prep.ag$sp_uptake_rate <- tbl.prep.ag$n15.perg / tbl.prep.ag$day.no
tbl.prep.ag2 <- ddply(tbl.prep.ag, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.ag2$sample.type2 <- "Total aboveground"
tbl.prep.bg <- ddply(CN_mass_data[(!CN_mass_data$time %in% "t0") &
(CN_mass_data$sample.type2 %in% c("fine roots", "coarse roots", "rhizomes")), ],
.(day.no, species, new.core.id), summarise,
n15.xs = sum(n15xs_MOM * n_core, na.rm = T) / sum(n_core, na.rm = TRUE),
g = sum(g_core / pot.m2, na.rm = T),
n15.mg = sum(n15_g, na.rm = T) * 1000 / pot.m2, # mg per m2
n15.perg = n15.mg * pot.m2 / sum(g_core, na.rm = T)) # mg 15N per g dw
# uptake_rate = mean(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / day.no,
# sp_uptake_rate = sum(ap(d15n) * n_pct * g_core, na.rm = T) * 1000 / sum(g_core, na.rm = T) / day.no
tbl.prep.bg$uptake_rate <- tbl.prep.bg$n15.mg / tbl.prep.bg$day.no
tbl.prep.bg$sp_uptake_rate <- tbl.prep.bg$n15.perg / tbl.prep.bg$day.no
tbl.prep.bg2 <- ddply(tbl.prep.bg, .(species), summarise,
n15xs = mean(n15.xs, na.rm = TRUE),
n15xs.se = se(n15.xs),
n15 = mean(n15.mg),
n15.se = se(n15.mg),
n15.pg = mean(n15.perg),
n15.pg.se = se(n15.perg),
uptake = mean(uptake_rate),
uptake.se = se(uptake_rate),
uptake.pg = mean(sp_uptake_rate),
uptake.pg.se = se(sp_uptake_rate)
)
tbl.prep.bg2$sample.type2 <- "Total belowground"
tbl.final <- plyr::rbind.fill(list(tbl.prep.for.stats, tbl.prep.ag2, tbl.prep.bg2, tbl.prep.all2))
# write.csv(tbl.final, "C:/RDATA/greenhouse/output/GRO/table_15n_180524.csv")
for (i in 1:length(unique(tbl.prep$sample.type2))) {
print(unique(tbl.prep$sample.type2)[i])
print(paste("n15.mg P-value: ", t.test(n15.mg ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
print(paste("del n-15 P-value: ", t.test(n15.xs ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
print(paste("uptake_rate P-value: ", t.test(uptake_rate ~ species, data = tbl.prep[tbl.prep$sample.type2 %in% unique(tbl.prep$sample.type2)[i], ])$p.value))
}
setDat <- tbl.prep.all # tbl.prep.all, tbl.prep.ag, tbl.prep.bg
paste("del n-15 P-value: ", t.test(n15.xs ~ species, data = setDat)$p.value)
paste("n15.mg P-value: ", t.test(n15.mg ~ species, data = setDat)$p.value)
paste("uptake_rate P-value: ", t.test(uptake_rate ~ species, data = setDat)$p.value)
# average recovery across species
ddply(tbl.prep.ag, .(), summarise,
n15 = mean(n15.mg) / (spike * 1e3 / pot.m2),
n15.se = se(n15.mg) / (spike * 1e3 / pot.m2)
)
ddply(tbl.prep.bg, .(), summarise,
n15 = mean(n15.mg) / (spike * 1e3 / pot.m2),
n15.se = se(n15.mg) / (spike * 1e3 / pot.m2)
)
# difference between above and belowground recovery by species
ag$type <- "aboveground"
bgd$type <- "belowground"
bgag <- rbind(ag[(!ag$time %in% "t0"), c(1:3, 7:8)], bgd[(!bgd$time %in% "t0"), ])
Anova(aov2 <- aov(n15 ~ type * species, data = bgag), type = 2)
TukeyHSD(aov2)
## ----Estimating belowground production, include = FALSE, echo = FALSE----
### differences in total nitrogen interception between species
summary(aov.sp1 <- stats::aov(tot_N_int ~ species, data = master))
car::Anova(aov.sp2 <- stats::aov(tot_N_int ~ species, data = master))
stats::TukeyHSD(aov.sp2)
stats::t.test(tot_N_int ~ species, data = master)
plyr::ddply(master, plyr::.(species), plyr::summarise,
Nint = mean(tot_N_int, na.rm = TRUE),
Nint.se = se(tot_N_int))
bgp.out <- plyr::ddply(master[master$time %in% c("t3", "t4"), ], plyr::.(species), plyr::summarise,
bg_n_uptake = mean(bg_n_est), # mg/m2/d
bg_n_uptake.se = se(bg_n_est),
ag_n_uptake = mean(n_uptake_biomass),
ag_n_uptake.se = se(n_uptake_biomass),
bg_production = mean(bgp_est), # g/m2/day
bg_production.se = se(bgp_est),
ag_production = mean(prodn_rate), # g/m2/day
ag_production.se = se(prodn_rate),
new_bg_biomass = mean(bgp_biomass_est), # g/m2
new_bg_biomass.se = se(bgp_biomass_est),
root_shoot_ratio = mean(bgp_est / prodn_rate),
root_shoot_ratio.se = se(bgp_est / prodn_rate)
)
# write.csv(bgp.out, file = "C:/RDATA/greenhouse/output/GRO/bgp_estimates.csv", row.names = FALSE)
# stats
stats::t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DS higher than SA
stats::t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DISP is more productive than SPAL
# stats for table 2
t.test(bg_n_est ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DS higher than SA
t.test(n_uptake_biomass ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(bgp_est ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ]) # DISP is more productive than SPAL
t.test(prodn_rate ~ species, data = master[master$time %in% c("t3", "t4"), ])
t.test(DEA.m2 ~ species, data = dea)
t.test(IV.m2 ~ species, data = dea)
## ----Results section E, include = FALSE, echo = FALSE, message = FALSE----
# Results e. Denitrification enzyme activity and total N interception ----
t.test(nmolHr_mgDay(DEA.m2) ~ species, data = dea)
t.test(nmolHr_mgDay(IV.m2) ~ species, data = dea)
dea.x <- ddply(dea, .(species), numcolwise(mean))
dea.se <- ddply(dea, .(species), numcolwise(se))
# dea.x
# dea.se
# dea.se[, 6:7] / dea.x[, 6:7]
tableCols <- c("species", "bg_n_est", "n_uptake_biomass", "bgp_est", "prodn_rate")
ddply(master[master$time %in% c("t3", "t4"), tableCols], .(species), numcolwise(mean))
# ddply(master[master$time %in% c("t3", "t4"), tableCols], .(species), numcolwise(se))
# nmolHr_mgDay(dea.x[, c("DEA.m2", "IV.m2")])
# nmolHr_mgDay(dea.se[, c("DEA.m2", "IV.m2")])
## ----Figure 1 - Allometry, fig.width = 6, fig.height = 4, message = FALSE, echo = FALSE----
# png(filename = paste0("C:/RDATA/greenhouse/output/GRO/Figure1_", todaysDate, ".png"), width = 90, height = 90, units = "mm", res = 1000)
# tiff(file = paste0("C:/RDATA/greenhouse/output/GRO/Figure1_", todaysDate, ".png"), width = 90, height = 90, units = "mm", res = 800)
graphics::par(fig = c(0, 1, 0, 1), mar = c(4, 4, 0.5, 0.5))
graphics::plot(sample ~ height_cm, data = allometry[(allometry$spp %in% "DISP"), ],
cex = pointSize / 2, pch = 19, col = fig2Col,
ylab = "Total mass (g)", xlab = "Height (cm)", xlim = c(0, 80), ylim = c(0, 0.65), xaxt = "n",
las = 1, tcl = 0.25, tck = 0.01, bty = "n", yaxs = "i", xaxs = "i")
graphics::abline(h = 0)
graphics::abline(v = 0)
graphics::axis(side = 1, tcl = 0.25, tck = 0.01, at = axTicks(1), labels = graphics::axTicks(1))
graphics::points(sample ~ height_cm, data = allometry[(allometry$spp %in% "SPAL"), ], cex = pointSize / 2, pch = 17)
x <- allometry[(allometry$spp %in% "DISP"), "height_cm"]
x.spal <- allometry[(allometry$spp %in% "SPAL"), "height_cm"]
y.pred2 <- (x * CSP.coef[1, 3] + CSP.coef[1, "(Intercept)"])^(1/CSP.coef$lam[1]) # DISP
y.pred3 <- (x.spal * CSP.coef[2, 3] + CSP.coef[2, "(Intercept)"])^(1/CSP.coef$lam[2]) # SPAL
graphics::lines(x = x[order(y.pred2)], y = y.pred2[order(y.pred2)], lty = 1, lwd = 2, col = fig2Col)
graphics::lines(x = x.spal[order(y.pred3)], y = y.pred3[order(y.pred3)], lty = 1, lwd = 2)
graphics::text(x = 20, y = 0.52, cex = 0.95, expression(italic("Spartina")))
graphics::text(x = 65, y = 0.3, cex = 0.95, expression(italic("Distichlis")))
# dev.off()
## ----Figure 2 (JEMBE), fig.width = 6, fig.height = 4, message = FALSE, echo=FALSE----
# Figure 2 - Aboveground biomass over time ----------------------------------
ggplot2::ggplot(ddHgt4[ddHgt4$cohort > 0, ], ggplot2::aes(y = mass / pot.m2, x = as.Date(session), colour = species, shape = species)) +
ggplot2::geom_point(size = pointSize, position = pd) + ggplot2::theme_classic() +
ggplot2::geom_errorbar(ggplot2::aes(ymin = (mass - mass.se) / pot.m2, ymax = (mass + mass.se) / pot.m2), width = 0, position = pd) +
ggplot2::scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(ddHgt4$species[ddHgt4$cohort > 0])[1], unique(ddHgt4$species[ddHgt4$cohort > 0])[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ggplot2::scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(ddHgt4$species[ddHgt4$cohort > 0])[1], unique(ddHgt4$species[ddHgt4$cohort > 0])[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ggplot2::scale_x_date(breaks = as.Date(unique(ddHgt4$session))[c(1, 2, 4, 6, 8)], labels = scales::date_format("%b-%d")) +
ggplot2::ylim(0, 600) + ggplot2::facet_grid(. ~ cohort, labeller = label_parsed) +
ggplot2::labs(y = expression("Biomass (g "%.%m^-2~")"), x = "") +
ggplot2::theme(legend.position = c(0.125, 0.87), legend.text.align = 0,
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Figure 3 (JEMBE), fig.width = 6, fig.height = 4, message = FALSE, echo=FALSE----
# Figure 3A (JEMBE) - NAPP ------------------------------------------
# set dodge width for points & error bars
pd3 <- position_dodge(0.8)
Fig3A <- ggplot(dd.master, aes(y = prodn_rate, x = as.Date(session), shape = species, colour = species)) + geom_point(size = pointSize) + theme_classic() +
# facet_wrap(~ species) +
geom_errorbar(aes(ymin = prodn_rate - prodn_rate.se, ymax = prodn_rate + prodn_rate.se), width = 0) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylim(0, 10.2) +
ylab(expression("Primary production (g"%.%m^{-2}%.%"d"^{-1}~")")) + xlab("") + ggplot2::scale_x_date(breaks = as.Date(unique(dd.master$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(5,5),
legend.background = element_rect(fill = NA, colour = NA)) +
annotate("text", x = as.Date(unique(dd.master$session)[3]), y = 10, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[4]), y = 8, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[1]), y = 8, label = "A")
#Fig3A
# Figure 3B - Aboveground N uptake ----------------------------------------
Fig3B <- ggplot(dd.master, aes(y = n_uptake_biomass, x = as.Date(session), shape = species, colour = species)) + geom_point(size = pointSize, position = pd3) + theme_classic() +
# facet_wrap(~ species) +
geom_errorbar(aes(ymin = n_uptake_biomass - n_uptake_biomass.se, ymax = n_uptake_biomass + n_uptake_biomass.se), width = 0, position = pd3) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(dd.master$species)[1], unique(dd.master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
ylim(0, 150) +
ylab(expression("N uptake (mg N"%.%d^{-1}%.%"m"^{-2}*")")) + xlab("") +
scale_x_date(breaks = as.Date(unique(dd.master$session)), labels = scales::date_format("%b-%d")) +
theme(legend.text.align = 0, legend.position = c(5,5),
legend.background = element_rect(fill = NA, colour = NA)) +
annotate("text", x = as.Date(unique(dd.master$session)[3]), y = 147, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[4]), y = 125, label = "*") +
annotate("text", x = as.Date(unique(dd.master$session)[1]), y = 150, label = "B")
Fig3B
## ----Figures 4-5 (JEMBE), fig.width = 6, fig.height = 4, echo=FALSE------
# Figure 4 - 15N recoveries ------------------------------------------------
ggplot(mgd[!mgd$time %in% "t0", ], aes(x = as.Date(session), y = recovery, colour = species, shape = species)) +
geom_point(size = pointSize, position = pd2) + ylab(expression(" "^15~"N recovery")) + xlab("") +
geom_errorbar(aes(ymin = recovery - recovery.se, ymax = recovery + recovery.se), width = 0, position = pd2) + theme_classic() +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = unique(mgd$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = unique(mgd$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
facet_grid(. ~ type) + scale_y_continuous(labels = scales::percent, lim = c(0, 0.85)) +
scale_x_date(breaks = as.Date(unique(mgd[!mgd$time %in% "t0", "session"])), labels = scales::date_format("%b-%d")) +
theme(legend.position = c(0.15, 0.9), legend.text.align = 0, axis.title.x=element_blank(),
legend.background = element_rect(fill = NA, colour = NA), axis.text.x=element_text(angle=45,hjust=1))
## ----Fig5, fig.width = 6, fig.height = 4, echo=FALSE---------------------
# Figure 5 - 15N uptake per gram fine roots ---------------------------------
rt2 <- plyr::ddply(master, plyr::.(species, timeDiff), plyr::summarise,
xbar = mean(n15_pgFR, na.rm = TRUE),
se = se(n15_pgFR),
session = mean(session))
ggplot(rt2, aes(x = timeDiff, y = xbar, colour = species, shape = species)) + geom_point(size = pointSize) + theme_classic() +
ylim(0, 0.15) +
geom_errorbar(aes(ymin = xbar - se, ymax = xbar + se), width = 0) +
xlab("Experiment duration (days)") +
ylab(expression("N uptake (mg "^{15}~N%.%d^{-1}%.%"gdw"^{-1}~")")) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = unique(rt2$species), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(rt2$species)), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_x_continuous(breaks = unique(rt2$timeDiff), lim = c(7, 28)) +
theme(legend.text.align = 0, legend.position = c(0.8, 0.8),
legend.background = element_rect(fill = NA, colour = NA)) +
geom_smooth(method = "lm", se = FALSE)
## ----Figure6, fig.width = 6, fig.height = 4, echo=FALSE------------------
# Figure 6 - N uptake estimated by 2 methods --------------------------------
ggplot(master[master$cluster %in% "3-4 weeks", ], aes(x = n_uptake_15n, y = n_uptake_biomass, colour = species, shape = species)) + geom_point(size = pointSize) + theme_classic() +
ylim(0, 150) + # facet_wrap(~ cluster, scales = "free_x") +
ylab(expression("N uptake (mg N "%.%d^{-1}%.%"m"^{-2}*")")) +
xlab(expression(""^15*"N uptake (mg "^15*"N"%.%d^{-1}%.%"m"^{-2}*")")) +
scale_colour_grey(start = grayColor, end = 0.1, name = "", breaks = c(unique(master$species)[1], unique(master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
scale_shape_manual(values = c(16, 17), name = "", breaks = c(unique(master$species)[1], unique(master$species)[2]), labels = c(expression(italic(Distichlis)), expression(italic(Spartina)))) +
theme(legend.text.align = 0, legend.position = c(0.2, 0.9),
legend.background = element_rect(fill = NA, colour = NA)) +
geom_smooth(data = subset(master, cluster %in% "3-4 weeks"), aes(group = 1), method = "lm", se = FALSE, colour = "black") +
geom_text(data = data.frame(cluster = "3-4 weeks", n_uptake_biomass = 50, n_uptake_15n = 3, species = "SA"), label =
"paste(italic(R) ^ 2, \" = 0.45\")", parse = TRUE, colour = "black")
## ----Allometry model table, results = "asis", message = FALSE, echo = FALSE----
knitr::kable(CSP.coef)
## ----Table2, results = "asis", echo = FALSE------------------------------
kable(tbl.final[, 1:6])
## ----Table 3, results = "asis", echo = FALSE-----------------------------
Species <- c("Belowground uptake (mg N)", "Aboveground uptake (mg N)",
"NBPP (g dw)", "NAPP (g dw)",
"DEA (mg N)", "Unamended N2O+N2 flux (mg N)")
Distichlis <- c("62.2 (20.5)", "115.5*** (7.7)", "7.2 (5.9)", "7.5*** (0.7)", "99.4 (32.4)", "39.4 (23.0)")
Spartina <- c("39.6 (7.1)", "49.2*** (9.8)", "5.9 (1.3)", "2.9*** (0.5)", "65.4 (23.6)", "8.9 (1.1)"
)
t3 <- rbind(Species, Distichlis)
t3.2 <- rbind(t3, Spartina)
knitr::kable(t3.2)
Try the NitrogenUptake2016 package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.