Sarah Pennington
For this assignment we will focus on springtail abundance and how
these are affected by temperature and mite predation.
In the attached dataset we will use the columns total.prey (total
Collembola count), temp.factor (the temperature treatment), and
predators (the predator treatment).
1. First use an appropriate ‘standard’ model for counts (i.e., not
zero-inflated) to ask whether the abundance of Collembola is affected by
temperature, by the presence of one or more predators, and whether the
effect of predators depends on temperature. Use the function glmmTMB()
to fit the model, because later we will compare this model to others fit
with this function.
col <- read_csv("collembola_revised_2025.csv")
## New names:
## Rows: 60 Columns: 16
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (2): predators, temp.factor dbl (14): ...1, Pred1, Pred2, TEMP, Litter_loss,
## Cmic, Fol_pop, Pro_pop, Pre...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
str(col)
## spc_tbl_ [60 × 16] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ...1 : num [1:60] 1 2 3 4 5 6 7 8 9 10 ...
## $ predators : chr [1:60] "none" "none" "none" "none" ...
## $ Pred1 : num [1:60] 0 0 0 0 0 0 0 0 0 0 ...
## $ Pred2 : num [1:60] 0 0 0 0 0 0 0 0 0 0 ...
## $ TEMP : num [1:60] 0 0 0 0 0 1 1 1 1 1 ...
## $ Litter_loss : num [1:60] 0.416 0.448 0.637 0.454 0.466 ...
## $ Cmic : num [1:60] 28658 21208 17997 8766 13371 ...
## $ Fol_pop : num [1:60] 225 382 0 415 445 ...
## $ Pro_pop : num [1:60] 163 0 33 263 0 0 150 0 0 27 ...
## $ Pred_1_pop : num [1:60] NA NA NA NA NA NA NA NA NA NA ...
## $ Pred_2_pop : num [1:60] NA NA NA NA NA NA NA NA NA NA ...
## $ lipid_protein_Folsomia : num [1:60] 0.425 0.556 NA 0.378 0.472 ...
## $ bodylength_mean_Folsomia: num [1:60] 1301 1763 NA 1492 1632 ...
## $ total.prey : num [1:60] 388 382 33 678 445 ...
## $ temp : num [1:60] 0 0 0 0 0 1 1 1 1 1 ...
## $ temp.factor : chr [1:60] "12_15" "12_15" "12_15" "12_15" ...
## - attr(*, "spec")=
## .. cols(
## .. ...1 = col_double(),
## .. predators = col_character(),
## .. Pred1 = col_double(),
## .. Pred2 = col_double(),
## .. TEMP = col_double(),
## .. Litter_loss = col_double(),
## .. Cmic = col_double(),
## .. Fol_pop = col_double(),
## .. Pro_pop = col_double(),
## .. Pred_1_pop = col_double(),
## .. Pred_2_pop = col_double(),
## .. lipid_protein_Folsomia = col_double(),
## .. bodylength_mean_Folsomia = col_double(),
## .. total.prey = col_double(),
## .. temp = col_double(),
## .. temp.factor = col_character()
## .. )
## - attr(*, "problems")=<externalptr>
col$temp.factor <- as.factor(col$temp.factor)
col$predators <- as.factor(col$predators)
#First use an appropriate ‘standard’ model for counts (i.e., not zero-inflated) to ask whether the abundance of Collembola is affected by temperature, by the presence of one or more predators, and whether the effect of predators depends on temperature.
plot(col$total.prey ~ col$temp.factor, ylab = 'total prey', xlab = 'temp')
plot(col$total.prey ~ col$predators, ylab = 'total prey', xlab = 'predators')
#Overdispersed?
model <- glm(total.prey ~ temp.factor * predators, data = col, family= poisson)
summary(model)
##
## Call:
## glm(formula = total.prey ~ temp.factor * predators, family = poisson,
## data = col)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.69441 0.02594 219.512 < 2e-16 ***
## temp.factor17_20 0.61442 0.03220 19.080 < 2e-16 ***
## temp.factor22_25 0.04023 0.03632 1.108 0.268
## predatorsHA+HM 0.63887 0.03207 19.924 < 2e-16 ***
## predatorsHM 0.33291 0.03399 9.794 < 2e-16 ***
## predatorsnone 0.25936 0.03453 7.512 5.84e-14 ***
## temp.factor17_20:predatorsHA+HM -0.91079 0.04324 -21.062 < 2e-16 ***
## temp.factor22_25:predatorsHA+HM -0.63179 0.04813 -13.127 < 2e-16 ***
## temp.factor17_20:predatorsHM -1.08955 0.04791 -22.743 < 2e-16 ***
## temp.factor22_25:predatorsHM -1.67557 0.06540 -25.621 < 2e-16 ***
## temp.factor17_20:predatorsnone 0.69416 0.04119 16.853 < 2e-16 ***
## temp.factor22_25:predatorsnone 0.68510 0.04572 14.986 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 45181 on 59 degrees of freedom
## Residual deviance: 33693 on 48 degrees of freedom
## AIC: 34028
##
## Number of Fisher Scoring iterations: 7
dispersion <- sum(residuals(model, type = "pearson")^2)/(length(model$y)- length(model$coefficients))
print(dispersion)
## [1] 641.9313
#Yep, overdispersed. For 'standard' model
model <- glm.nb(total.prey ~ temp.factor * predators, data = col)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
#
summary(model)
##
## Call:
## glm.nb(formula = total.prey ~ temp.factor * predators, data = col,
## init.theta = 2512782.316, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.69441 0.02594 219.499 < 2e-16 ***
## temp.factor17_20 0.61442 0.03220 19.079 < 2e-16 ***
## temp.factor22_25 0.04023 0.03633 1.107 0.268
## predatorsHA+HM 0.63887 0.03207 19.923 < 2e-16 ***
## predatorsHM 0.33291 0.03399 9.794 < 2e-16 ***
## predatorsnone 0.25936 0.03453 7.511 5.86e-14 ***
## temp.factor17_20:predatorsHA+HM -0.91079 0.04325 -21.060 < 2e-16 ***
## temp.factor22_25:predatorsHA+HM -0.63179 0.04813 -13.126 < 2e-16 ***
## temp.factor17_20:predatorsHM -1.08955 0.04791 -22.741 < 2e-16 ***
## temp.factor22_25:predatorsHM -1.67557 0.06540 -25.620 < 2e-16 ***
## temp.factor17_20:predatorsnone 0.69416 0.04119 16.852 < 2e-16 ***
## temp.factor22_25:predatorsnone 0.68510 0.04572 14.985 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2512782) family taken to be 1)
##
## Null deviance: 45171 on 59 degrees of freedom
## Residual deviance: 33685 on 48 degrees of freedom
## AIC: 34022
##
## Number of Fisher Scoring iterations: 1
## Error in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, : invalid 'nsmall' argument
#
#Use the function glmmTMB() to fit the model, because later we will compare this model to others fit with this function. (exclude zero inflated part)
#Trying to figure out the negative binomial
mod <- glmmTMB(total.prey ~ temp.factor * predators, data =col, family = "nbinom2")
#
summary(mod)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * predators
## Data: col
##
## AIC BIC logLik deviance df.resid
## 734.5 761.7 -354.2 708.5 47
##
##
## Dispersion parameter for nbinom2 family (): 0.188
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.69440 1.03150 5.520 3.38e-08 ***
## temp.factor17_20 0.61436 1.45864 0.421 0.674
## temp.factor22_25 0.04023 1.45876 0.028 0.978
## predatorsHA+HM 0.63890 1.45867 0.438 0.661
## predatorsHM 0.33289 1.45870 0.228 0.819
## predatorsnone 0.25936 1.45872 0.178 0.859
## temp.factor17_20:predatorsHA+HM -0.91071 2.06281 -0.441 0.659
## temp.factor22_25:predatorsHA+HM -0.63180 2.06293 -0.306 0.759
## temp.factor17_20:predatorsHM -1.08947 2.06290 -0.528 0.597
## temp.factor22_25:predatorsHM -1.67560 2.06338 -0.812 0.417
## temp.factor17_20:predatorsnone 0.69425 2.06276 0.337 0.736
## temp.factor22_25:predatorsnone 0.68512 2.06287 0.332 0.740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Plot fitted effects and perform likelihood ratio tests on the
relevant terms. To perform marginal tests you will want to compare pairs
of models using the function anova(). Previously we used the function
Anova() (the capitalized and uncapitalized functions are different) to
automate this process, but this function returns a less accurate test
for glmmTMB (a Wald test). So, for each term you want to test, compare a
model with this term to a model without this term.
#Plot fitted effects and perform likelihood ratio tests on the relevant terms
#trying with the which handles zero-inflated models well
effect_plot <- allEffects(mod)
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
## Error in model.frame.default(formula = total.prey ~ temp.factor * predators, : 'data' must be a data.frame, environment, or list
#Error in model.frame.default(formula = total.prey ~ temp.factor * predators, : 'data' must be a data.frame, environment, or list
#problem with col?
#TRoubleshooting
col2 <- read_excel("collembola_revised_2025_1.xlsx")
col2$temp.factor <- as.factor(col2$temp.factor)
col2$predators <- as.factor(col2$predators)
str(col2)
## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
## $ predators : Factor w/ 4 levels "HA","HA+HM","HM",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ total.prey : num [1:60] 388 382 33 678 445 ...
## $ temp.factor: Factor w/ 3 levels "12_15","17_20",..: 1 1 1 1 1 2 2 2 2 2 ...
col2 <- as.data.frame(col2)
mod.z2 <- glmmTMB(total.prey ~ temp.factor * predators, data =col, family = "nbinom1")
summary(mod.z2)
## Family: nbinom1 ( log )
## Formula: total.prey ~ temp.factor * predators
## Data: col
##
## AIC BIC logLik deviance df.resid
## 706.0 733.2 -340.0 680.0 47
##
##
## Dispersion parameter for nbinom1 family (): 1.7e+03
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2244 0.5201 11.968 <2e-16 ***
## temp.factor17_20 0.2076 0.6684 0.311 0.7561
## temp.factor22_25 -2.1635 1.1083 -1.952 0.0509 .
## predatorsHA+HM 0.9027 0.6099 1.480 0.1389
## predatorsHM 0.7172 0.6173 1.162 0.2453
## predatorsnone 0.5747 0.6228 0.923 0.3561
## temp.factor17_20:predatorsHA+HM -1.2306 0.9094 -1.353 0.1760
## temp.factor22_25:predatorsHA+HM 0.5668 1.3008 0.436 0.6630
## temp.factor17_20:predatorsHM -2.1722 1.0503 -2.068 0.0386 *
## temp.factor22_25:predatorsHM -0.7619 1.5420 -0.494 0.6213
## temp.factor17_20:predatorsnone -1.0500 0.9598 -1.094 0.2740
## temp.factor22_25:predatorsnone 0.4456 1.3708 0.325 0.7451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
effect_plotz <- allEffects(mod.z2)
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
## Error in model.frame.default(formula = total.prey ~ temp.factor * predators, : 'data' must be a data.frame, environment, or list
plot(effect_plotz)#More errors Warning: overriding variance function for
## Error: object 'effect_plotz' not found
#effects/dev.resids: computed variances may be incorrectError in
#model.frame.default(formula = total.prey ~ temp.factor * predators, : 'data' #must be a data.frame, environment, or list
#Moving on to the likelihood ratio tests
model_tempz <- glmmTMB(total.prey ~ predators, data =col, family = "nbinom2")
model_predz <- glmmTMB(total.prey ~ temp.factor, data =col, family = "nbinom2")
anova(model_tempz,mod) #chi 0.9722
## Data: col
## Models:
## model_tempz: total.prey ~ predators, zi=~0, disp=~1
## mod: total.prey ~ temp.factor * predators, zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model_tempz 5 720.73 731.2 -355.36 710.73
## mod 13 734.47 761.7 -354.24 708.47 2.2555 8 0.9722
anova(model_predz, mod) #chi 0.435
## Data: col
## Models:
## model_predz: total.prey ~ temp.factor, zi=~0, disp=~1
## mod: total.prey ~ temp.factor * predators, zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model_predz 4 720.27 728.64 -356.13 712.27
## mod 13 734.47 761.70 -354.24 708.47 3.793 9 0.9245
(PT. 2) Plot fitted effects and perform likelihood ratio tests on
the relevant terms.
#Plot fitted effects and perform likelihood ratio tests on the relevant terms
#I am just assuming we are still working with the "standard" model too?
effect_plot <- allEffects(model)
plot(effect_plot)
plot(ggeffect(model, terms = c('temp.factor', 'predators')))
model_temp <- glm.nb(total.prey ~ predators, data = col)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
#
model_pred <- glm.nb(total.prey ~ temp.factor, data = col)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
#
anova(model_temp, model) # chi 0
## Likelihood ratio tests of Negative Binomial Models
##
## Response: total.prey
## Model theta Resid. df 2 x log-lik. Test df
## 1 predators 3022113 56 -39322.31
## 2 temp.factor * predators 2512782 48 -33996.48 1 vs 2 8
## LR stat. Pr(Chi)
## 1
## 2 5325.827 0
#
anova(model_pred, model) # chi 0
## Likelihood ratio tests of Negative Binomial Models
##
## Response: total.prey
## Model theta Resid. df 2 x log-lik. Test df
## 1 temp.factor 2860757 57 -43554.53
## 2 temp.factor * predators 2512782 48 -33996.48 1 vs 2 9
## LR stat. Pr(Chi)
## 1
## 2 9558.049 0
#
How do you interpret the results at this stage?
#With the glmmTMB model, I couldn't figure out how to get the effects plot put #together through regular ggeffects. Looking at the likelihood ratio tests, #the effect of predators and tempature on total prey isn't significant for #this model.
#With the glm.nb model, I could make an effects plot. I can see in the effects #plot and in the likelihood tests that predators and temperature have a #significant impact on total prey.
2. A large proportion of the data are zeros, and it may be the case
that processes controlling abundance are different from processes
controlling ‘extra’ zeros – if, in fact, there are extra zeros.
Use glmmTMB to fit zero-inflated count model(s). You should decide
how many to fit, and which kind of probability distribution(s) to use.
Use AIC to do model selection, to determine whether incorporating zero
inflation improves model fit. Using the best model (i.e., the lowest
AIC), perform marginal likelihood ratio tests on the predictors, again
using anova() to compare pairs of models. How have the results changed
from #1?
#no zero-inflation
mod.p <- glmmTMB(total.prey ~ temp.factor * predators, data =col, family = "poisson")
dispersion <- sum(residuals(mod.p, type = "pearson")^2)/(nrow(col) - length(coef(mod.p)))
print(dispersion) #540.574
## [1] 540.574
#From lecture: If you want to test for zero-inflation, fit a zero-inflated Poisson model.
mod.pz = glmmTMB(total.prey ~ temp.factor * predators, data =col, ziformula = ~ 1, family = "poisson")
dispersion2 <- sum(residuals(mod.z, type = "pearson")^2)/(nrow(col) - length(coef(mod.pz)))
## Error: object 'mod.z' not found
print(dispersion2) #1.889824
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'print': object 'dispersion2' not found
#If there is overdispersion, fit a zero-inflated negative binomial model.
mod.bz = glmmTMB(total.prey ~ temp.factor * predators, data =col, ziformula = ~ 1, family = "nbinom2")
summary(mod.bz)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * predators
## Zero inflation: ~1
## Data: col
##
## AIC BIC logLik deviance df.resid
## 673.6 702.9 -322.8 645.6 46
##
##
## Dispersion parameter for nbinom2 family (): 2.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.91754 0.34454 17.175 <2e-16 ***
## temp.factor17_20 0.61443 0.48694 1.262 0.2070
## temp.factor22_25 1.42652 0.76909 1.855 0.0636 .
## predatorsHA+HM 0.41574 0.46205 0.900 0.3682
## predatorsHM 0.10977 0.46219 0.238 0.8123
## predatorsnone 0.03622 0.46223 0.078 0.9375
## temp.factor17_20:predatorsHA+HM -0.68765 0.67112 -1.025 0.3055
## temp.factor22_25:predatorsHA+HM -1.50727 0.91886 -1.640 0.1009
## temp.factor17_20:predatorsHM -0.17328 0.75423 -0.230 0.8183
## temp.factor22_25:predatorsHM -1.45253 1.07754 -1.348 0.1777
## temp.factor17_20:predatorsnone 1.20499 0.69970 1.722 0.0850 .
## temp.factor22_25:predatorsnone 0.21510 0.96061 0.224 0.8228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6191 0.2707 -2.287 0.0222 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Dispersion parameter for nbinom2 family (): 2.12
#Just neg binomial with no zero-inflation
mod.b = glmmTMB(total.prey ~ temp.factor * predators, data =col, family = "nbinom2")
summary(mod.b)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * predators
## Data: col
##
## AIC BIC logLik deviance df.resid
## 734.5 761.7 -354.2 708.5 47
##
##
## Dispersion parameter for nbinom2 family (): 0.188
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.69440 1.03150 5.520 3.38e-08 ***
## temp.factor17_20 0.61436 1.45864 0.421 0.674
## temp.factor22_25 0.04023 1.45876 0.028 0.978
## predatorsHA+HM 0.63890 1.45867 0.438 0.661
## predatorsHM 0.33289 1.45870 0.228 0.819
## predatorsnone 0.25936 1.45872 0.178 0.859
## temp.factor17_20:predatorsHA+HM -0.91071 2.06281 -0.441 0.659
## temp.factor22_25:predatorsHA+HM -0.63180 2.06293 -0.306 0.759
## temp.factor17_20:predatorsHM -1.08947 2.06290 -0.528 0.597
## temp.factor22_25:predatorsHM -1.67560 2.06338 -0.812 0.417
## temp.factor17_20:predatorsnone 0.69425 2.06276 0.337 0.736
## temp.factor22_25:predatorsnone 0.68512 2.06287 0.332 0.740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Dispersion parameter for nbinom2 family (): 0.188
#Compare models
AIC(mod.p)
## [1] 34027.98
## [1] 34027.98
AIC(mod.pz)
## [1] 6478.168
## [1] 6478.168
AIC(mod.bz)
## [1] 673.5606
## [1] 673.5606
AIC(mod.b)
## [1] 734.4738
## [1] 734.4738
#The negative binomial model that accounts for zero inflationwas the best
#Using the best model (i.e., the lowest AIC), perform marginal likelihood ratio tests on the predictors, again using anova() to compare pairs of models. How have the results changed from #1?
#mod.bz = total
mod.bz.temp = glmmTMB(total.prey ~ predators, data =col, ziformula = ~ 1, family = "nbinom2")
mod.bz.pred = glmmTMB(total.prey ~ temp.factor, data =col, ziformula = ~ 1, family = "nbinom2")
anova(mod.bz.temp, mod.bz) #chi 0.03853
## Data: col
## Models:
## mod.bz.temp: total.prey ~ predators, zi=~1, disp=~1
## mod.bz: total.prey ~ temp.factor * predators, zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz.temp 6 673.84 686.41 -330.92 661.84
## mod.bz 14 673.56 702.88 -322.78 645.56 16.281 8 0.03853 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod.bz.pred, mod.bz) #chi 0.0752
## Data: col
## Models:
## mod.bz.pred: total.prey ~ temp.factor, zi=~1, disp=~1
## mod.bz: total.prey ~ temp.factor * predators, zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz.pred 5 671.18 681.65 -330.59 661.18
## mod.bz 14 673.56 702.88 -322.78 645.56 15.623 9 0.0752 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# How have the results changed from #1?
### Based on the best selected model, the negative binomial with the zero-inflation, tempature did have a significant impact on total prey, but predaotrs still did not.
If your results look the way mine do it seems like the three
treatments with predators (HA, HM, HA+HM) may not be very different in
their effect on the prey. Let’s test the effect of predators, but
treating all treatments with predators as the same. Create a new column
that recodes the predator treatment so that it only has two levels:
predators or no predators. Fit a model using the new predator predictor
instead of the old one, and plot fitted effects and perform likelihood
ratio tests as before (you only need to use the ‘best’ probability
distribution model determined in #2). How do these results compare to
what you saw previously? Why do you think the results have changed? How
do you interpret these patterns?
col$pred <- ifelse(col$predators %in% c("HA", "HA+HM", "HM"), "predators", "no predators")
mod.bz2 = glmmTMB(total.prey ~ temp.factor * pred, data =col, ziformula = ~ 1, family = "nbinom2")
summary(mod.bz2)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~1
## Data: col
##
## AIC BIC logLik deviance df.resid
## 665.4 682.2 -324.7 649.4 52
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.486 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000533 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006364 **
## predpredators 0.1672 0.3752 0.446 0.655916
## temp.factor17_20:predpredators -1.5221 0.6039 -2.521 0.011717 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088126 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6191 0.2707 -2.287 0.0222 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.bz.temp2 = glmmTMB(total.prey ~ pred, data =col, ziformula = ~ 1, family = "nbinom2")
mod.bz.pred2 = glmmTMB(total.prey ~ temp.factor, data =col, ziformula = ~ 1, family = "nbinom2")
anova(mod.bz.temp2, mod.bz2) # 0.01125
## Data: col
## Models:
## mod.bz.temp2: total.prey ~ pred, zi=~1, disp=~1
## mod.bz2: total.prey ~ temp.factor * pred, zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz.temp2 4 670.45 678.83 -331.22 662.45
## mod.bz2 8 665.44 682.20 -324.72 649.44 13.005 4 0.01125 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod.bz.pred2, mod.bz2) #0.008333
## Data: col
## Models:
## mod.bz.pred2: total.prey ~ temp.factor, zi=~1, disp=~1
## mod.bz2: total.prey ~ temp.factor * pred, zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz.pred2 5 671.18 681.65 -330.59 661.18
## mod.bz2 8 665.44 682.20 -324.72 649.44 11.739 3 0.008333 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###When taken separately by type of predator, the effect was not significant. But when considering all possible predator treatments together, the effect on total prey was significant.
4. Finally, we have not considered that ‘extra’ zeros could
themselves vary across the experimental treatments. Use an appropriate
zero-inflated model and allow the extra zeros to vary with temperature,
predators, and the interaction between these treatments. Use likelihood
ratio tests to test these terms. The various packages for plotting
fitted effects do not (to my knowledge) have a convenient way to plot
zero-inflated effects (although this could be done ‘manually’ by
extracting model predictions). Instead, look at the model coefficients
returned by summary() to interpret what’s going on with the extra zeros.
It may be easier to interpret the coefficients if you use a model that
removes the non-significant terms. What is the pattern in the extra
zeros according to this model? How does this pattern differ from the
patterns in the count model?
# Use an appropriate zero-inflated model and allow the extra zeros to vary with temperature, predators, and the interaction between these treatments.
mod.bz3 = glmmTMB(total.prey ~ temp.factor * pred, data =col, ziformula = ~ temp.factor * pred, family = "nbinom2")
summary(mod.bz3)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~temp.factor * pred
## Data: col
##
## AIC BIC logLik deviance df.resid
## 656.7 684.0 -315.4 630.7 47
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.485 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000533 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006366 **
## predpredators 0.1671 0.3752 0.445 0.655969
## temp.factor17_20:predpredators -1.5221 0.6039 -2.520 0.011721 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088148 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.15 6435.58 -0.003 0.998
## temp.factor17_20 18.74 6435.58 0.003 0.998
## temp.factor22_25 19.55 6435.58 0.003 0.998
## predpredators 16.51 6435.58 0.003 0.998
## temp.factor17_20:predpredators -16.80 6435.58 -0.003 0.998
## temp.factor22_25:predpredators -16.22 6435.58 -0.003 0.998
mod.bz3.temp = glmmTMB(total.prey ~ temp.factor * pred, data =col, ziformula = ~ pred, family = "nbinom2")
summary(mod.bz3.temp)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~pred
## Data: col
##
## AIC BIC logLik deviance df.resid
## 667.4 686.3 -324.7 649.4 51
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.486 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000532 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006364 **
## predpredators 0.1672 0.3752 0.446 0.655912
## temp.factor17_20:predpredators -1.5221 0.6039 -2.521 0.011717 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088127 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6931 0.5477 -1.266 0.206
## predpredators 0.0984 0.6301 0.156 0.876
mod.bz3.pred = glmmTMB(total.prey ~ temp.factor * pred, data =col, ziformula = ~ temp.factor, family = "nbinom2")
summary(mod.bz3.pred)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~temp.factor
## Data: col
##
## AIC BIC logLik deviance df.resid
## 651.5 672.4 -315.7 631.5 50
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.485 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000533 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006367 **
## predpredators 0.1671 0.3752 0.445 0.655998
## temp.factor17_20:predpredators -1.5221 0.6039 -2.520 0.011723 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088158 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.945 1.026 -2.869 0.00412 **
## temp.factor17_20 2.326 1.128 2.061 0.03929 *
## temp.factor22_25 3.564 1.128 3.158 0.00159 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod.bz3.temp, mod.bz3) # 0.0009122
## Data: col
## Models:
## mod.bz3.temp: total.prey ~ temp.factor * pred, zi=~pred, disp=~1
## mod.bz3: total.prey ~ temp.factor * pred, zi=~temp.factor * pred, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz3.temp 9 667.42 686.27 -324.71 649.42
## mod.bz3 13 656.75 683.98 -315.37 630.75 18.67 4 0.0009122 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod.bz3.pred, mod.bz3) # 0.8644
## Data: col
## Models:
## mod.bz3.pred: total.prey ~ temp.factor * pred, zi=~temp.factor, disp=~1
## mod.bz3: total.prey ~ temp.factor * pred, zi=~temp.factor * pred, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.bz3.pred 10 651.49 672.43 -315.74 631.49
## mod.bz3 13 656.75 683.98 -315.38 630.75 0.7371 3 0.8644
#What is the pattern in the extra zeros according to this model? How does this pattern differ from the patterns in the count model?
summary(mod.bz3)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~temp.factor * pred
## Data: col
##
## AIC BIC logLik deviance df.resid
## 656.7 684.0 -315.4 630.7 47
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.485 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000533 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006366 **
## predpredators 0.1671 0.3752 0.445 0.655969
## temp.factor17_20:predpredators -1.5221 0.6039 -2.520 0.011721 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088148 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.15 6435.58 -0.003 0.998
## temp.factor17_20 18.74 6435.58 0.003 0.998
## temp.factor22_25 19.55 6435.58 0.003 0.998
## predpredators 16.51 6435.58 0.003 0.998
## temp.factor17_20:predpredators -16.80 6435.58 -0.003 0.998
## temp.factor22_25:predpredators -16.22 6435.58 -0.003 0.998
summary(mod.bz3.pred)
## Family: nbinom2 ( log )
## Formula: total.prey ~ temp.factor * pred
## Zero inflation: ~temp.factor
## Data: col
##
## AIC BIC logLik deviance df.resid
## 651.5 672.4 -315.7 631.5 50
##
##
## Dispersion parameter for nbinom2 family (): 1.94
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9538 0.3221 18.485 < 2e-16 ***
## temp.factor17_20 1.8194 0.5253 3.464 0.000533 ***
## temp.factor22_25 1.6416 0.6017 2.728 0.006367 **
## predpredators 0.1671 0.3752 0.445 0.655998
## temp.factor17_20:predpredators -1.5221 0.6039 -2.520 0.011723 *
## temp.factor22_25:predpredators -1.2089 0.7089 -1.705 0.088158 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.945 1.026 -2.869 0.00412 **
## temp.factor17_20 2.326 1.128 2.061 0.03929 *
## temp.factor22_25 3.564 1.128 3.158 0.00159 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###The temperature appears to determine the pattern of zeros in the #model.The coeffcients of the zero-inflation model represent the log-odds the #coefficients can be interpreted as the log of observing an extra zero. #According to the model with ~temp.factor as the zero inflation, the intercept #is-2.945 or
odds <- exp( -2.945)
#The odds are 0.05 or about 1:20.