For all case studies, follow the general strategy of the workflow chapter (Chapter 12): fit a base model derived from the question, check its residuals, and adjust the model to the problems you find (Table B.10).
C.1 Hurricanes
In Jung et al. (2014), the authors claim that “Female hurricanes are deadlier than male hurricanes”.
Specifically, they analyze the number of hurricane fatalities, and claim that there is an effect of the femininity of the name on the number of fatalities, correcting for several possible confounders. They interpret the result as causal (including mediators), claiming that giving only male names to hurricanes would considerably reduce death toll.
The original model from the paper fits a negative binomial, using mgcv. I suppose the reason is mainly that glmmTMB was not available at the time, and implementations of the negative binomial, in particular MASS::glm.nb and lme4::glmer.nb often had convergence problems.
Confirm that you get the same results as in the paper. It makes sense to translate their model to glmmTMB. Note that the nb parameterization of mgcv corresponds to nbinom2 in glmmTMB. You will get different results when choosing nbinom1
Inspect the fitted model for potential problems, in particular perform a residual analysis of the model, including residuals against all predictors, and improve the model if you find problems.
Forget what they did. Go back to start, do a causal analysis like we did (Chapter 7), and do your own model, diagnosing all residual problems that we discussed (Section 2.4, Section 9.3, Chapter 10, Chapter 11). Do you think there is an effect of femininity?
TipSolution
This is the model fit by Jung et al. (2014), fit with glmmTMB
Note that the model above scales only NDAM (and the paper does not say whether they scaled), but as we are looking for main effects in the presence of interactions, we should at least center the predictors (Section 2.10.3); here, we also scale them, which makes the effects comparable.
now main effect is n.s.; it’s a bit dodgy, but if you read in the main paper, they do not claim a significant main effect, they mainly argue via ANOVA and significance at high values of NDAM, so let’s run an ANOVA. car::Anova() gives type II tests, which, unlike type III, do not depend on the centering (see Section 3.2.4 and Caution 3.1):
In the ANOVA we see that MasFem still n.s. but interactions, and if you would calculate effect of MasFem at high NDAM, it is significant. Something like that is argued in the paper. We can emulate this by changing NDAM centering to high NDAM, which gives us a p-value for the main effect of MasFem at high values of NDAM
Now we see the significant main effect that they report. Note, however, that the significant difference is only there for high NDAM, i.e. what we do here is to project the effect of the interaction on the main effect. An alternative to do the same thing would be an effects plot (Section 2.3), or to specifically use predict() (Section 2.11) to calculate differences and CIs at high NDAM values.
OK, this means we can replicate the results of the paper, even if concentrating the entire analysis exclusive on high NDAM seems a bit cherry-picking. Choosing the evaluation point after looking at the data inflates the false-positive risk, and at the maximum NDAM the estimate rests on very few storms. Another way to phrase the result is that we don’t find a main effect of MasFem. However, to be fair: the current results do say that there is a significant difference at high NDAM, and such a difference, if it existed, would be important.
But we haven’t done residual checks yet. Let’s do that for the original model of the paper (originalModelGAM):
res <-simulateResiduals(originalModelGAM)plot(res)
No significant deviation in the general DHARMa plot, but residuals ~ NDAM looks funny, which was also pointed out by Bob O’Hara in a blog post after publication of the paper. Let’s try to correct this - transforming NDAM with ^0.2 (found by trying a few powers; 0.1 to 0.2 fit best) does a great job:
All gone, only damage is doing the effect. As the interactions are n.s., the main effects would hardly change with a different centering.
What if we would have fit our own model? First of all, note that if hurricane names were given randomly, we wouldn’t have to worry about confounders (Section 7.2.1). However, this is not the case, hurricanes were only named randomly after 1978 or so.
plot(MasFem ~ Year, data = hurricanes)
So, we could either take the earlier data out, which would remove half of our data, or we have to worry about confounding with variables that change over time. The most obvious thing would be to take time itself (Year) in the model, to correct for temporal confounding.
Do we need other variables that are not confounders? There are two reasons to add them (the first is the idea of precision variables, Tip 8.1):
they have strong effects on the response - not adding them could lead to residual problems and increase residual variance, which increases uncertainties and cost power
we want to fit interactions.
I added NDAM to the model, because we saw earlier that it has a strong effect. I think it’s not unreasonable to check for an interaction as well.
As we have several observations per year, a conservative approach would be to add a RE on year. Note that we use year both as a fixed effect (to remove temporal trends) and a random intercept, which is perfectly fine, however.
newModel =glmmTMB(alldeaths ~ sMasFem * sNDAM02 + Year + (1|Year),data = hurricanes, family = nbinom2)summary(newModel)
Family: nbinom2 ( log )
Formula: alldeaths ~ sMasFem * sNDAM02 + Year + (1 | Year)
Data: hurricanes
AIC BIC logLik -2*log(L) df.resid
630.8 648.4 -308.4 616.8 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Year (Intercept) 2.571e-07 0.0005071
Number of obs: 92, groups: Year, 49
Dispersion parameter for nbinom2 family (): 1.11
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.542287 12.730846 -0.200 0.842
sMasFem 0.073207 0.119273 0.614 0.539
sNDAM02 1.309624 0.118106 11.089 <2e-16 ***
Year 0.002426 0.006423 0.378 0.706
sMasFem:sNDAM02 0.179874 0.117191 1.535 0.125
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The variance of the year random effect is estimated as practically zero (a boundary fit), so there is no extra year-to-year variation beyond the fixed effect of Year, and removing the random effect would not change the results. The results remain that there is no effect of MasFem!
The hurricanes are also a time series. The check for temporal autocorrelation and an AR1 model for these data are shown in the hurricane exercise of the chapter on correlation structures (Section 11.2.1); they show only a weak, negative autocorrelation that could also be a statistical fluke.
C.2 Researcher degrees of freedom — skin color and red cards
Silberzahn et al. (2018) published a “meta-analysis” in Advances in Methods and Practices in Psychological Science, where they had provided 29 teams with the same data set to answer one research question: “[W]hether soccer players with dark skin tone are more likely than those with light skin tone to receive red cards from referees”.
Spoiler: They found that the “[a]nalytic approaches varied widely across the teams, and the estimated effect sizes ranged from 0.89 to 2.93 (Mdn = 1.31) in odds-ratio units”, highlighting that different approaches in data analysis can yield significant variation in the results.
The data is in
library(EcoData)?redCards
CautionTasks
Do a re-analysis of the data as if you were the 30th team to contribute the results to the meta-analysis. You can find the data in the ecodata package, dataset redCards.
Response variable: ‘redCards’ (+‘yellowReds’?).
primary predictors: ‘rater1’, ‘rater2’
Multiple variables, potentially accounting for confounding, offsetting, grouping, … are included in the data.
The rater variable contains ratings of “two independent raters blind to the research question who, based on their profile photo, categorized players on a 5-point scale ranging from (1) very light skin to (5) very dark skin. Note that ‘rater1’ and ‘rater2’ in EcoData::redCards are already rescaled to the range 0 … 1, as described in the paper (”This variable was rescaled to be bounded by 0 (very light skin) and 1 (very dark skin) prior to the final analysis, to ensure consistency of effect sizes across the teams of analysts. The raw ratings were rescaled to 0, .25, .50, .75, and 1 to create this new scale.”)
When you’re done, have a look at the other modeling teams. Do you understand the models they fit? Note that the results are displayed in terms of odds ratios. Are your results within the range of estimates from the 29 teams in Silberzahn et al. (2018)?
TipSolution
Base model
There is no single “correct” analysis - that is the whole point of Silberzahn et al. (2018). The following is therefore not a fully worked-out solution, but one defensible base model to start from. Each row is a player-referee dyad, so red cards are a count out of the number of games the two shared; we use that as exposure via a log offset, take the mean of the two skin-tone ratings rescaled to 0-1 as the predictor, and add random effects for player and referee to respect the non-independence.
library(EcoData); library(glmmTMB); library(DHARMa)d = redCardsd$skin = (d$rater1 + d$rater2)/2# mean of the two ratings (already on the 0..1 scale)d = d[!is.na(d$skin) & d$games >0, ]m =glmmTMB(redCards ~ skin + (1|playerShort) + (1|refNum),offset =log(games), family = poisson, data = d)summary(m)
exp(fixef(m)$cond["skin"]) # rate ratio, darkest (1) vs lightest (0) skin
skin
1.390424
res =simulateResiduals(m); plot(res); testDispersion(res)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.0768, p-value = 0.008
alternative hypothesis: two.sided
The exponentiated skin coefficient is our effect in rate-ratio units (comparable to the paper’s odds ratios: red cards are rare, about 0.4% per player-referee game, and for rare events a rate ratio is close to an odds ratio). Compare it to the 0.89-2.93 range of odds ratios reported across the 29 teams. This bare model gives a rate ratio of about 1.4, which is inside the range of the 29 teams and close to their median (1.31). The study shows, however, that the “researcher degrees of freedom” (choice of exposure, random-effect structure, which rater, count vs. binomial, and - crucially - which confounders are included) can move the result substantially. From this base model, you would go on to address the slight overdispersion (Section 9.3) and to add plausible confounders (see Chapter 7; league, playing position, player skill, …) and see how much the estimate changes.
Together with the orientation of the Maze, this makes 3 possible directional confounders, and the main predictor (if the Ant went to the pheromone in the first visit).
The treatments and orientations are nearly balanced by the setup, but the main predictor is a behavior of the ants, and ants do not behave randomly. Its distribution across the design variables is not balanced:
Most ants (75%) went to the pheromone on their first visit, but this share differs between treatments (83% in LL, where the pheromone was on the left side in both runs, against 71 to 74% in the others). The main predictor is therefore weakly but visibly associated with the design variables (see Section 3.2.3.2 for why this matters even in balanced designs).
Adding an RE on colony is logical, and then let’s run the model:
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: went.phero ~ ant_group_main + directionConst + directionPhero +
Orientation + (1 | Colony)
Data: dat
AIC BIC logLik -2*log(L) df.resid
2192.5 2225.7 -1090.2 2180.5 1869
Scaled residuals:
Min 1Q Median 3Q Max
-4.6870 -1.0451 0.4836 0.6876 1.0318
Random effects:
Groups Name Variance Std.Dev.
Colony (Intercept) 0.9284 0.9636
Number of obs: 1875, groups: Colony, 15
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.9505 0.3330 5.857 4.72e-09 ***
ant_group_mainSecondvisit_1st_to_phero 0.1311 0.1248 1.050 0.29352
directionConstTRUE -1.1465 0.2409 -4.760 1.94e-06 ***
directionPheroright -0.5680 0.1984 -2.863 0.00420 **
Orientationright -0.3493 0.1222 -2.859 0.00425 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) a__S_1 dCTRUE drctnP
ant_g_S_1__ -0.252
drctnCnTRUE -0.413 -0.062
drctnPhrrgh -0.441 -0.001 0.306
Orinttnrght -0.182 -0.093 -0.074 0.085
Surprisingly, we find large effects of the other variables. Because of these large effects, we also test for interactions with the experimental treatment
We now find that there is an interaction with the main predictor, and there could be effects. We can also look at this visually.
plot(allEffects(fit2))
The results are difficult to interpret. The large effects of the design variables suggest that the ants did not behave randomly with respect to the maze (a bias of the setup). The main predictor is a behavior of the same ants and weakly associated with these design variables (see the table above), so its effect is hard to separate from them.
For our education, we can also look at the residual plots. I will use fit1, because there was a misfit:
res <-simulateResiduals(fit1, plot = T)
As we would expect significant interactions, we would probably see something if we plot residuals against predictors or their interactions, but I want to show you something else:
We will not see dispersion problems in a 0/1 binomial, but actually, this is a k/n binomial, just that the data are not prepared as such.
Either way, in DHARMa, you can aggregate residuals by a grouping variable.
res2 <-recalculateResiduals(res, group = dat$Colony)
Now, we essentially check k/n data, and we see that there is overdispersion, which is caused by the misfit.
plot(res2)
testDispersion(res2)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.0233, p-value = 0.832
alternative hypothesis: two.sided
Let’s do the same for model 2, which included the interactions.
res <-simulateResiduals(fit2, plot = T)
res2 <-recalculateResiduals(res, group = dat$Colony)plot(res2)
Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
: Fitting terminated with step failure - check results carefully
testDispersion(res2)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.0422, p-value = 0.736
alternative hypothesis: two.sided
Look at the Snails data set in the EcoData package, and find out which environmental and/or seasonal predictors i) explain the total abundance and ii) the infection rate of Bulinus spp.
The snails data set in the EcoData package includes observations on the distribution of freshwater snails and their infection rates ( schistosomiasis (a parasite)).
The first scientific question is that their abundance depends on the water conditions. The second scientific question is that their infection rate depends on the water conditions and seasonal factors
The data also contains data on other environmental (and seasonal factors). You should consider if it is useful to add them to the analysis.
Number of infected individuals: BP_pos_tot, BF_pos_tot, BT_pos_tot
Total abundance of Bulinus spp.: Bulinus_tot
Total number of infected Bulinus spp.: Bulinus_pos_tot
CautionTasks
Model the total abundance of Bulinus spp. (Bulinus_tot)
Model the infection rate of Bulinus spp. (Bulinus_pos_tot) (k/n binomial)
Optional: Model the species individually (BP_tot, BF_tot, BT_tot)
Optional: Fit a multivariate (joint) species distribution model (see Section 11.6)
TipSolution
Total abundance of Bulinus (GLMM)
Prepare+scale data:
library(lme4)library(glmmTMB)library(DHARMa)data = EcoData::snailsdata$sTemp_Water =scale(data$Temp_Water)data$spH =scale(data$pH)data$swater_speed_ms =scale(data$water_speed_ms)data$swater_depth =scale(data$water_depth)data$sCond =scale(data$Cond)data$swmo_prec =scale(data$wmo_prec)data$syear =scale(data$year)data$sLat =scale(data$Latitude)data$sLon =scale(data$Longitude)data$sTemp_Air =scale(data$Temp_Air)# remove rows with NAs in the variables used (see text)rows =rownames(model.matrix(Bulinus_tot~sTemp_Water + spH + sLat + sLon + sCond + seas_wmo+ swmo_prec + swater_speed_ms + swater_depth +sTemp_Air+ syear + duration + locality + site_irn + coll_date, data = data))data = data[rows, ]
We remove the rows with NAs in the variables we use beforehand (with model.matrix()), so that all models use the same rows and can be compared (Section 6.1). We do the same for the seed bank data.
Our hypothesis is that the abundance of Bulinus species depends on the water characteristics, e.g. site_type, Temp_water, pH, Cond, swmo_prec, water_speed_ms, and water_depth. We will set the length of the collection duration as an offset.
As the sites are nested within localities, we will set a nested random intercept (see Tip 4.2) on site_irn within locality. Also potential confounders are collection date (coll_date), the season (wet or dry months, seas_wmo), year, and maybe other environmental factors such as the air temperature?.
DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
Does not look great -> dispersion problems -> switch to -> negative binomial distribution:
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.50462, p-value = 0.016
alternative hypothesis: two.sided
testZeroInflation(res)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.0533, p-value = 0.072
alternative hypothesis: two.sided
So far, we have not included a spatial trend, so there could be spatial autocorrelation. Let’s check for it. Note that we simulate unconditionally here (simulateREs = "unconditional"): with the default conditional simulations, the fitted random effects of the sites would absorb the spatial pattern, and the test could not detect it (Tip 4.4).
## Spatialres =simulateResiduals(model4, simulateREs ="unconditional")res2 =recalculateResiduals(res, group =c(data$site_irn))groupLocations =aggregate(cbind(data$sLat, data$sLon ), list( data$site_irn), mean)testSpatialAutocorrelation(res2, x = groupLocations$V1, y = groupLocations$V2)
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.295785, expected = -0.011364, sd = 0.066951, p-value =
4.482e-06
alternative hypothesis: Distance-based autocorrelation
Significant! Following the rule from Section 11.1.1, we first try to remove a large-scale spatial trend, here with a quadratic surface in latitude and longitude:
res =simulateResiduals(model5, simulateREs ="unconditional")res2 =recalculateResiduals(res, group =c(data$site_irn))testSpatialAutocorrelation(res2, x = groupLocations$V1, y = groupLocations$V2)
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.067369, expected = -0.011364, sd = 0.066933, p-value =
0.2395
alternative hypothesis: Distance-based autocorrelation
AIC(model4, model5)
df AIC
model4 25 14138.24
model5 30 14131.57
After de-trending, there is no significant spatial autocorrelation left, and the AIC also prefers the de-trended model, so we don’t need a spatial correlation structure. If autocorrelation had remained, we would add one, e.g. exp(0 + numFac | group) (see Section 11.3). Final residual check:
res =simulateResiduals(model5, plot =TRUE)
Residuals do not look perfect but I would say that we can stop here now.
TipSolution
Prevalence of Bulinus (GLMM)
We use the data prepared and scaled in the solution for the total abundance (run that one first; it also loads lme4, glmmTMB and DHARMa). Let’s start directly with all potential confounders:
DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
The dispersion test is n.s. (if k/n data are overdispersed, use family = betabinomial or an observation-level random effect, Section 9.3.2).
Check for spatial autocorrelation (again with unconditional simulations, as for the abundance model):
## Spatialres =simulateResiduals(model1, simulateREs ="unconditional")res2 =recalculateResiduals(res, group =c(data$site_irn))groupLocations =aggregate(cbind(data$sLat, data$sLon ), list( data$site_irn), mean)testSpatialAutocorrelation(res2, x = groupLocations$V1, y = groupLocations$V2)
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.159982, expected = -0.011364, sd = 0.067031, p-value =
0.01058
alternative hypothesis: Distance-based autocorrelation
Significant! As before, we first try to remove a large-scale spatial trend:
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.098452, expected = -0.011364, sd = 0.067040, p-value =
0.1014
alternative hypothesis: Distance-based autocorrelation
AIC(model1, model2)
df AIC
model1 21 1323.300
model2 26 1328.605
After de-trending, the spatial autocorrelation is no longer significant, so we don’t need a spatial correlation structure. Note, however, that the AIC does not support the trend terms, i.e. the spatial signal was weak; keeping model1 would also be defensible.
TipSolution
Multivariate (joint) species distribution model
The species models are connected by their response to latent variable (unobserved environment). For that, we will transform our dataset with respect to species from wide (sp1, sp2, sp3) to long format (species abundances in one column and a second column telling us the group (species)). In the model then, we will separate the species and their responses by using ~0+Species + Species:(predictors).
The latent variable structure is set by the rr(…) object in the formula. We again use the data prepared in the solution for the total abundance (run that one first):
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: logSBDensity ~ sSeedMass + sSeedShape + sSeedN + (1 | Site) +
(1 | Species)
Data: data
REML criterion at convergence: 7531.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.6261 -0.4945 -0.0823 0.4542 3.1525
Random effects:
Groups Name Variance Std.Dev.
Species (Intercept) 3.0724 1.7528
Site (Intercept) 0.4705 0.6859
Residual 3.6790 1.9181
Number of obs: 1729, groups: Species, 152; Site, 17
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.5314 0.2271 42.4393 6.743 3.23e-08 ***
sSeedMass -0.3316 0.1273 175.6509 -2.604 0.00999 **
sSeedShape -0.3175 0.1546 155.2402 -2.054 0.04168 *
sSeedN -0.1818 0.1289 165.3141 -1.410 0.16028
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) sSdMss sSdShp
sSeedMass -0.037
sSeedShape -0.034 0.168
sSeedN -0.031 0.039 0.077
Altitude and humidity may influence both which species (and thus which seed traits) occur at a site and the seed bank density, so we add them as potential confounders:
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: logSBDensity ~ sSeedMass + sSeedShape + sSeedN + sAltitude +
sHum + (1 | Site) + (1 | Species)
Data: data
REML criterion at convergence: 7510.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.6468 -0.5089 -0.0720 0.4388 3.2140
Random effects:
Groups Name Variance Std.Dev.
Species (Intercept) 3.03405 1.7419
Site (Intercept) 0.06929 0.2632
Residual 3.68099 1.9186
Number of obs: 1729, groups: Species, 152; Site, 17
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.65385 0.16664 113.38609 9.925 < 2e-16 ***
sSeedMass -0.32886 0.12662 176.33088 -2.597 0.0102 *
sSeedShape -0.32165 0.15372 155.68190 -2.092 0.0380 *
sSeedN -0.18051 0.12816 165.92065 -1.408 0.1609
sAltitude -0.56452 0.10125 17.10538 -5.575 3.27e-05 ***
sHum 0.14280 0.09738 14.04198 1.466 0.1646
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) sSdMss sSdShp sSeedN sAlttd
sSeedMass -0.052
sSeedShape -0.044 0.168
sSeedN -0.042 0.039 0.077
sAltitude -0.058 0.007 -0.012 -0.006
sHum -0.004 -0.008 -0.003 0.004 0.512
Adding altitude and humidity barely changes the seed-trait estimates: seed mass and seed shape have negative effects on seed bank density, seed number has no clear effect.
Question:
What about the germination temperatures, T50/T95?
Answer:
They could be mediators (Seed Shape -> T95 -> SBDensity). For the total effect of the seed traits we leave them out; we would include them only to estimate the direct effect (Section 7.3).
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00691219 (tol = 0.002, component 1)
See ?lme4::convergence and ?lme4::troubleshooting.
summary(model1)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SBPA ~ sSeedMass + sSeedShape + sSeedN + sAltitude + sHum + (1 |
Site) + (sAltitude | Species)
Data: data
AIC BIC logLik -2*log(L) df.resid
1368.9 1423.4 -674.4 1348.9 1719
Scaled residuals:
Min 1Q Median 3Q Max
-2.5936 -0.2658 -0.1017 0.2001 3.3347
Random effects:
Groups Name Variance Std.Dev. Corr
Species (Intercept) 8.16725 2.8578
sAltitude 4.38865 2.0949 -0.13
Site (Intercept) 0.06798 0.2607
Number of obs: 1729, groups: Species, 152; Site, 17
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.6285 0.4066 -6.465 1.01e-10 ***
sSeedMass -1.3670 0.6431 -2.126 0.0335 *
sSeedShape -0.5813 0.3100 -1.875 0.0608 .
sSeedN -2.1359 1.3165 -1.622 0.1047
sAltitude -1.3935 0.3364 -4.142 3.44e-05 ***
sHum -0.1025 0.1308 -0.784 0.4330
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) sSdMss sSdShp sSeedN sAlttd
sSeedMass 0.221
sSeedShape 0.033 0.062
sSeedN 0.336 0.060 -0.066
sAltitude 0.132 0.142 0.079 0.012
sHum 0.023 0.002 -0.003 -0.011 0.229
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00691219 (tol = 0.002, component 1)
See ?lme4::convergence and ?lme4::troubleshooting.
Model did not converge, but there is a trick which often helps (see Section 4.4.1). The default optimizer in lme4 is not the best optimizer, changing it to ‘bobyqa’ often helps with convergence issues
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: SBPA ~ sSeedMass + sSeedShape + sSeedN + sAltitude + sHum + (1 |
Site) + (sAltitude | Species)
Data: data
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik -2*log(L) df.resid
1368.9 1423.4 -674.4 1348.9 1719
Scaled residuals:
Min 1Q Median 3Q Max
-2.5938 -0.2658 -0.1018 0.2001 3.3353
Random effects:
Groups Name Variance Std.Dev. Corr
Species (Intercept) 8.16321 2.8571
sAltitude 4.38183 2.0933 -0.13
Site (Intercept) 0.06786 0.2605
Number of obs: 1729, groups: Species, 152; Site, 17
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.6284 0.4066 -6.465 1.01e-10 ***
sSeedMass -1.3693 0.6437 -2.127 0.0334 *
sSeedShape -0.5816 0.3100 -1.876 0.0606 .
sSeedN -2.1419 1.3182 -1.625 0.1042
sAltitude -1.3920 0.3362 -4.141 3.46e-05 ***
sHum -0.1022 0.1307 -0.781 0.4345
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) sSdMss sSdShp sSeedN sAlttd
sSeedMass 0.222
sSeedShape 0.033 0.062
sSeedN 0.337 0.061 -0.066
sAltitude 0.132 0.142 0.079 0.012
sHum 0.023 0.002 -0.003 -0.011 0.229
Success, it converged! Note that silencing the warning by switching the optimizer is no evidence that the estimates are good: check that estimates agree across optimizers (allFit(model1)) and that the fit is not singular (isSingular(model1)).
Residual checks:
Check residuals:
res =simulateResiduals(model1, plot=TRUE)
Residuals look good!
Bonus: With phylogenetic correlation structure:
(This bonus part uses phyl.upd2 and the function as.theta.vcov() from the SBDensity bonus solution above, so run that one first.)
Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
Scaled residual values: 0.5748496 0.5526574 0.05684348 0.6166831 0.1098403 0.2653392 0.2043353 0.531952 0.3925722 0.2693407 0.07212984 0.2086342 0.4148425 0.5437874 0.5376287 0.2171818 0.7868062 0.3029527 0.5870516 0.7049796 ...
Conditional simulations:
pred =predict(model6, re.form =NULL, type ="response")simulations =sapply(1:1000, function(i) rbinom(length(pred),1, pred))res =createDHARMa(simulations, model.frame(model6)[,1], pred)plot(res)
C.6 Snouter
CautionTasks
Fit one of the (continuous) responses in the snouter dataset, e.g. snouter1.1, against the predictors rain + djungle (see ?snouter); the coordinates are X and Y. Check for spatial autocorrelation and proceed to fitting a spatial model if needed (Section 11.3). See the data set’s help for details on the variables.
library(EcoData)str(snouter)
TipSolution
Fit a base model, then check its residuals for spatial autocorrelation. We use snouter1.1 (a continuous response) here; the coordinates in this dataset are X and Y:
testSpatialAutocorrelation(res, x = snouter$X, y = snouter$Y)
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = 0.12556787, expected = -0.00090334, sd = 0.00141608, p-value
< 2.2e-16
alternative hypothesis: Distance-based autocorrelation
If the spatial test is significant, we first de-trend with a 2D spline (de-trend before adding a correlation structure, Section 11.1.1):
library(mgcv)fit2 =gam(snouter1.1~ rain + djungle +te(X, Y), data = snouter)res2 =simulateResiduals(fit2)testSpatialAutocorrelation(res2, x = snouter$X, y = snouter$Y)
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.03162372, expected = -0.00090334, sd = 0.00141607, p-value
< 2.2e-16
alternative hypothesis: Distance-based autocorrelation
This changes residuals somewhat, but there is still significant spatial autocorrelation left. Thus, we have to move to a spatial model.
# TODO
C.7 Volcano Island
CautionTasks
Analyze ?volcanoisland in the EcoData package. Perform an appropriate causal statistical analysis (see Chapter 7; including residual checks and everything) to understand the predictors that determine the value of
windObs (numeric)
lizards (1/0)
beetles2 (counts)
survivedOf20 (k out of 20)
There are various predictors, including plot, year, x, y that can be considered. Assume that all the issues that we talked about (random effects, Chapter 4; overdispersion, Section 9.3.2; zero-inflation, Section 9.3.3; spatial autocorrelation, Section 11.3) could appear in this data.
Hint: it is possible that the possibility to observe a species depends on wind, and it is possible that lizards eat beetles, so you can also consider these variables (or predictions for these variables) as predictors of the others.
Hint 2: binomial k/n data in R is specified as a response that is cbind(k, n-k), so you would write glm(cbind(survivedOf20, 20 - survivedOf20) ~ predictor, data = data, family = “binomial”)
Hint 3: lizards may be present but not observed (imperfect detection, see Hint 1). Occupancy models, which separate the true presence from the detection process, are introduced with the same data set in the chapter on occupancy models of Hartig (2024).
This is an open, integrative exercise, so the following is not a fully worked-out solution but a base model to start from for one of the responses (beetle counts). The other responses (windObs -> Gaussian / gls, lizards -> binomial, survivedOf20 -> k/n binomial) follow exactly the same recipe.
library(EcoData); library(glmmTMB); library(DHARMa)d = volcanoisland# beetles: counts -> start with Poisson and plot as RE, then diagnosem0 =glmmTMB(beetles2 ~ windObs + lizards + (1|plot), family = poisson, data = d)res =simulateResiduals(m0); plot(res)
testDispersion(res); testZeroInflation(res)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 2.6752, p-value < 2.2e-16
alternative hypothesis: two.sided
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.1159, p-value < 2.2e-16
alternative hypothesis: two.sided
# if over-dispersed / zero-inflated -> negative binomial + ZI termm1 =glmmTMB(beetles2 ~ windObs + lizards + (1|plot),family = nbinom2, ziformula =~1, data = d)# then check for residual spatial autocorrelation on the coordinatesres1 =simulateResiduals(m1)grp =recalculateResiduals(res1, group = d$plot)testSpatialAutocorrelation(grp,x =aggregate(x ~ plot, d, mean)$x,y =aggregate(y ~ plot, d, mean)$y)
DHARMa Moran's I test for distance-based autocorrelation
data: grp
observed = 0.205367, expected = -0.010101, sd = 0.017216, p-value <
2.2e-16
alternative hypothesis: Distance-based autocorrelation
Work through the diagnostic ladder (Table B.10) in order - functional form -> distribution -> dispersion -> zero-inflation -> correlation - adding structure only when a residual check demands it. Note the hint that lizards (predators) may belong in the beetle model, and that detectability may depend on windObs.
