library(tidyverse)
library(ggformula)
library(lme4)
library(car)
library(broom)
library(cowplot)
options(scipen = 999)
library(here)
library(bbmle)
library(latex2exp)
library(emmeans)
library(glue)
(licor <- read_rds(here("cleaned_data", "photosynthesis data.rds")))
(spad <- read_rds(here("cleaned_data", "SPAD data.rds")))
(treatment <- read_rds(here("cleaned_data", "BACE_Treatments.rds")))
spad
contains data collected using the SPAD-502 chlorophyll meter produced by Spectrum Technologies, INC (Plainfield, IL). I measured chlorophyll content on three leaves of each plant in the data frame, which is why each plant.id
is repeated three times. SPAD is measured in SPAD units which is calculated as a ratio of intensities of two wavelengths of light (see manual for more detail). Each SPAD data point is a mean of three measurements, averaged using the instrument itself.
plot.id
: Which plot in the BACE experiment the plant was grown inplant.id
: The ID number of the plant (1-11) in each plotSPAD
: Mean chlorophyll content in SPAD units for a leaf (n = 3 locations on the leaf)Treatment
: Which drought treatment receivedlicor
contains data collected using a LI-COR LI6400 portable photosynthesis meter. Data import and cleaning are done in LI-COR data wrangling.Rmd.
plot.id
: Which plot in the BACE experiment the plant was grown inplant.id
: The ID number of the plant (1-11) in each plotleaf
: Which of three leaves (a, b, c) the data were collected from...
: The rest of the columns are explained in documentation for the LI6400 (LICOR System Variables.pdf)spad.lmer <- lmer(SPAD ~ Treatment + (1|plot.id) + (1|unique.id), data = spad)
plot(spad.lmer)
shapiro.test(augment(spad.lmer)$.resid)
Shapiro-Wilk normality test
data: augment(spad.lmer)$.resid
W = 0.98416, p-value = 0.004362
tidy(spad.lmer)
binding factor and character vector, coercing into character vectorbinding character and factor vector, coercing into character vector
fixef(spad.lmer)
(Intercept) Treatment75% Treatment50%
54.363441 -1.308996 -1.120912
summary(spad.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: SPAD ~ Treatment + (1 | plot.id) + (1 | unique.id)
Data: spad
REML criterion at convergence: 2137.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.90994 -0.62705 -0.04793 0.60822 2.45291
Random effects:
Groups Name Variance Std.Dev.
unique.id (Intercept) 69.721054022779626 8.349913414
plot.id (Intercept) 0.000000000003814 0.000001953
Residual 120.134444096770594 10.960585938
Number of obs: 270, groups: unique.id, 90; plot.id, 9
Fixed effects:
Estimate Std. Error t value
(Intercept) 54.363 1.882 28.890
Treatment75% -1.309 2.683 -0.488
Treatment50% -1.121 2.707 -0.414
Correlation of Fixed Effects:
(Intr) Trt75%
Treatmnt75% -0.701
Treatmnt50% -0.695 0.488
It looks like the effect of plot is very small relative to the effect of plant.
Remove random effect of plot, then both random effects. Then compare all to null model.
spad.lmer.2 <- lmer(SPAD ~ Treatment + (1|unique.id), data = spad)
spad.lm <- lm(SPAD ~ Treatment, data = spad)
spad.null <- lm(SPAD ~ 1, data = spad)
AICtab(spad.lmer, spad.lmer.2, spad.lm, spad.null)
dAIC df
spad.lmer.2 0.0 5
spad.lmer 2.0 6
spad.null 34.5 2
spad.lm 38.0 4
spad.posthoc <- emmeans(spad.lmer.2, pairwise ~ Treatment)
Anova(spad.lmer.2)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: SPAD
Chisq Df Pr(>Chisq)
Treatment 0.2788 2 0.8699
spad.lmer.2
wins. This model includes treatment as a main effect and unique plant ID as a random effect. However,the main effect of drought treatment is not significant.
spad.errbar <- spad.posthoc$emmeans %>% tidy()
# summarize data by plant for plots
spad.plotdata <- spad %>%
group_by(Treatment, unique.id) %>%
summarize(SPAD = mean(SPAD))
spad.bar <- ggplot(spad.plotdata, aes(x = Treatment, y = SPAD)) +
geom_jitter(alpha = 0.5, width = 0.2) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high,x = Treatment, y = estimate),
width = 0.2,
data = spad.errbar) +
geom_col(aes(x = Treatment, y = estimate), alpha = 0.6, data = spad.errbar) +
geom_label(aes(x = 3, y = 90, label ="p = 0.870"), data = data.frame()) +
ylab("Chlorophyll Content (SPAD units)") +
xlab("% Ambient Rainfall") +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
Ignoring unknown aesthetics: y
spad.bar
To start, I need to do a little bit of data wrangling since I took multiple measurements on each leaf. These are just technical replicates and represent only instrument variation, so I’m going to average them. I’m going to try using summarize_if
to get means for all the columns, then remove some columns where averages don’t make sense (e.g. Obs
which is just observation number)
#ammend treatment data to licor data. Spray data isn't relevant since photosynthesis done before spraying.
licor2 <- left_join(licor, treatment, by = c("plot.id" = "Plot", "plant.id" = "Plant")) %>%
select(-Spray) %>%
select(Treatment, everything()) %>%
#some of the Treatment entries are "NA" because different plants were used for photosynthesis measurement than were used for spraying and chemistry. Fill those in.
fill(Treatment)
licor.avg <- licor2 %>%
group_by(Treatment, plot.id, plant.id, leaf) %>%
summarise_if(is.double, mean) %>%
select(-Obs, -Date.Time, -FTime, -`EBal?`) %>%
#create unique.id column
mutate(unique.id = paste0(plot.id, plant.id)) %>%
select(Treatment, plot.id, plant.id, unique.id, leaf, everything())
At this point leaf
represents biological replicates, plant.id
and plot.id
are random effects and Treatment
is our only fixed effect.
I’m mostly interested in photo
which is photosynthetic rate measured in units of \(\frac{\mu \textrm{mol CO}_2 /\textrm{s}}{\textrm{m}^2}\).
photo.lmer <- lmer(Photo ~ Treatment + (1|plot.id) + (1|unique.id), data = licor.avg)
tidy(photo.lmer)
binding factor and character vector, coercing into character vectorbinding character and factor vector, coercing into character vector
summary(photo.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: Photo ~ Treatment + (1 | plot.id) + (1 | unique.id)
Data: licor.avg
REML criterion at convergence: 271.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.3434 -0.5237 -0.0152 0.4390 1.7885
Random effects:
Groups Name Variance Std.Dev.
unique.id (Intercept) 1.9733 1.4048
plot.id (Intercept) 0.3694 0.6078
Residual 0.9062 0.9519
Number of obs: 80, groups: unique.id, 27; plot.id, 9
Fixed effects:
Estimate Std. Error t value
(Intercept) 6.1003 0.6132 9.949
Treatment75% -1.9435 0.8671 -2.241
Treatment50% -3.4175 0.8682 -3.936
Correlation of Fixed Effects:
(Intr) Trt75%
Treatmnt75% -0.707
Treatmnt50% -0.706 0.499
Model diagnostics
plot(photo.lmer)
augment(photo.lmer) %>%
ggplot(aes(sample = .resid)) +
geom_qq() +
stat_qqline()+
coord_flip() +
theme_bw()
Looks great
Model competition
#remove plot effect
photo.lmer2 <- lmer(Photo ~ Treatment + (1|unique.id), data = licor.avg)
# Only fixed effects
photo.lm <- lm(Photo ~ Treatment, data = licor.avg)
#null model
photo.null <- lm(Photo ~ 1, data = licor.avg)
AICtab(photo.lm, photo.lmer, photo.lmer2, photo.null)
dAIC df
photo.lmer2 0.0 5
photo.lmer 1.6 6
photo.lm 38.0 4
photo.null 74.7 2
photo.lmer2
wins. This is the same specification as the SPAD model—a random effect of plant, but not plot.
Anova(photo.lmer2)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Photo
Chisq Df Pr(>Chisq)
Treatment 20.646 2 0.00003287 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
photo.posthoc <- emmeans(photo.lmer2, pairwise ~ Treatment)
photo.posthoc$contrasts
contrast estimate SE df t.ratio p.value
100% - 75% 1.943454 0.7528227 23.90 2.582 0.0419
100% - 50% 3.414586 0.7540121 24.04 4.529 0.0004
75% - 50% 1.471132 0.7540121 24.04 1.951 0.1463
P value adjustment: tukey method for comparing a family of 3 estimates
photo.lmer2.noint <- lmer(Photo ~ -1+ Treatment + (1|unique.id), data = licor.avg)
fixef(photo.lmer2.noint)
Treatment100% Treatment75% Treatment50%
6.100260 4.156806 2.685674
Significant effect of drought treatment. 100 = a; 75 = b; 50 = b
letters_data <- licor.avg %>%
group_by(Treatment) %>%
summarise(max.p = max(Photo), max.c = max(Cond))
photo.errbar <- photo.posthoc$emmeans %>% tidy()
# photo.errbar
#summarize data for plots
licor.plotdata <- licor.avg %>%
group_by(Treatment, unique.id) %>%
summarize(Photo = mean(Photo),
Cond = mean(Cond))
#jitter plot with mean ± CI
photo.bar <- ggplot(licor.plotdata, aes(x = Treatment, y = Photo)) +
geom_jitter(alpha = 0.5, width = 0.2) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high,x = Treatment, y = estimate),
width = 0.2,
data = photo.errbar) +
geom_col(aes(x = Treatment, y = estimate), alpha = 0.6, data = photo.errbar) +
geom_text(aes(y = max.p, label = c("a", "b", "b")), vjust = -1, data = letters_data, size = 4) +
geom_label(aes(x = 3, y = 10, label = "p < 0.001"), data = data.frame()) +
ylim(0, 10) +
xlab("% Ambient Rainfall") +
ylab(TeX("Net Assimilation Rate $(\\mu$ mol $CO_{2}$ $m^{-2}$ $s^{-1})$")) +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
Ignoring unknown aesthetics: y
photo.bar
cond.lmer <- lmer(log(Cond) ~ Treatment + (1|plot.id) + (1|unique.id), data = licor.avg)
lmer(log(Cond) ~ -1 + Treatment + (1|plot.id) + (1|unique.id), data = licor.avg) %>% fixef() %>% exp()
Treatment100% Treatment75% Treatment50%
0.05532703 0.02877802 0.01721222
tidy(cond.lmer)
binding factor and character vector, coercing into character vectorbinding character and factor vector, coercing into character vector
Model diagnostics
plot(cond.lmer)
augment(cond.lmer) %>%
ggplot(aes(sample = .resid)) +
geom_qq() +
stat_qqline()+
coord_flip() +
theme_bw()
Some slight pattern in residuals plot, qq plot looks fine.
Model competition
#remove plot effect
cond.lmer2 <- lmer(log(Cond) ~ Treatment + (1|unique.id), data = licor.avg)
cond.lmer3 <- lmer(log(Cond) ~ 1 + (1|unique.id), data = licor.avg)
# fixed efects only
cond.lm <- lm(log(Cond) ~ Treatment, data = licor.avg)
#null model
cond.null <- lm(log(Cond) ~ 1, data = licor.avg)
AICtab(cond.lm, cond.lmer, cond.lmer2, cond.null, cond.lmer3)
dAIC df
cond.lmer 0.0 6
cond.lmer2 2.9 5
cond.lmer3 10.1 3
cond.lm 42.7 4
cond.null 71.1 2
cond.lmer
is the best model, which takes plot and plant into account as random effects and drought treatment as a fixed effect.
Anova(cond.lmer)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(Cond)
Chisq Df Pr(>Chisq)
Treatment 6.6552 2 0.03588 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
cond.posthoc <- emmeans(cond.lmer, pairwise ~ Treatment)
cond.posthoc$contrasts
contrast estimate SE df t.ratio p.value
100% - 75% 0.6536496 0.4534570 5.99 1.441 0.3803
100% - 50% 1.1676417 0.4537173 6.00 2.574 0.0929
75% - 50% 0.5139921 0.4537173 6.00 1.133 0.5306
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates
Treatment is significant. Anova of cond.lmer2
shows treatment as MUCH more significant. Must mean there is an interaction between plot and treatment. Post-hoc tests show no significant pairwise differences.
This could also be done by modeling with a log-link rather than log transforming. However, the advantage of log transforming is that it takes care of the linearity assumption and the heterosketasticity assumption which are both violated by these data (I think). Using a log-link only takes care of the lineartiy assumption because the it transforms the means instead of taking the means of transformed data. More info here: https://stats.stackexchange.com/questions/47840/linear-model-with-log-transformed-response-vs-generalized-linear-model-with-log
# build models
cond.glmer <- glmer(Cond ~ Treatment + (1|plot.id) + (1|unique.id), family = gaussian(link = "log"), data = licor.avg)
cond.glmer2 <- glmer(Cond ~ Treatment + (1|unique.id), family = gaussian(link = "log"), data = licor.avg)
cond.glmer3 <- glmer(Cond ~ Treatment + (1|plot.id), family = gaussian(link = "log"), data = licor.avg)
cond.glm <- glm(Cond ~ Treatment, family = gaussian(link = "log"), data = licor.avg)
cond.glmer4 <- glmer(Cond ~ 1 + (1|plot.id) + (1|unique.id), family = gaussian(link = "log"), data = licor.avg)
cond.gnull <- glm(Cond ~ 1, family = gaussian(link = "log"), data = licor.avg)
#model competition
AICtab(cond.lmer, cond.glmer, cond.glmer2, cond.glmer3, cond.glm, cond.glmer4, cond.gnull)
#plot diagnostics
Anova(cond.glmer2)
plot(cond.glmer2)
augment(cond.glmer2) %>%
ggplot(aes(sample = .resid)) +
geom_qq() +
stat_qqline()+
coord_flip() +
theme_bw()
#remake bar plot
test <- emmeans(cond.glmer2, pairwise~Treatment)
test2 <- confint(test)$emmeans
test2
ggplot(licor.avg, aes(x = Treatment, y = Cond)) +
geom_jitter(alpha = 0.5, width = 0.2) +
geom_errorbar(aes(ymin = exp(asymp.LCL), ymax = exp(asymp.UCL),x = Treatment, y = exp(emmean)),
width = 0.2,
data = test2) +
geom_col(aes(x = Treatment, y = exp(emmean)), alpha = 0.6, data = test2) +
geom_label(aes(x = 3, y = 0.15, label = paste0("p = ", Anova(cond.glmer2)$`Pr(>Chisq)` %>% round(3))), data = data.frame()) +
ylim(0, 0.15) +
xlab("% Ambient Rainfall") +
ylab(TeX("Stomatal conductance $($mol $H_{2}O$ $m^{-2}$ $s^{-1})$")) +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
Results in tighter error bars and non-significant ANOVA (in fact, the best model is a null model with only random effects). However, diagnostics on residuals look much worse.
cond.errbar <- cond.posthoc$emmeans %>% tidy()
#jitter plot with mean ± CI
cond.bar <- ggplot(licor.plotdata, aes(x = Treatment, y = Cond)) +
geom_jitter(alpha = 0.5, width = 0.2, height = 0) +
geom_errorbar(aes(ymin = exp(conf.low), ymax = exp(conf.high),x = Treatment, y = exp(estimate)),
width = 0.2,
data = cond.errbar) +
geom_col(aes(x = Treatment, y = exp(estimate)), alpha = 0.6,data = cond.errbar) +
geom_text(aes(y = max.c, label = c("a", "a", "a")), vjust = -1, data = letters_data) +
geom_label(aes(x = 3, y = 0.15, label = paste0("p = ", Anova(cond.lmer)$`Pr(>Chisq)` %>% round(3))), data = data.frame()) +
ylim(0, 0.15) +
xlab("% Ambient Rainfall") +
ylab(TeX("Stomatal conductance $($mol $H_{2}O$ $m^{-2}$ $s^{-1})$")) +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
Ignoring unknown aesthetics: y
fig1plots <- plot_grid(photo.bar, cond.bar, spad.bar,
ncol = 2,
nrow = 2,
align = "v",
axis = "l",
labels = "AUTO",
label_x = 0)
save_plot(here("figs", "photoplots.png"), fig1plots, ncol = 2, nrow = 2, base_width = 3)