In troyhill/LUMCONplants: Allometry and Plant Dynamics in Louisiana Universities Marine Consortium Salt Marshes
### This section loads libraries, sets some local variables
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
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(gridExtra)){
install.packages("gridExtra", 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')
}
if(!require(multcomp)){
install.packages("multcomp", repos='http://cran.us.r-project.org')
}
if(!require(sandwich)){
install.packages("sandwich", repos='http://cran.us.r-project.org')
}
if(!require(car)){
install.packages("car", repos='http://cran.us.r-project.org')
}
if(!require(Rmisc)){
install.packages("Rmisc", repos='http://cran.us.r-project.org')
}
if(!require(tidyverse)){
install.packages("tidyverse", repos='http://cran.us.r-project.org')
}
if(!require(lattice)){
install.packages("lattice", repos='http://cran.us.r-project.org')
}
if(!require(nlme)){
install.packages("nlme", repos='http://cran.us.r-project.org')
}
if(!require(lubridate)){
install.packages("lubridate", repos='http://cran.us.r-project.org')
}
# if(!require(svglite)){
# install.packages("svglite", repos='http://cran.us.r-project.org')
# }
library(ggplot2)
library(gridExtra)
library(plyr)
library(LUMCONplants)
library(multcomp) # for glht
library(sandwich) # for HC3 estimator
library(car) # for Levene's test
library(Rmisc)
library(tidyverse)
library(lattice)
library(nlme)
library(lubridate)
# library(svglite) # doesn't install properly on linux
### define some constants
plotSize <- 0.25 * 0.25 # size of quadrat, m2
coreTube <- pi*(6.9/2)^2 # area of plastic coring tube used prior to August 2016: units = cm2
coreTube_metal <- pi*(5.8/2)^2 # area of metal coring device used beginning in August 2016: units = cm2
growingMonths <- 4:10 # growing season defined as April:October
bagMass <- 5.4 # grams; bags used for belowground cores
tckSize <- 0.01
getDemographics <- function(data = stemDat, status = "LIVE", groupBreaks = c(0, 25, 50, 100, 200), groupLabels = c("0-25 cm", "25-50 cm", "50-100 cm", ">100 cm"), quadrat_m2 = 1) {
df1 <- transform(data[(data$status %in% status), ], group = cut(hgt,
breaks = groupBreaks,
labels = groupLabels))
res <- do.call(data.frame, aggregate(hgt ~ time + group, df1,
FUN = function(x) c(count = sum(!is.na(x)) / quadrat_m2, mean = mean(x, na.rm = TRUE) / quadrat_m2)))
dNew <- data.frame(group = levels(df1$group))
df2 <- merge(res, dNew, all = TRUE)
# add missing bins, set to 0 hgt/dens
missedLevs <- expand.grid(time = unique(res$time), group = unique(res$group), hgt.count = 0, hgt.mean = 0)
missingLevs <- missedLevs[which(!paste0(missedLevs$time, missedLevs$group) %in% paste0(res$time, res$group)), ]
df2 <- rbind(df2, missingLevs)
#df2$site <- targReg
invisible(df2)
}
regionLabel <- function(data, siteCol = "site") {
# function takes a dataset and the name of the column with site names (e.g., "LUM1"), and
# adds a column with marsh names ("LUM", "TB-A", "TB-B")
data$region <- as.character(NA)
for (i in 1:nrow(data)) {
if (data[i, siteCol] %in% paste0("LUM", 1:3)) {
data$region[i] <- "LUMCON"
} else if (data[i, siteCol] %in% c("TB1", "TB2")) {
data$region[i] <- "Bay LaFleur"
} else if (data[i, siteCol] %in% c("TB3", "TB4")) {
data$region[i] <- "Lake Barre"
}
}
data
}
### a function to make correlation matrices with significance stars
### from http://myowelt.blogspot.com/2008/04/beautiful-correlation-tables-in-r.html
### usage: kable(corstarsl(as.matrix(data)))
corstarsl <- function(x){
x <- as.matrix(x)
R <- Hmisc::rcorr(x)$r
p <- Hmisc::rcorr(x)$P
## define notions for significance levels; spacing is important.
mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " ")))
## trunctuate the matrix that holds the correlations to two decimal
R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1]
## build a new matrix that includes the correlations with their apropriate stars
Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
diag(Rnew) <- paste(diag(R), " ", sep="")
rownames(Rnew) <- colnames(x)
colnames(Rnew) <- paste(colnames(x), "", sep="")
## remove upper triangle
Rnew <- as.matrix(Rnew)
Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
Rnew <- as.data.frame(Rnew)
## remove last column and return the matrix (which is now a data frame)
Rnew <- cbind(Rnew[1:length(Rnew)-1])
return(Rnew)
}
### This section parameterizes allometry models and uses the models to estimate unknown stem masses
stemDat$moYr <- paste0(stemDat$mo, "-", substr(stemDat$yr, 3, 4))
stemDat$time <- as.yearmon(stemDat$moYr, "%b-%y")
stemDat <- seasonLabel(data = stemDat)
stemDat$seas <- substr(stemDat$season, 1, 4)
stemDat$plantYr <- as.numeric(paste0(20, substr(stemDat$season, 6, 7))) # avoids splitting winter samples (Jan 2018 is in winter 2017)
### separate plot As for allometry work
stm <- stemDat[stemDat$plot %in% "A", ]
### allometry models
# nls(mass_tot ~ I(a * hgt^b), start = list(a = 0.2, b = 0.2), data = stm[(stm$site %in% "LUM1") & (stm$status %in% "LIVE") & (stm$seas %in% "fall"), ])
# models <- dlply(stm, .(site, season, status), function(x)
# nls(mass_tot ~ I(a * hgt^b), start = list(a = 0.2, b = 0.2), data = x))
# ldply(models, coef)
models.hiRes <- ddply(stm, .(site, seas), function(x)
getParams(dataset = x))
models.reg <- ddply(stm, .(region, seas, plantYr), function(x)
getParams(dataset = x))
models <- ddply(stm, .(region, seas), function(x)
getParams(dataset = x))
dead.allom <- getParams(dataset = stm)
dead.allom$region <- "all"
dead.allom$seas <- "all"
models <- rbind.fill(models, dead.allom)
exponents_vs_coefs <- ggplot(aes(y = exp.live, x = coef.live, colour = seas), data = models) + geom_point() + theme_classic() + labs(y = "exponent", x = "coefficient")
# spring = low coefficients, high exponents
allom_by_seas_plantYr <- ggplot(aes(y = mass_tot, x = hgt, colour = site),
data = stm[(stm$region %in% "LUMCON") & (stm$status %in% "LIVE"), ]) +
geom_point() + theme_classic() + labs(y = "mass", x = "hgt") +
facet_grid(seas ~ plantYr)
allom_by_seas <- ggplot(aes(y = mass_tot, x = hgt, colour = site),
data = stm[(stm$region %in% "LUMCON") & (stm$status %in% "LIVE"), ]) +
geom_point() + theme_classic() + labs(y = "mass", x = "hgt") +
facet_grid(seas ~ .)
# Apply allometry to un-weighed stems -------------------------------------
### Generate estimates from season-region models and pooled data
for (i in 1:nrow(stemDat)) {
# identify region and season
targRegion <- stemDat$region[i]
targetSeas <- stemDat$seas[i]
# find seasonal allometry parameters for the region
if (stemDat$status[i] %in% "LIVE") {
coefficient <- models$coef.live[(models$region %in% targRegion) & (models$seas %in% targetSeas)]
exponent <- models$exp.live[(models$region %in% targRegion) & (models$seas %in% targetSeas)]
} else if (stemDat$status[i] %in% "DEAD") {
coefficient <- models$coef.dead[(models$region %in% "all") & (models$seas %in% "all")]
exponent <- models$exp.dead[(models$region %in% "all") & (models$seas %in% "all")]
}
# apply allometry to get mass estimate
stemDat$mass.est1[i] <- coefficient * (stemDat$hgt[i] ^ exponent) # season-region models for live, pooled model for dead
if (stemDat$status[i] %in% "LIVE") {
stemDat$mass.est2[i] <- models$coef.live[(models$region %in% "all")] * (stemDat$hgt[i] ^ models$exp.live[(models$region %in% "all")]) # single model for all data (always used for dead stems)
} else {
stemDat$mass.est2[i] <- NA
}
# define final mass to be used
# set dummy variable if plant mass was actually measured
if (!is.na(stemDat$mass_tot[i])) {
stemDat$modeled[i] <- 0
stemDat$mass[i] <- stemDat$mass_tot[i]
} else {
stemDat$modeled[i] <- 1
stemDat$mass[i] <- stemDat$mass.est1[i]
}
}
### Identify differences between seasons in the LUMCON region (pooled)
subData <- stm[(stm$status %in% "LIVE") & !is.na(stm$hgt) & !is.na(stm$mass_tot) & (stm$region %in% "LUMCON"), ]
subData$seas <- ordered(subData$seas, c("sprg", "sumr", "fall", "wint"))
full.nls <- nls(mass_tot ~ I(a * hgt^b), data = subData, start = list(a = 0.2, b = 0.2))
summary(full.nls)
subData$resid <- resid(full.nls)
# Residuals have unequal variances among seasons (based on Levene Test and boxplots).
# To make comparisons I follow the procedure described in Herberich et al. 2010
leveneTest(resid ~ seas, data = subData)
boxplot(resid ~ seas, data = subData)
# ANOVA of residuals by seas
modSeas <- aov(resid ~ seas, data = subData)
summary(modSeas) # F = 104.9 p<<0.001 DEAD: F = 3.2; p = 0.02
# glht sets up multiple contrasts. vcov = vcovHC specifies heteroscedasticity-consistent estimator of covariance matrix (Herberich et al. 2010)
modSeas_glht <- glht(modSeas, linfct = mcp(seas = "Tukey"), vcov = vcovHC)
coef(modSeas_glht)
summary(modSeas_glht)
plot(confint(modSeas_glht)) # all seasons are significantly different. Dead stem allometry: no seasonal differences
### This section explores trends in stem size class distributions
### get demographics by plot
int_df <- ddply(stemDat, .(region, site, plot), function(x)
getDemographics(x, quadrat_m2 = plotSize))
### mean by marsh
final_df <- ddply(int_df, .(region, site, group, time), summarise,
stmDens = mean(hgt.count, na.rm = TRUE),
stmDens.se = se(hgt.count),
hgt = mean(hgt.mean, na.rm = TRUE),
hgt.se = se(hgt.mean))
stm.demographics <- ggplot(final_df[final_df$region %in% "LUMCON",],
aes(y = stmDens, x = as.numeric(time))) + #, col = region)) +
geom_point(size = 1.2) + geom_errorbar(aes(ymin = (stmDens - stmDens.se), ymax = (stmDens + stmDens.se)), width = 0) +
geom_line(lwd = 0.7) +
theme_classic() + theme(legend.title = element_blank()) +
ylab(expression("Live stem density (m"^-2~")")) + xlab("") +
facet_wrap("group", nrow = 2) +
theme_bw() %+replace% theme(strip.background = element_blank(), legend.title = element_blank())
### add in plot-level mass data
mass_df <- ddply(stemDat, .(region, site, plot, time), summarise,
mass = mean(mass, na.rm = TRUE),
mass.se = se(mass))
mass_df2 <- join_all(list(mass_df, int_df))
ggplot(mass_df2, aes(x = log(hgt.count), y = log(mass))) + geom_point() + theme_classic()
# identify seasonality
for (i in 1:length(unique(final_df$group)[1:3])) {
test <- ts(final_df[(final_df$region %in% "LUMCON") &
(final_df$group %in% unique(final_df$group)[i]), "stmDens"],
start = c(12, 5), deltat = 1/12)
test2 <- stl(test, s.window = 12)
print(apply(test2$time.series, 2, var) / var(test))
# evaluate change as the difference between start and end values of the trend component
print(noquote(paste0("trend for stems in ", unique(final_df$group)[i], " size class: ",
round(as.numeric(test2$time.series[nrow(test2$time.series), 2] - test2$time.series[1, 2])))))
}
# no obvious trends
Allometry data
allom_by_seas
Allometry models differ significantly by season, and seasonality has become more pronounced since Hill and Roberts (2017).
sub <- stemDat[(stemDat$region %in% "LUMCON") & (stemDat$status %in% "LIVE"), ]
fig2Col <- "gray40" # rgb(169, 169, 169, 140, maxColorValue = 255)
pointSize <- 0.3
tckSize <- 0.018
lwdSize <- 1.5
# png(filename = "allometry.png", width = 120, height = 120, units = "mm", res = 200)
par(mar = c(2, 2, 0.5, 0.5))
plot(sub$mass_tot ~ sub$hgt,
cex = pointSize, pch = 19, col = fig2Col,
ylab = "", xlab = "", xlim = c(0, 172), ylim = c(0, 15.1), xaxt = "n", yaxt = "n",
las = 1, tcl = 0.25, tck = tckSize, bty = "n", yaxs = "i", xaxs = "i")
abline(h = 0, lwd = lwdSize)
abline(v = 0, lwd = 1.5*lwdSize)
axis(side = 2, tcl = 0.25, tck = tckSize, at = axTicks(2), labels = FALSE, lwd = 1)
mtext(side = 2, text = "Mass (g)", line = 1.2)
mtext(side = 2, at = axTicks(2), text = axTicks(2), las = 1, line = 0.35)
axis(side = 1, tcl = 0.25, tck = tckSize, at = axTicks(1), labels = FALSE, lwd = 1)
mtext(text = axTicks(1), at = axTicks(1), side = 1, line = 0.2)
mtext(text = "Height (cm)", side = 1, line = 1.1)
x <- sub$hgt
y.pred2 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "sprg")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "sprg")]
y.pred3 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "sumr")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "sumr")]
y.pred4 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "fall")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "fall")]
y.pred5 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "wint")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "wint")]
lines(x = x[order(y.pred2)], y = y.pred2[order(y.pred2)], lty = 1, lwd = 2)
lines(x = x[order(y.pred3)], y = y.pred3[order(y.pred3)], lty = 2, lwd = 2)
lines(x = x[order(y.pred4)], y = y.pred4[order(y.pred4)], lty = 3, lwd = 2)
lines(x = x[order(y.pred5)], y = y.pred5[order(y.pred5)], lty = 4, lwd = 2)
legend(x = 4, y = 13, legend = c("Spring", "Summer", "Fall", "Winter"),
lty = c(1:4), lwd = 2, bty = 'n', cex = 0.85)
# dev.off()
exponents_vs_coefs
The relationship between allometry exponents and coefficients remains intact and is consistent with Hill and Roberts (2017).
stm.demographics
NAPP
### This section explores NAPP using the bioDat object
bioDat$moYr <- paste0(bioDat$mo, "-", substr(bioDat$yr, 3, 4))
bioDat$monthYear <- zoo::as.yearmon(bioDat$moYr, "%b-%y")
bioDat <- seasonLabel(data = bioDat)
marsh.biomass <- ddply(bioDat, .(region, site, yr, season, monthYear), summarise,
above = mean(mass, na.rm = T),
above.se = se(mass),
stemDensity = mean(stems, na.rm = T),
stemDensity.se = se(stems),
mass.dead = mean(mass.dead, na.rm = T), # stem biomass only
mass.dead.se = se(mass.dead), # stem biomass only
litter = mean(litter, na.rm = TRUE),
litter.se = se(litter),
above.dead = mean(mass.dead, na.rm = T) + mean(litter, na.rm = TRUE), # stems + litter
above.dead.se = sqrt(se(mass.dead)^2 + se(litter)^2),
live.bg = mean(live.bg.gm2, na.rm = TRUE),
live.bg.se = se(live.bg.gm2),
dead.bg = mean(live.bg.gm2, na.rm = TRUE),
dead.bg.se = se(dead.bg.gm2),
tot.bg = mean(tot.bg.gm2, na.rm = TRUE),
tot.bg.se = se(tot.bg.gm2)
)
### plot biomass over time
abv <- ggplot(marsh.biomass[marsh.biomass$region %in% "LUMCON", ], aes(y = above, x = monthYear)) + geom_pointrange(aes(ymax = above + above.se, ymin = above - above.se), size = 0.25) + theme_classic() + facet_grid(site ~ .) + labs(x = "", y = expression("Aboveground biomass (g "%.%m^-2*")"))
### NAPP
# "dead" column in napp is the sum of standing dead and litter
napp.full <- nappCalc(marsh.biomass, summarize = TRUE, EOS = TRUE, EOS_window = 0, annualReset = FALSE)
napp <- regionLabel(napp.full$summary)
napp.trend <- ggplot(napp[grepl(x = napp$site, pattern = "LUM"), ], aes(x = year, y = napp.smalley)) +
geom_bar(stat = "identity") + theme_classic() + facet_grid(. ~ site) +
labs(x = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), y = "")
### correlations between NAPP estimates
NAPP_columns <- c(5:6, 8:9, 17) # columns with NAPP estimates of interest
napp.melt <- melt(napp[grepl(x = napp$site, pattern = "LUM"), c(1, NAPP_columns)], id.vars = c("napp.smalley", "site"))
levels(napp.melt$variable) <- c("Milner and hughes 1968", "Valiela et al. 1975", "Peak biomass", "End-of-season")
napp.correlations <- ggplot(napp.melt, aes(x = napp.smalley, y = value, col = variable)) +
geom_abline(slope = 1, intercept = 0, size = 0.7) +
labs(x = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), y = expression("NAPP (g "%.%m^-2%.%yr^-1~")")) +
theme_classic() + theme(legend.title = element_blank()) +
geom_point(aes(shape = as.factor(site))) + ylim(0, 3800) + xlim(0, 3400) +
geom_smooth(data = napp.melt, method = lm, se = TRUE, linetype = 1, size = 0.7, alpha = 0.2)
### This section explores environmental drivers of biomass and NAPP using the plot-level biomass data (bioDat object) and soil/water data (nutDat)
### NAPP
dd.nut <- ddply(nutDat[nutDat$month.number %in% growingMonths, ], .(year, region, site),
numcolwise(mean, na.rm = TRUE))
napp.nuts <- join_all(list(napp, dd.nut))
### salinity relationships
rsqr.salinity <- round(summary(lm(napp.smalley ~ salinity.psu.bay, data = napp.nuts))$r.squared, 2)
rsqr.salinity.lum <- round(summary(lm(napp.smalley ~ salinity.psu.bay, data = napp.nuts[napp.nuts$region %in% "LUMCON", ]))$r.squared, 2)
napp.salt <- ggplot(napp.nuts, aes(x = salinity.psu.bay, y = napp.smalley)) + geom_point(aes(colour = site)) + theme_classic() + geom_smooth(method = "lm") +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") # all data - looks nice
napp.salt.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = salinity.psu.bay, y = napp.smalley, colour = site)) + geom_point() + theme_classic() + geom_smooth(method = "lm") +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") # weak!
### NO3 relationship
rsqr.no3 <- round(summary(lm(napp.smalley ~ dissolved.no3.no2.bay, data = napp.nuts))$r.squared, 2)
rsqr.no3.lum <- round(summary(lm(napp.smalley ~ dissolved.no3.no2.bay, data = napp.nuts[napp.nuts$region %in% "LUMCON", ]))$r.squared, 2)
napp.no3 <- ggplot(napp.nuts, aes(x = I(log(dissolved.no3.no2.bay)), y = I(log(napp.smalley)))) + geom_point(aes(colour = site)) + theme_classic() + geom_smooth(method = "lm") +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") # all data - looks pretty nice
napp.no3.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = dissolved.no3.no2.bay, y = napp.smalley, colour = site)) + geom_point() + theme_classic() + geom_smooth(method = "lm") +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") # weak
abv
napp.trend
napp.correlations
# ggsave("napp_correlations.png", width = 5.5, height = 4, units = "in")
# table with correlations between NAPP measures (for LUM only, all years)
# kable(corstarsl(as.matrix(napp[grepl(x = napp$site, pattern = "LUM"), NAPP_columns]))) # each site-year counted as independent obs
driver.fig <- arrangeGrob(napp.salt, napp.no3, napp.salt.lum, napp.no3.lum, nrow = 2)
plot(driver.fig)
The relationships between NAPP and bay salinity/NO3 are reasonably strong when Terrebonne Bay sites are included (r2 = 0.46 and 0.23 for salinity and NO3, respectively) but disappear when the dataset is reduced to LUMCON sites (r2 = 0.03 and 0.04 for salinity and NO3, respectively).
NDVI and Aboveground Biomass Relationship
#creating new dataset with 2017-2018 data only (since that is what NDVI data we have)
LUM<-subset(bioDat, region=="LUMCON")
LUM1718<-subset(LUM, yr > 2016)
#creating new dataset of aboveground biomass with means.
mean.mass<-summarySE(data=LUM1718, measurevar="mass", groupvars= c("site", "yr", "mo"))
#this is the NDVI data
# ndvi <- read_csv("ndvi_lum.csv")
# usethis::use_data(ndvi, overwrite = TRUE)
#Merge ndvi and lum data
mydata<-merge(ndvi, mean.mass, by=c("site","yr","mo"))
plot(avg.ndvi~mass, data=mydata)
rsqr.ndvi<- round(summary(lm(mass ~ avg.ndvi, data = mydata))$r.squared, 2)
mass.ndvi <- ggplot(mydata, aes(x = avg.ndvi, y = mass)) + geom_point(aes(colour = site)) + theme_classic() +
geom_smooth(method = "lm") + labs(x = ("NDVI"), y = expression("Aboveground biomass (g "%.%m^-2*")")) # all data - looks nice
##making text larger for figure
exp_vs_coefs <- ggplot(aes(y = exp.live, x = coef.live, colour = seas), data = models) + geom_point(size=3) + theme_classic() + labs(y = "exponent", x = "coefficient") +
theme(text = element_text(size=14))
# spring = low coefficients, high exponents
##salt figure
napp.salt.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = salinity.psu.bay, y = napp.smalley, colour = site)) + geom_point(size=3) + theme_classic() +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") + theme(text = element_text(size=14)) # weak!
##no3 figure
napp.no3.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = dissolved.no3.no2.bay, y = napp.smalley, colour = site)) + geom_point(size=3) + theme_classic() +
labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") + theme(text = element_text(size=14)) # weak
#stem density figure
live.stem<- summarySE(data=LUM, measurevar = "stems", groupvars=c("site", "yr", "mo", "monthYear"))
dead.stem<- summarySE(data=LUM, measurevar = "stems.dead", groupvars=c("site", "yr", "mo", "monthYear"))
mass.ndvi
#rsq is 0.33
ggplot(live.stem,
aes(y = stems, x = as.numeric(monthYear))) + #, col = region)) +
geom_point(size = 1.2) + geom_errorbar(aes(ymin = (stems - se), ymax = (stems + se)), width = 0) +
geom_line(lwd = 0.7) +
theme_classic() + theme(legend.title = element_blank()) +
ylab(expression("Live stem density (m"^-2~")")) + xlab("") +
# facet_wrap("group", nrow = 2) +
theme_bw() %+replace% theme(strip.background = element_blank(), legend.title = element_blank())
troyhill/LUMCONplants documentation built on May 14, 2019, 9:39 a.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.
### This section loads libraries, sets some local variables knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) 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(gridExtra)){ install.packages("gridExtra", 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') } if(!require(multcomp)){ install.packages("multcomp", repos='http://cran.us.r-project.org') } if(!require(sandwich)){ install.packages("sandwich", repos='http://cran.us.r-project.org') } if(!require(car)){ install.packages("car", repos='http://cran.us.r-project.org') } if(!require(Rmisc)){ install.packages("Rmisc", repos='http://cran.us.r-project.org') } if(!require(tidyverse)){ install.packages("tidyverse", repos='http://cran.us.r-project.org') } if(!require(lattice)){ install.packages("lattice", repos='http://cran.us.r-project.org') } if(!require(nlme)){ install.packages("nlme", repos='http://cran.us.r-project.org') } if(!require(lubridate)){ install.packages("lubridate", repos='http://cran.us.r-project.org') } # if(!require(svglite)){ # install.packages("svglite", repos='http://cran.us.r-project.org') # } library(ggplot2) library(gridExtra) library(plyr) library(LUMCONplants) library(multcomp) # for glht library(sandwich) # for HC3 estimator library(car) # for Levene's test library(Rmisc) library(tidyverse) library(lattice) library(nlme) library(lubridate) # library(svglite) # doesn't install properly on linux ### define some constants plotSize <- 0.25 * 0.25 # size of quadrat, m2 coreTube <- pi*(6.9/2)^2 # area of plastic coring tube used prior to August 2016: units = cm2 coreTube_metal <- pi*(5.8/2)^2 # area of metal coring device used beginning in August 2016: units = cm2 growingMonths <- 4:10 # growing season defined as April:October bagMass <- 5.4 # grams; bags used for belowground cores tckSize <- 0.01 getDemographics <- function(data = stemDat, status = "LIVE", groupBreaks = c(0, 25, 50, 100, 200), groupLabels = c("0-25 cm", "25-50 cm", "50-100 cm", ">100 cm"), quadrat_m2 = 1) { df1 <- transform(data[(data$status %in% status), ], group = cut(hgt, breaks = groupBreaks, labels = groupLabels)) res <- do.call(data.frame, aggregate(hgt ~ time + group, df1, FUN = function(x) c(count = sum(!is.na(x)) / quadrat_m2, mean = mean(x, na.rm = TRUE) / quadrat_m2))) dNew <- data.frame(group = levels(df1$group)) df2 <- merge(res, dNew, all = TRUE) # add missing bins, set to 0 hgt/dens missedLevs <- expand.grid(time = unique(res$time), group = unique(res$group), hgt.count = 0, hgt.mean = 0) missingLevs <- missedLevs[which(!paste0(missedLevs$time, missedLevs$group) %in% paste0(res$time, res$group)), ] df2 <- rbind(df2, missingLevs) #df2$site <- targReg invisible(df2) } regionLabel <- function(data, siteCol = "site") { # function takes a dataset and the name of the column with site names (e.g., "LUM1"), and # adds a column with marsh names ("LUM", "TB-A", "TB-B") data$region <- as.character(NA) for (i in 1:nrow(data)) { if (data[i, siteCol] %in% paste0("LUM", 1:3)) { data$region[i] <- "LUMCON" } else if (data[i, siteCol] %in% c("TB1", "TB2")) { data$region[i] <- "Bay LaFleur" } else if (data[i, siteCol] %in% c("TB3", "TB4")) { data$region[i] <- "Lake Barre" } } data } ### a function to make correlation matrices with significance stars ### from http://myowelt.blogspot.com/2008/04/beautiful-correlation-tables-in-r.html ### usage: kable(corstarsl(as.matrix(data))) corstarsl <- function(x){ x <- as.matrix(x) R <- Hmisc::rcorr(x)$r p <- Hmisc::rcorr(x)$P ## define notions for significance levels; spacing is important. mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " "))) ## trunctuate the matrix that holds the correlations to two decimal R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1] ## build a new matrix that includes the correlations with their apropriate stars Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x)) diag(Rnew) <- paste(diag(R), " ", sep="") rownames(Rnew) <- colnames(x) colnames(Rnew) <- paste(colnames(x), "", sep="") ## remove upper triangle Rnew <- as.matrix(Rnew) Rnew[upper.tri(Rnew, diag = TRUE)] <- "" Rnew <- as.data.frame(Rnew) ## remove last column and return the matrix (which is now a data frame) Rnew <- cbind(Rnew[1:length(Rnew)-1]) return(Rnew) }
### This section parameterizes allometry models and uses the models to estimate unknown stem masses stemDat$moYr <- paste0(stemDat$mo, "-", substr(stemDat$yr, 3, 4)) stemDat$time <- as.yearmon(stemDat$moYr, "%b-%y") stemDat <- seasonLabel(data = stemDat) stemDat$seas <- substr(stemDat$season, 1, 4) stemDat$plantYr <- as.numeric(paste0(20, substr(stemDat$season, 6, 7))) # avoids splitting winter samples (Jan 2018 is in winter 2017) ### separate plot As for allometry work stm <- stemDat[stemDat$plot %in% "A", ] ### allometry models # nls(mass_tot ~ I(a * hgt^b), start = list(a = 0.2, b = 0.2), data = stm[(stm$site %in% "LUM1") & (stm$status %in% "LIVE") & (stm$seas %in% "fall"), ]) # models <- dlply(stm, .(site, season, status), function(x) # nls(mass_tot ~ I(a * hgt^b), start = list(a = 0.2, b = 0.2), data = x)) # ldply(models, coef) models.hiRes <- ddply(stm, .(site, seas), function(x) getParams(dataset = x)) models.reg <- ddply(stm, .(region, seas, plantYr), function(x) getParams(dataset = x)) models <- ddply(stm, .(region, seas), function(x) getParams(dataset = x)) dead.allom <- getParams(dataset = stm) dead.allom$region <- "all" dead.allom$seas <- "all" models <- rbind.fill(models, dead.allom) exponents_vs_coefs <- ggplot(aes(y = exp.live, x = coef.live, colour = seas), data = models) + geom_point() + theme_classic() + labs(y = "exponent", x = "coefficient") # spring = low coefficients, high exponents allom_by_seas_plantYr <- ggplot(aes(y = mass_tot, x = hgt, colour = site), data = stm[(stm$region %in% "LUMCON") & (stm$status %in% "LIVE"), ]) + geom_point() + theme_classic() + labs(y = "mass", x = "hgt") + facet_grid(seas ~ plantYr) allom_by_seas <- ggplot(aes(y = mass_tot, x = hgt, colour = site), data = stm[(stm$region %in% "LUMCON") & (stm$status %in% "LIVE"), ]) + geom_point() + theme_classic() + labs(y = "mass", x = "hgt") + facet_grid(seas ~ .) # Apply allometry to un-weighed stems ------------------------------------- ### Generate estimates from season-region models and pooled data for (i in 1:nrow(stemDat)) { # identify region and season targRegion <- stemDat$region[i] targetSeas <- stemDat$seas[i] # find seasonal allometry parameters for the region if (stemDat$status[i] %in% "LIVE") { coefficient <- models$coef.live[(models$region %in% targRegion) & (models$seas %in% targetSeas)] exponent <- models$exp.live[(models$region %in% targRegion) & (models$seas %in% targetSeas)] } else if (stemDat$status[i] %in% "DEAD") { coefficient <- models$coef.dead[(models$region %in% "all") & (models$seas %in% "all")] exponent <- models$exp.dead[(models$region %in% "all") & (models$seas %in% "all")] } # apply allometry to get mass estimate stemDat$mass.est1[i] <- coefficient * (stemDat$hgt[i] ^ exponent) # season-region models for live, pooled model for dead if (stemDat$status[i] %in% "LIVE") { stemDat$mass.est2[i] <- models$coef.live[(models$region %in% "all")] * (stemDat$hgt[i] ^ models$exp.live[(models$region %in% "all")]) # single model for all data (always used for dead stems) } else { stemDat$mass.est2[i] <- NA } # define final mass to be used # set dummy variable if plant mass was actually measured if (!is.na(stemDat$mass_tot[i])) { stemDat$modeled[i] <- 0 stemDat$mass[i] <- stemDat$mass_tot[i] } else { stemDat$modeled[i] <- 1 stemDat$mass[i] <- stemDat$mass.est1[i] } }
### Identify differences between seasons in the LUMCON region (pooled) subData <- stm[(stm$status %in% "LIVE") & !is.na(stm$hgt) & !is.na(stm$mass_tot) & (stm$region %in% "LUMCON"), ] subData$seas <- ordered(subData$seas, c("sprg", "sumr", "fall", "wint")) full.nls <- nls(mass_tot ~ I(a * hgt^b), data = subData, start = list(a = 0.2, b = 0.2)) summary(full.nls) subData$resid <- resid(full.nls) # Residuals have unequal variances among seasons (based on Levene Test and boxplots). # To make comparisons I follow the procedure described in Herberich et al. 2010 leveneTest(resid ~ seas, data = subData) boxplot(resid ~ seas, data = subData) # ANOVA of residuals by seas modSeas <- aov(resid ~ seas, data = subData) summary(modSeas) # F = 104.9 p<<0.001 DEAD: F = 3.2; p = 0.02 # glht sets up multiple contrasts. vcov = vcovHC specifies heteroscedasticity-consistent estimator of covariance matrix (Herberich et al. 2010) modSeas_glht <- glht(modSeas, linfct = mcp(seas = "Tukey"), vcov = vcovHC) coef(modSeas_glht) summary(modSeas_glht) plot(confint(modSeas_glht)) # all seasons are significantly different. Dead stem allometry: no seasonal differences
### This section explores trends in stem size class distributions ### get demographics by plot int_df <- ddply(stemDat, .(region, site, plot), function(x) getDemographics(x, quadrat_m2 = plotSize)) ### mean by marsh final_df <- ddply(int_df, .(region, site, group, time), summarise, stmDens = mean(hgt.count, na.rm = TRUE), stmDens.se = se(hgt.count), hgt = mean(hgt.mean, na.rm = TRUE), hgt.se = se(hgt.mean)) stm.demographics <- ggplot(final_df[final_df$region %in% "LUMCON",], aes(y = stmDens, x = as.numeric(time))) + #, col = region)) + geom_point(size = 1.2) + geom_errorbar(aes(ymin = (stmDens - stmDens.se), ymax = (stmDens + stmDens.se)), width = 0) + geom_line(lwd = 0.7) + theme_classic() + theme(legend.title = element_blank()) + ylab(expression("Live stem density (m"^-2~")")) + xlab("") + facet_wrap("group", nrow = 2) + theme_bw() %+replace% theme(strip.background = element_blank(), legend.title = element_blank()) ### add in plot-level mass data mass_df <- ddply(stemDat, .(region, site, plot, time), summarise, mass = mean(mass, na.rm = TRUE), mass.se = se(mass)) mass_df2 <- join_all(list(mass_df, int_df)) ggplot(mass_df2, aes(x = log(hgt.count), y = log(mass))) + geom_point() + theme_classic() # identify seasonality for (i in 1:length(unique(final_df$group)[1:3])) { test <- ts(final_df[(final_df$region %in% "LUMCON") & (final_df$group %in% unique(final_df$group)[i]), "stmDens"], start = c(12, 5), deltat = 1/12) test2 <- stl(test, s.window = 12) print(apply(test2$time.series, 2, var) / var(test)) # evaluate change as the difference between start and end values of the trend component print(noquote(paste0("trend for stems in ", unique(final_df$group)[i], " size class: ", round(as.numeric(test2$time.series[nrow(test2$time.series), 2] - test2$time.series[1, 2]))))) } # no obvious trends
Allometry data
allom_by_seas
Allometry models differ significantly by season, and seasonality has become more pronounced since Hill and Roberts (2017).
sub <- stemDat[(stemDat$region %in% "LUMCON") & (stemDat$status %in% "LIVE"), ] fig2Col <- "gray40" # rgb(169, 169, 169, 140, maxColorValue = 255) pointSize <- 0.3 tckSize <- 0.018 lwdSize <- 1.5 # png(filename = "allometry.png", width = 120, height = 120, units = "mm", res = 200) par(mar = c(2, 2, 0.5, 0.5)) plot(sub$mass_tot ~ sub$hgt, cex = pointSize, pch = 19, col = fig2Col, ylab = "", xlab = "", xlim = c(0, 172), ylim = c(0, 15.1), xaxt = "n", yaxt = "n", las = 1, tcl = 0.25, tck = tckSize, bty = "n", yaxs = "i", xaxs = "i") abline(h = 0, lwd = lwdSize) abline(v = 0, lwd = 1.5*lwdSize) axis(side = 2, tcl = 0.25, tck = tckSize, at = axTicks(2), labels = FALSE, lwd = 1) mtext(side = 2, text = "Mass (g)", line = 1.2) mtext(side = 2, at = axTicks(2), text = axTicks(2), las = 1, line = 0.35) axis(side = 1, tcl = 0.25, tck = tckSize, at = axTicks(1), labels = FALSE, lwd = 1) mtext(text = axTicks(1), at = axTicks(1), side = 1, line = 0.2) mtext(text = "Height (cm)", side = 1, line = 1.1) x <- sub$hgt y.pred2 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "sprg")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "sprg")] y.pred3 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "sumr")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "sumr")] y.pred4 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "fall")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "fall")] y.pred5 <- models$coef.live[(models$region %in% "LUMCON") & (models$seas %in% "wint")] * x^models$exp.live[(models$region %in% "LUMCON") & (models$seas %in% "wint")] lines(x = x[order(y.pred2)], y = y.pred2[order(y.pred2)], lty = 1, lwd = 2) lines(x = x[order(y.pred3)], y = y.pred3[order(y.pred3)], lty = 2, lwd = 2) lines(x = x[order(y.pred4)], y = y.pred4[order(y.pred4)], lty = 3, lwd = 2) lines(x = x[order(y.pred5)], y = y.pred5[order(y.pred5)], lty = 4, lwd = 2) legend(x = 4, y = 13, legend = c("Spring", "Summer", "Fall", "Winter"), lty = c(1:4), lwd = 2, bty = 'n', cex = 0.85) # dev.off()
exponents_vs_coefs
The relationship between allometry exponents and coefficients remains intact and is consistent with Hill and Roberts (2017).
stm.demographics
NAPP
### This section explores NAPP using the bioDat object bioDat$moYr <- paste0(bioDat$mo, "-", substr(bioDat$yr, 3, 4)) bioDat$monthYear <- zoo::as.yearmon(bioDat$moYr, "%b-%y") bioDat <- seasonLabel(data = bioDat) marsh.biomass <- ddply(bioDat, .(region, site, yr, season, monthYear), summarise, above = mean(mass, na.rm = T), above.se = se(mass), stemDensity = mean(stems, na.rm = T), stemDensity.se = se(stems), mass.dead = mean(mass.dead, na.rm = T), # stem biomass only mass.dead.se = se(mass.dead), # stem biomass only litter = mean(litter, na.rm = TRUE), litter.se = se(litter), above.dead = mean(mass.dead, na.rm = T) + mean(litter, na.rm = TRUE), # stems + litter above.dead.se = sqrt(se(mass.dead)^2 + se(litter)^2), live.bg = mean(live.bg.gm2, na.rm = TRUE), live.bg.se = se(live.bg.gm2), dead.bg = mean(live.bg.gm2, na.rm = TRUE), dead.bg.se = se(dead.bg.gm2), tot.bg = mean(tot.bg.gm2, na.rm = TRUE), tot.bg.se = se(tot.bg.gm2) ) ### plot biomass over time abv <- ggplot(marsh.biomass[marsh.biomass$region %in% "LUMCON", ], aes(y = above, x = monthYear)) + geom_pointrange(aes(ymax = above + above.se, ymin = above - above.se), size = 0.25) + theme_classic() + facet_grid(site ~ .) + labs(x = "", y = expression("Aboveground biomass (g "%.%m^-2*")")) ### NAPP # "dead" column in napp is the sum of standing dead and litter napp.full <- nappCalc(marsh.biomass, summarize = TRUE, EOS = TRUE, EOS_window = 0, annualReset = FALSE) napp <- regionLabel(napp.full$summary) napp.trend <- ggplot(napp[grepl(x = napp$site, pattern = "LUM"), ], aes(x = year, y = napp.smalley)) + geom_bar(stat = "identity") + theme_classic() + facet_grid(. ~ site) + labs(x = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), y = "") ### correlations between NAPP estimates NAPP_columns <- c(5:6, 8:9, 17) # columns with NAPP estimates of interest napp.melt <- melt(napp[grepl(x = napp$site, pattern = "LUM"), c(1, NAPP_columns)], id.vars = c("napp.smalley", "site")) levels(napp.melt$variable) <- c("Milner and hughes 1968", "Valiela et al. 1975", "Peak biomass", "End-of-season") napp.correlations <- ggplot(napp.melt, aes(x = napp.smalley, y = value, col = variable)) + geom_abline(slope = 1, intercept = 0, size = 0.7) + labs(x = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), y = expression("NAPP (g "%.%m^-2%.%yr^-1~")")) + theme_classic() + theme(legend.title = element_blank()) + geom_point(aes(shape = as.factor(site))) + ylim(0, 3800) + xlim(0, 3400) + geom_smooth(data = napp.melt, method = lm, se = TRUE, linetype = 1, size = 0.7, alpha = 0.2)
### This section explores environmental drivers of biomass and NAPP using the plot-level biomass data (bioDat object) and soil/water data (nutDat) ### NAPP dd.nut <- ddply(nutDat[nutDat$month.number %in% growingMonths, ], .(year, region, site), numcolwise(mean, na.rm = TRUE)) napp.nuts <- join_all(list(napp, dd.nut)) ### salinity relationships rsqr.salinity <- round(summary(lm(napp.smalley ~ salinity.psu.bay, data = napp.nuts))$r.squared, 2) rsqr.salinity.lum <- round(summary(lm(napp.smalley ~ salinity.psu.bay, data = napp.nuts[napp.nuts$region %in% "LUMCON", ]))$r.squared, 2) napp.salt <- ggplot(napp.nuts, aes(x = salinity.psu.bay, y = napp.smalley)) + geom_point(aes(colour = site)) + theme_classic() + geom_smooth(method = "lm") + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") # all data - looks nice napp.salt.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = salinity.psu.bay, y = napp.smalley, colour = site)) + geom_point() + theme_classic() + geom_smooth(method = "lm") + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") # weak! ### NO3 relationship rsqr.no3 <- round(summary(lm(napp.smalley ~ dissolved.no3.no2.bay, data = napp.nuts))$r.squared, 2) rsqr.no3.lum <- round(summary(lm(napp.smalley ~ dissolved.no3.no2.bay, data = napp.nuts[napp.nuts$region %in% "LUMCON", ]))$r.squared, 2) napp.no3 <- ggplot(napp.nuts, aes(x = I(log(dissolved.no3.no2.bay)), y = I(log(napp.smalley)))) + geom_point(aes(colour = site)) + theme_classic() + geom_smooth(method = "lm") + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") # all data - looks pretty nice napp.no3.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = dissolved.no3.no2.bay, y = napp.smalley, colour = site)) + geom_point() + theme_classic() + geom_smooth(method = "lm") + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") # weak
abv
napp.trend
napp.correlations
# ggsave("napp_correlations.png", width = 5.5, height = 4, units = "in")
# table with correlations between NAPP measures (for LUM only, all years) # kable(corstarsl(as.matrix(napp[grepl(x = napp$site, pattern = "LUM"), NAPP_columns]))) # each site-year counted as independent obs
driver.fig <- arrangeGrob(napp.salt, napp.no3, napp.salt.lum, napp.no3.lum, nrow = 2) plot(driver.fig)
The relationships between NAPP and bay salinity/NO3 are reasonably strong when Terrebonne Bay sites are included (r2 = 0.46 and 0.23 for salinity and NO3, respectively) but disappear when the dataset is reduced to LUMCON sites (r2 = 0.03 and 0.04 for salinity and NO3, respectively).
NDVI and Aboveground Biomass Relationship
#creating new dataset with 2017-2018 data only (since that is what NDVI data we have) LUM<-subset(bioDat, region=="LUMCON") LUM1718<-subset(LUM, yr > 2016) #creating new dataset of aboveground biomass with means. mean.mass<-summarySE(data=LUM1718, measurevar="mass", groupvars= c("site", "yr", "mo")) #this is the NDVI data # ndvi <- read_csv("ndvi_lum.csv") # usethis::use_data(ndvi, overwrite = TRUE) #Merge ndvi and lum data mydata<-merge(ndvi, mean.mass, by=c("site","yr","mo")) plot(avg.ndvi~mass, data=mydata) rsqr.ndvi<- round(summary(lm(mass ~ avg.ndvi, data = mydata))$r.squared, 2) mass.ndvi <- ggplot(mydata, aes(x = avg.ndvi, y = mass)) + geom_point(aes(colour = site)) + theme_classic() + geom_smooth(method = "lm") + labs(x = ("NDVI"), y = expression("Aboveground biomass (g "%.%m^-2*")")) # all data - looks nice ##making text larger for figure exp_vs_coefs <- ggplot(aes(y = exp.live, x = coef.live, colour = seas), data = models) + geom_point(size=3) + theme_classic() + labs(y = "exponent", x = "coefficient") + theme(text = element_text(size=14)) # spring = low coefficients, high exponents ##salt figure napp.salt.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = salinity.psu.bay, y = napp.smalley, colour = site)) + geom_point(size=3) + theme_classic() + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creek salinity (psu)") + theme(text = element_text(size=14)) # weak! ##no3 figure napp.no3.lum <- ggplot(napp.nuts[napp.nuts$region %in% "LUMCON", ], aes(x = dissolved.no3.no2.bay, y = napp.smalley, colour = site)) + geom_point(size=3) + theme_classic() + labs(y = expression("Smalley NAPP (g "%.%m^-2%.%yr^-1~")"), x = "Creekwater NO3 (uM)") + theme(text = element_text(size=14)) # weak #stem density figure live.stem<- summarySE(data=LUM, measurevar = "stems", groupvars=c("site", "yr", "mo", "monthYear")) dead.stem<- summarySE(data=LUM, measurevar = "stems.dead", groupvars=c("site", "yr", "mo", "monthYear"))
mass.ndvi
#rsq is 0.33
ggplot(live.stem, aes(y = stems, x = as.numeric(monthYear))) + #, col = region)) + geom_point(size = 1.2) + geom_errorbar(aes(ymin = (stems - se), ymax = (stems + se)), width = 0) + geom_line(lwd = 0.7) + theme_classic() + theme(legend.title = element_blank()) + ylab(expression("Live stem density (m"^-2~")")) + xlab("") + # facet_wrap("group", nrow = 2) + theme_bw() %+replace% theme(strip.background = element_blank(), legend.title = element_blank())
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.