Week 12 GLMMs
Load libraries and data
library(here)
library(dplyr)
library(lme4)
library(sjPlot)
library(ggplot2)
library(tidyr)
library(patchwork)
library(MuMIn)
library(lmerTest)
library(effects)
library(RVAideMemoire)
ohia <- read.csv(here("Week_12", "Data", "ohia_merged.csv"))We will focus on data from the greenhouse experiment. Seedlings were planted in germination trays and watered normally for four weeks. After four weeks seedlings were placed into one of two treatments: a control treatment with normal daily watering, or a ‘press drought’ treatment in which seedlings were watered every other day, to 80% capacity. After six weeks of the press drought half of the seedlings (randomly chosen at the start of the experiment) were harvested for phenotypic measurements (if they survived), while the remaining seedlings (from both treatments) were placed into terminal drought (no more watering). The attached file contains data on each seedling in the experiment. The column definitions are as follows:
ID Seedling experimental identification number. BLOCK Experimental block, due to variation in germination timing. POPULATION Source of seeds, collected on Oahu Island. DROUGHT Experimental water treatment: control (watered every other day to 100% field capacity) and drought (watered every other day to 80% field capacity). FATE Harvested after press drought treatment, or subjected to terminal drought with no water until the seedling died. DEAD Binary mortality (1 = dead; 0 = alive). HEIGHT Seedling height at harvest (cm). SHOOT Dry shoot biomass at harvest (g). ROOT Dry root biomass at harvest (g). LONGEVITY Number of days in terminal drought with no water before dying. TERMINAL_GWC Soil gravimetric water content (% field capacity) at the time of mortality in terminal drought. MAR Historical mean annual rainfall at the site from which the seeds were collected (cm) MET Historical mean evapotranspiration at the site from which the seeds were collected (mm)
Pop Drought & Height
- Using the seedlings that were harvested after the press drought treatment, analyze how seedling height varies among treatments and populations. Specifically, set up a model to test whether populations differ in height, whether the press drought treatment affects height, and whether populations differ in how press drought affects height.
# harvested after drought treatment - need to be living & harvest fate
harvestohia <- ohia %>%
filter(Dead == 0,
Fate == "harvest",
!Population %in% c("Palikea", "Waimano")) %>% #these two pops have 0 in drought
mutate(Population = as.factor(Population),
Drought = as.factor(Drought),
Block = as.factor(Block))
ggplot(harvestohia, aes(Drought, Height)) + geom_boxplot() + facet_wrap(~Population)
## Models
mod1 = lmer(Height ~ Population*Drought + (1|Block), data = harvestohia)
mod2 = lmer(Height ~ MAR*Drought + (1|Population) + (1|Block), data = harvestohia)
mod3 = lmer(Height ~ MET*Drought + (1|Population) + (1|Block), data = harvestohia)Summary of each mod
Height ~ Population*Drought + (1|Block)
summary(mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Height ~ Population * Drought + (1 | Block)
## Data: harvestohia
##
## REML criterion at convergence: 421.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4117 -0.5495 -0.1005 0.4236 4.1822
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.2107 0.459
## Residual 1.1923 1.092
## Number of obs: 153, groups: Block, 9
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.03897 0.62218 7.40139 1.670
## PopulationKuaokala 1.23739 0.80753 9.29240 1.532
## PopulationKuliouou 0.62103 0.91448 6.91632 0.679
## PopulationManana1 1.29791 0.80470 10.08769 1.613
## PopulationManana2 -0.22956 0.95789 16.65682 -0.240
## PopulationMaumae 0.52810 0.79569 10.04002 0.664
## PopulationMokuleia3 1.28711 1.05653 23.62711 1.218
## PopulationPauoa1 -0.11891 0.71888 122.71734 -0.165
## PopulationPoamoho2 1.07000 0.76112 9.24936 1.406
## PopulationPoamoho3 0.89512 0.75618 8.87835 1.184
## PopulationPTAA 0.82619 0.53280 122.58156 1.551
## PopulationPTAB 0.79853 0.86421 5.54023 0.924
## PopulationWaahila1 -0.27313 0.95077 18.16735 -0.287
## PopulationWaahila2 2.52687 0.95077 18.16735 2.658
## Droughtdrought 0.28333 0.89155 122.58156 0.318
## PopulationKuaokala:Droughtdrought 0.94972 1.20044 121.75711 0.791
## PopulationKuliouou:Droughtdrought -1.44333 1.49185 122.58156 -0.967
## PopulationManana1:Droughtdrought -0.04583 1.11444 122.58156 -0.041
## PopulationManana2:Droughtdrought 0.01667 1.26084 122.58156 0.013
## PopulationMaumae:Droughtdrought -0.27772 1.08773 122.65513 -0.255
## PopulationMokuleia3:Droughtdrought -0.75333 1.27651 122.58156 -0.590
## PopulationPauoa1:Droughtdrought -0.68109 1.30732 124.59874 -0.521
## PopulationPoamoho2:Droughtdrought 0.49305 1.05607 123.09481 0.467
## PopulationPoamoho3:Droughtdrought -1.04135 1.17296 123.46102 -0.888
## PopulationPTAA:Droughtdrought -0.25801 0.99419 122.58156 -0.260
## PopulationPTAB:Droughtdrought -0.58750 1.03749 122.58156 -0.566
## PopulationWaahila1:Droughtdrought -0.38417 1.27271 124.96141 -0.302
## PopulationWaahila2:Droughtdrought -3.28333 1.22081 122.58156 -2.689
## Pr(>|t|)
## (Intercept) 0.13654
## PopulationKuaokala 0.15875
## PopulationKuliouou 0.51914
## PopulationManana1 0.13757
## PopulationManana2 0.81352
## PopulationMaumae 0.52183
## PopulationMokuleia3 0.23516
## PopulationPauoa1 0.86889
## PopulationPoamoho2 0.19247
## PopulationPoamoho3 0.26724
## PopulationPTAA 0.12356
## PopulationPTAB 0.39395
## PopulationWaahila1 0.77715
## PopulationWaahila2 0.01594 *
## Droughtdrought 0.75118
## PopulationKuaokala:Droughtdrought 0.43040
## PopulationKuliouou:Droughtdrought 0.33521
## PopulationManana1:Droughtdrought 0.96726
## PopulationManana2:Droughtdrought 0.98947
## PopulationMaumae:Droughtdrought 0.79890
## PopulationMokuleia3:Droughtdrought 0.55617
## PopulationPauoa1:Droughtdrought 0.60330
## PopulationPoamoho2:Droughtdrought 0.64141
## PopulationPoamoho3:Droughtdrought 0.37638
## PopulationPTAA:Droughtdrought 0.79567
## PopulationPTAB:Droughtdrought 0.57224
## PopulationWaahila1:Droughtdrought 0.76327
## PopulationWaahila2:Droughtdrought 0.00815 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itHeight ~ MAR*Drought + (1|Population) + (1|Block)
summary(mod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Height ~ MAR * Drought + (1 | Population) + (1 | Block)
## Data: harvestohia
##
## REML criterion at convergence: 510.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6253 -0.6370 -0.1249 0.4152 5.0893
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2021 0.4496
## Block (Intercept) 0.1523 0.3902
## Residual 1.2620 1.1234
## Number of obs: 153, groups: Population, 14; Block, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.946e+00 3.932e-01 1.215e+01 4.948 0.000325 ***
## MAR -8.985e-04 1.409e-03 1.299e+01 -0.638 0.534873
## Droughtdrought 3.076e-02 3.407e-01 1.393e+02 0.090 0.928188
## MAR:Droughtdrought -6.055e-04 1.330e-03 1.404e+02 -0.455 0.649591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MAR Drghtd
## MAR -0.813
## Droghtdrght -0.350 0.316
## MAR:Drghtdr 0.280 -0.360 -0.826Height ~ MET*Drought + (1|Population) + (1|Block)
summary(mod3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Height ~ MET * Drought + (1 | Population) + (1 | Block)
## Data: harvestohia
##
## REML criterion at convergence: 513
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5552 -0.5951 -0.1426 0.3954 5.1915
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.17938 0.4235
## Block (Intercept) 0.09965 0.3157
## Residual 1.27600 1.1296
## Number of obs: 153, groups: Population, 14; Block, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.339e+00 5.719e-01 1.163e+01 4.089 0.0016 **
## MET -6.700e-04 6.133e-04 1.226e+01 -1.092 0.2956
## Droughtdrought -2.925e-01 5.227e-01 1.344e+02 -0.560 0.5767
## MET:Droughtdrought 2.456e-04 5.969e-04 1.342e+02 0.411 0.6814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MET Drghtd
## MET -0.931
## Droghtdrght -0.400 0.383
## MET:Drghtdr 0.362 -0.394 -0.930Anovas
Height ~ Population*Drought + (1|Block)
anova(mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Population 40.890 3.1454 13 11.70 2.6381 0.0532 .
## Drought 1.439 1.4394 1 124.98 1.2072 0.2740
## Population:Drought 23.901 1.8386 13 122.68 1.5420 0.1117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The press drought treatment does not impact height, and populations do not appear to differ in how press drought impacts height given the p-values of 0.11 and 0.27. There is a small evidence that height varies among population given p-val = 0.05, but all are insignificant.
Height ~ MAR*Drought + (1|Population) + (1|Block)
anova(mod2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## MAR 1.03790 1.03790 1 10.497 0.8224 0.3848
## Drought 0.01029 0.01029 1 139.346 0.0082 0.9282
## MAR:Drought 0.26162 0.26162 1 140.366 0.2073 0.6496Height ~ MET*Drought + (1|Population) + (1|Block)
anova(mod3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## MET 1.18985 1.18985 1 8.943 0.9325 0.3596
## Drought 0.39951 0.39951 1 134.442 0.3131 0.5767
## MET:Drought 0.21601 0.21601 1 134.224 0.1693 0.6814Plots
Height ~ Population*Drought + (1|Block)
plot_model(mod1, type = 're')
Height ~ MAR*Drought + (1|Population) + (1|Block)
plot_model(mod2, type = 're')## [[1]]
##
## [[2]]
Height ~ MET*Drought + (1|Population) + (1|Block)
plot_model(mod3, type = 're')## [[1]]
##
## [[2]]
GLMMs Pop & Mortality
- Using the seedlings that were fated to be harvested after the press drought treatment, analyze how survival during the press drought is affected by the various potential predictors. Specifically, set up a model to test whether populations vary in probability of mortality (i.e., the response column DEAD), whether the press drought treatment affects mean mortality, and whether populations vary in how mortality is affected by the press drought treatment.
Next, set up a model to test whether historical rainfall at a site predicts mortality, and whether historical rainfall at a site predicts the effect of press drought on mortality. Finally, set up a model to test whether historical evapotranspiration at a site predicts mortality, and whether historical evapotranspiration at a site predicts the effect of press drought on mortality. Be sure to think about what the appropriate probability distribution is for the response variable. How do you interpret the results?
Models
ohiamortality <- ohia %>%
filter(Fate == "harvest")
gmod1= glmer(Dead ~ Population*Drought + (1|Block), data = ohiamortality, family = binomial)
gmod2 = glmer(Dead ~ MAR*Drought + (1|Population) + (1|Block), data = ohiamortality, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?gmod3 = glmer(Dead ~ MET*Drought + (1|Population) + (1|Block), data = ohiamortality, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## See ?lme4::convergence and ?lme4::troubleshooting.Summary of each mod
Dead ~ Population*Drought + (1|Block)
summary(gmod1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Dead ~ Population * Drought + (1 | Block)
## Data: ohiamortality
##
## AIC BIC logLik -2*log(L) df.resid
## 535.1 670.7 -234.5 469.1 418
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2091 -0.6655 0.2718 0.5879 2.3327
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.967 0.9834
## Number of obs: 451, groups: Block, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.93190 1.14924 -0.811 0.41743
## PopulationKuaokala 1.97864 1.43008 1.384 0.16648
## PopulationKuliouou 1.14368 1.55938 0.733 0.46330
## PopulationManana1 1.08991 1.40317 0.777 0.43731
## PopulationManana2 3.99141 1.54160 2.589 0.00962 **
## PopulationMaumae 1.31834 1.35238 0.975 0.32965
## PopulationMokuleia3 3.62883 1.66244 2.183 0.02905 *
## PopulationPalikea 0.79511 1.68063 0.473 0.63614
## PopulationPauoa1 1.37825 1.07582 1.281 0.20015
## PopulationPoamoho2 1.76207 1.27935 1.377 0.16841
## PopulationPoamoho3 1.23487 1.37825 0.896 0.37027
## PopulationPTAA -0.16331 0.99846 -0.164 0.87008
## PopulationPTAB -0.08696 1.66056 -0.052 0.95823
## PopulationWaahila1 1.41064 1.35129 1.044 0.29652
## PopulationWaahila2 1.15402 1.48490 0.777 0.43706
## PopulationWaimano 1.79580 1.85813 0.966 0.33381
## Droughtdrought 2.22616 1.16249 1.915 0.05549 .
## PopulationKuaokala:Droughtdrought -1.14636 1.45864 -0.786 0.43192
## PopulationKuliouou:Droughtdrought 0.18883 1.65259 0.114 0.90903
## PopulationManana1:Droughtdrought -1.05260 1.40849 -0.747 0.45487
## PopulationManana2:Droughtdrought -1.91425 1.48887 -1.286 0.19855
## PopulationMaumae:Droughtdrought -1.11802 1.38086 -0.810 0.41814
## PopulationMokuleia3:Droughtdrought -3.45341 1.63037 -2.118 0.03416 *
## PopulationPalikea:Droughtdrought 12.78258 296.22395 0.043 0.96558
## PopulationPauoa1:Droughtdrought -0.81855 1.53911 -0.532 0.59484
## PopulationPoamoho2:Droughtdrought -1.59513 1.31535 -1.213 0.22524
## PopulationPoamoho3:Droughtdrought 0.76598 1.44012 0.532 0.59480
## PopulationPTAA:Droughtdrought -1.75383 1.40363 -1.249 0.21148
## PopulationPTAB:Droughtdrought -2.05073 1.47952 -1.386 0.16572
## PopulationWaahila1:Droughtdrought -1.84747 1.52274 -1.213 0.22503
## PopulationWaahila2:Droughtdrought -2.51259 1.61941 -1.552 0.12077
## PopulationWaimano:Droughtdrought 13.62222 229.45229 0.059 0.95266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 32 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itDead ~ MAR*Drought + (1|Population) + (1|Block)
summary(gmod2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Dead ~ MAR * Drought + (1 | Population) + (1 | Block)
## Data: ohiamortality
##
## AIC BIC logLik -2*log(L) df.resid
## 516.4 541.0 -252.2 504.4 445
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3968 -0.8109 0.3756 0.5843 1.7398
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.001837 0.04286
## Block (Intercept) 1.042499 1.02103
## Number of obs: 451, groups: Population, 16; Block, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.343131 0.506110 -0.678 0.4978
## MAR 0.003278 0.001339 2.448 0.0144 *
## Droughtdrought 0.607037 0.462166 1.313 0.1890
## MAR:Droughtdrought 0.001629 0.001636 0.995 0.3196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MAR Drghtd
## MAR -0.649
## Droghtdrght -0.448 0.523
## MAR:Drghtdr 0.383 -0.562 -0.875
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?MAR appears to be a significant predictor of mortality given the p-value of 0.014.
Dead ~ MET*Drought + (1|Population) + (1|Block)
summary(gmod3)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Dead ~ MET * Drought + (1 | Population) + (1 | Block)
## Data: ohiamortality
##
## AIC BIC logLik -2*log(L) df.resid
## 519.6 544.3 -253.8 507.6 445
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6142 -0.7679 0.3759 0.5790 1.6553
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.1484 0.3852
## Block (Intercept) 1.0214 1.0106
## Number of obs: 451, groups: Population, 16; Block, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1494229 0.7989206 -0.187 0.852
## MET 0.0006951 0.0007943 0.875 0.382
## Droughtdrought -0.2936294 0.7716470 -0.381 0.704
## MET:Droughtdrought 0.0014625 0.0008297 1.763 0.078 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MET Drghtd
## MET -0.870
## Droghtdrght -0.474 0.486
## MET:Drghtdr 0.451 -0.500 -0.957
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## See ?lme4::convergence and ?lme4::troubleshooting.Plots
Dead ~ Population*Drought + (1|Block)
plot(allEffects(gmod1))
Dead ~ MAR*Drought + (1|Population) + (1|Block)
plot(allEffects(gmod2))
Dead ~ MET*Drought + (1|Population) + (1|Block)
plot(allEffects(gmod3))
Random effects gmod3
plot_model(gmod3, type ='re')## [[1]]
##
## [[2]]
Anovas
Dead ~ Population*Drought + (1|Block)
anova(gmod1)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Population 15 23.418 1.5612 1.5612
## Drought 1 15.065 15.0652 15.0652
## Population:Drought 15 15.949 1.0633 1.0633The press drought treatment does not impact height, and populations do not appear to differ in how press drought impacts height given the p-values of 0.11 and 0.27. There is a small evidence that height varies among population given p-val = 0.05, but all are insignificant.
Dead ~ MAR*Drought + (1|Population) + (1|Block)
anova(gmod2)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## MAR 1 11.2469 11.2469 11.2469
## Drought 1 20.4295 20.4295 20.4295
## MAR:Drought 1 0.9959 0.9959 0.9959Dead ~ MET*Drought + (1|Population) + (1|Block)
anova(gmod3)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## MET 1 3.5991 3.5991 3.5991
## Drought 1 20.5293 20.5293 20.5293
## MET:Drought 1 3.1069 3.1069 3.1069AIC
AIC(gmod1, gmod2, gmod3)## df AIC
## gmod1 33 535.0576
## gmod2 6 516.3733
## gmod3 6 519.6130It appears that MAR, Historical mean annual rainfall at the site from which the seeds were collected , is the best predictor for seedling mortality given the lowest AIC value of 516.