BACI analysis for taxonomic richness (N0)
Valentin Stefan
21 Feb 2021
1 Load packages and helper functions
# Source scripts that contains various helper functions for reading and plotting
# data.
source(file = "analysis/baci_glmm/helper-functions.r")
load_install_packages()## glmmTMB AICcmodavg lsmeans DHARMa parallel data.table
## TRUE TRUE TRUE TRUE TRUE TRUE
## magrittr ggplot2 directlabels
## TRUE TRUE TRUE2 Prepare date
# Read data
baci_dt <- fread("data/variables_for_BACI.csv", stringsAsFactors = TRUE)
# Rename some columns:
setnames(baci_dt,
old = c("BACItime", "BACIpos", "Site", "Ind./m2", "% EPT", "Total organisms", "Total EPT organisms"),
new = c("period_ba", "treatment_ci", "site_f", "ind_m2", "ept_prc", "N", "N_EPT"))
# Create a factor column `year_f` where I rounded the year to integer and then
# converted to factor type.
baci_dt[, year_f := Year %>% round(digits = 0) %>% factor]
# This will be needed for testing for overdispersation
baci_dt[, obs := 1:.N]
# Variable of interest
varb <- "N0"3 Exploratory plots
3.1 Histogram of the explained variable
hist(baci_dt[[varb]],
main = paste("Histogram of" , varb),
xlab = varb)
3.2 Preliminary plot
preliminary_plot <- exploratory_plot(varb)
preliminary_plot
Fig. 1 - Observed values and trends. Thicker lines represent mean lines. Thinner lines stand for each site. There is only one control site which actually overlaps with the mean line for control.
We can observe already that there is no clear difference between “impact” and “control”. The two mean lines are mostly parallel.
4 GLMM fitting
The taxonomic richness (N0) represents counts. We can use the Poisson error distribution.
4.1 Fixed and random effects
Specify fixed and random effects and test for random effects structure.
The general model formula with fixed effects is varb ~ period_ba + treatment_ci + period_ba:treatment_ci Note the interaction term period_ba:treatment_ci which is the “BACI effect” (see Schwarz, 2015). Testing for its statistical significance is equivalent to testing for an environmental impact (Schwarz, 2015).
To the fixed effects model, random effects are added: site_f and year_f.
Below, it was tested if site_f:year_f interaction should be kept in the model.
# Create an empty list to be populated with models
models <- list()
# Fit model without the interaction of random effects (site_f and year_f)
models[[1]] <- glmmTMB(N0 ~ period_ba * treatment_ci + (1|site_f) + (1|year_f),
data = baci_dt, family = poisson(link = "log"))
# Model with interaction in the random effects structure
models[[2]] <- glmmTMB(N0 ~ period_ba * treatment_ci + (1|site_f) + (1|year_f) + (1|site_f:year_f),
data = baci_dt, family = poisson(link = "log"))Note that the warnings 'giveCsparse' has been deprecated; setting 'repr = "T"' for you are expected - see comment of Ben Bolker here.
Model summary:
summary(models[[1]])## Family: poisson ( log )
## Formula: N0 ~ period_ba * treatment_ci + (1 | site_f) + (1 | year_f)
## Data: baci_dt
##
## AIC BIC logLik deviance df.resid
## 957.3 976.5 -472.7 945.3 174
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site_f (Intercept) 0.06449 0.25395
## year_f (Intercept) 0.00209 0.04572
## Number of obs: 180, groups: site_f, 6; year_f, 5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.10736 0.26471 7.961 1.71e-15 ***
## period_babefore -0.05136 0.17014 -0.302 0.763
## treatment_ciimpact 0.31192 0.28845 1.081 0.280
## period_babefore:treatment_ciimpact 0.02585 0.17312 0.149 0.881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(models[[2]])## Family: poisson ( log )
## Formula:
## N0 ~ period_ba * treatment_ci + (1 | site_f) + (1 | year_f) +
## (1 | site_f:year_f)
## Data: baci_dt
##
## AIC BIC logLik deviance df.resid
## 959.3 981.7 -472.7 945.3 173
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site_f (Intercept) 6.449e-02 2.539e-01
## year_f (Intercept) 2.090e-03 4.572e-02
## site_f:year_f (Intercept) 5.396e-10 2.323e-05
## Number of obs: 180, groups: site_f, 6; year_f, 5; site_f:year_f, 30
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.10736 0.26471 7.961 1.71e-15 ***
## period_babefore -0.05136 0.17014 -0.302 0.763
## treatment_ciimpact 0.31192 0.28845 1.081 0.280
## period_babefore:treatment_ciimpact 0.02585 0.17312 0.149 0.881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The models are compared using the Akaike Information Criterion (AIC). The preferred model should be the one with a smaller AIC value. A model selection approach similar to the one proposed by Burnham & Anderson (2002) is used. The authors suggest the rule of thumb according to which models with AIC difference smaller or equal than 2 (delta-AIC ≤ 2) have substantial empirical support, those with 4 ≤ delta-AIC ≤ 7 have considerably less and those with delta-AIC > 10 have essentially none. On the other hand, they also suggest that the model with delta-AIC > 10 might still be used for inference if the sample size is large.
# AICc comparison
aictab(cand.set = models,
modnames = c("1-without site_f:year_f",
"2-with site_f:year_f"),
second.ord = TRUE)##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## 1-without site_f:year_f 6 957.80 0.00 0.75 0.75 -472.65
## 2-with site_f:year_f 7 959.96 2.17 0.25 1.00 -472.65Additionally to the AIC test, we can also run the likelihood ratio test (LRT) as per Bates (2015).
anova(models[[1]], models[[2]])## Data: baci_dt
## Models:
## models[[1]]: N0 ~ period_ba * treatment_ci + (1 | site_f) + (1 | year_f), zi=~0, disp=~1
## models[[2]]: N0 ~ period_ba * treatment_ci + (1 | site_f) + (1 | year_f) + , zi=~0, disp=~1
## models[[2]]: (1 | site_f:year_f), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## models[[1]] 6 957.31 976.47 -472.65 945.31
## models[[2]] 7 959.31 981.66 -472.65 945.31 0 1 1The AIC test favors that the less complex model (without site_f:year_f). The LRT test indicates that there are no significant differences between the two models. According to the rule of parsimony, we pick the simpler model for further analysis.
4.2 Overdispersion
Testing for overdispersion in the mixed model.
“Overdispersion: the occurrence of more variance in the data than predicted by a statistical model.” (Bolker, 2009). In case of overdispersed then we can fit a model with observation-level random effects see Harrison XA (2014) or Harrison XA (2015). This approach was also suggested by Ben Bolker here, where further references are mentioned.
The tests below indicate no significant overdispersion.
# Test for overdispersion with DHARMa residuals diagnostics
simulation <- DHARMa::simulateResiduals(fittedModel = models[[1]])
plot(simulation)
# A non-parametric test on the simulated residuals
DHARMa::testDispersion(simulation)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.0828, p-value = 0.552
## alternative hypothesis: two.sidedWe keep the model without any overdispersion adjustments for further analysis.
model <- models[[1]]4.3 Diagnostic plots
# Residuals
qqnorm(residuals(model), main = "Q-Q plot - residuals")
qqline(residuals(model), col = "red")
# inspecting the random effects (see also Bolker, 2009 - supp 1)
qqnorm(unlist(ranef(model)), main = "Q-Q plot, random effects")
qqline(unlist(ranef(model)), col = "red")
# fitted vs residuals
scatter.smooth(fitted(model),
residuals(model, type = "pearson"),
main = "fitted vs residuals",
xlab = "Fitted Values",
ylab = "Residuals")
abline(h = 0, col = "red")
# fitted vs observed
scatter.smooth(fitted(model), baci_dt[[varb]],
xlab = "Fitted Values",
ylab = "Observed")
abline(0, 1, col = "red")
4.4 Significance of BACI
We will test for statistical significance of BACI interaction term using the likelihood ratio test (LRT).
The nested models that we compare are:
- the model without the interaction between the fixed effects (the BACI effect)
- model with the BACI effect
# Model without the BACI interaction term
model.no.interaction <- glmmTMB(N0 ~ period_ba + treatment_ci + (1|site_f) + (1|year_f),
data = baci_dt, family = poisson(link = "log"))
# LRT
lrt <- anova(model, model.no.interaction)
lrt## Data: baci_dt
## Models:
## model.no.interaction: N0 ~ period_ba + treatment_ci + (1 | site_f) + (1 | year_f), zi=~0, disp=~1
## model: N0 ~ period_ba * treatment_ci + (1 | site_f) + (1 | year_f), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model.no.interaction 5 955.33 971.30 -472.67 945.33
## model 6 957.31 976.47 -472.65 945.31 0.0224 1 0.881The likelihood ratio (LRT) gives no significant p-value. There are no significant differences between the two models. Therefore, the interaction term (BACI effect) is not statistically significant (there is no environmental impact).
4.5 Model coefficients
Extract coefficients from final model.
Remove the intercept for extracting group means and their confidence intervals as suggested in Schielzeth H (2010).
final.model <- model
final.model.noIntercept <- update(final.model, . ~ . -1)To get the real proportions (Least-squares means/predicted marginal means/treatment means) one can do:
estimates <- lsmeans::lsmeans(final.model.noIntercept, ~ treatment_ci:period_ba, type = "response")
# https://stats.stackexchange.com/questions/192062/issue-calculating-adjusted-means-for-glmer-model
# Confidence level used: 0.95
estimates## treatment_ci period_ba rate SE df lower.CL upper.CL
## control after 8.23 2.18 174 4.88 13.9
## impact after 11.24 1.34 174 8.88 14.2
## control before 7.81 2.32 174 4.35 14.0
## impact before 10.95 1.47 174 8.41 14.3
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleAs indicated in Schwarz CJ (2015), the BACI effect is computed as BACI = avgCA - avgCB - (avgIA - avgIB):
# get the estimates only (without CI)
est <- predict(ref.grid(final.model.noIntercept), type = "response")
# give names to the estimates vector
names(est) <- c("CA","CB","IA","IB"); est ## CA CB IA IB
## 8.226456 7.814639 11.237692 10.954731baci <- est["CA"]-est["CB"]-(est["IA"]-est["IB"])
baci## CA
## 0.128856One can also get the BACI effect like:
contrast(regrid(estimates), list(baci=c(1,-1,-1,1)))## contrast estimate SE df t.ratio p.value
## baci 0.129 1.46 174 0.089 0.9296Or with asymptotic CI-s:
baci_ci <- confint(contrast(regrid(estimates), list(baci=c(1,-1,-1,1))))
baci_ci## contrast estimate SE df lower.CL upper.CL
## baci 0.129 1.46 174 -2.74 3
##
## Confidence level used: 0.95# Confidence level used: 0.95
# https://stats.stackexchange.com/questions/241523/testing-for-pairwise-proportion-differences4.6 Coefficient of determination (R^2)
Below is presented the coefficient of determination of the correlation between the fitted and observed values as suggested by Ben Bolker here.
r2_corr <- function(m) {
lmfit <- lm(model.response(model.frame(m)) ~ fitted(m))
summary(lmfit)$r.squared
}
R2 <- r2_corr(final.model)
R2## [1] 0.52938184.7 Interaction plot
This plot will confirm as well visually that there is no BACI effect. The confidence intervals of the means overlap and the two lines are almost parallel.
interaction_plot(estimates, varb)
Fig. 2 - Interaction plot for the model. Least square means values. Error bars show the ±95% CI.
We can compare the observed mean values (marked with x symbol in the graph below) with the estimated mean values (model coefficients, marked with a plain dot and connected with lines). The estimated means for N0 are close to the observed mean values.
interaction_plot(estimates, varb) +
# Add observed means as dots with bootstrapped CIs
stat_summary(data = baci_dt,
aes(x = period_ba,
y = get(varb),
group = treatment_ci,
color = treatment_ci),
fun.data = "mean_cl_boot",
# geom = "point",
size = 0.5,
shape = 4,
position = position_dodge(width = 0.1),
show.legend = TRUE)
Fig. 3 - Interaction plot depicting both the observed mean values (marked with x) and the estimated mean values (model coefficents, marked with plain dot). Error bars show the ±95% CI.
4.8 Model summary
The final GLMM model explained 52.94 % of the variance in N0.
There was no significant BACI period × treatment effect (LRT, p = 0.881). The BACI effect estimated from the model, which is the difference of the two changes (control after − control before) − (impact after − impact before]), was 0.13 ± 1.46 standard error.
The estimated mean N0 in the control sites varied from 7.81 to 8.23, and in the impact sites varied from 10.95 to 11.24, before and after the construction of the dam (Fig 2 & Fig 3).
5 References
Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious mixed models. arXiv preprint arXiv:1506.04967.
Bolker BM et al. (2009) Generalized linear mixed models: a practical guide for ecology and evolution. Trends in ecology & evolution 24:127–135 at http://www.sciencedirect.com/science/article/pii/S0169534709000196
Burnham, K. & Anderson, D., 2002. Model selection and multimodel inference: a practical information-theoretic approach, New York: Springer.
Harrison XA (2014) Using observation-level random effects to model overdispersion in count data in ecology and evolution. PeerJ 2:e616 https://doi.org/10.7717/peerj.616
Harrison XA (2015) A comparison of observation-level random effect and Beta-Binomial models for modelling overdispersion in Binomial data in ecology & evolution. PeerJ 3:e1114 at https://peerj.com/articles/1114.pdf
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of Memory and Language, 94, 305-315.
Schwarz CJ (2015) Analysis of BACI experiments. In Course Notes for Beginning and Intermediate Statistics. at http://people.stat.sfu.ca/~cschwarz/Stat-650/Notes/PDFbigbook-R/R-part013.pdf
Schielzeth H (2010) Simple means to improve the interpretability of regression coefficients. Methods in Ecology and Evolution, 1(2), 103-113.