GLMM FAQ

an informal GLMM FAQ list for the r-sig-mixed-models mailing list

Introduction

This is an informal FAQ list for the r-sig-mixed-models mailing list.

The most commonly used functions for mixed modeling in R are

• linear mixed models: aov(), nlme::lme¹, lme4::lmer; brms::brm
• generalized linear mixed models (GLMMs)
  • frequentist: MASS::glmmPQL, lme4::glmer; glmmTMB
  • Bayesian: MCMCglmm::MCMCglmm; brms::brm
• nonlinear mixed models: nlme::nlme, lme4::nlmer; brms::brm
• GNLMMs: brms::brm

Another quick-and-dirty way to search for mixed-model related packages on CRAN:

grep("l.?m[me][^t]",rownames(available.packages(repos = 'https://cran.us.r-project.org')),value=TRUE)

 [1] "blmeco"         "buildmer"       "cellVolumeDist"
 [4] "climextRemes"   "elementR"       "glmertree"     
 [7] "glmmboot"       "glmmEP"         "glmmfields"    
[10] "glmmLasso"      "glmmML"         "glmmSeq"       
[13] "glmmsr"         "glmmTMB"        "lamme"         
[16] "lme4"           "lmec"           "lmeInfo"       
[19] "lmem.qtler"     "lmeNB"          "lmeNBBayes"    
[22] "lmeresampler"   "lmerTest"       "lmeSplines"    
[25] "lmeVarComp"     "lmmot"          "lmmpar"        
[28] "lrmest"         "lsmeans"        "mailmerge"     
[31] "mlmm.gwas"      "mlmmm"          "mvglmmRank"    
[34] "nlmeODE"        "nlmeU"          "palmerpenguins"
[37] "tlmec"          "vagalumeR"     

There are some false positives here (e.g. palmerpenguins); see here if you’re interested in “regex golf”.

Other sources of help

• the mailing list is r-sig-mixed-models@r-project.org
  • sign up here
  • archives here
  • or Google search with the tag site:https://stat.ethz.ch/pipermail/r-sig-mixed-models/
• The source code of this document is available on GitHub; the rendered (HTML) version lives on GitHub pages.
• Searching on StackOverflow with the [r] [mixed-models] tags, or on CrossValidated with the [mixed-model] tag may be helpful (these sites also have an [lme4] tag).

DISCLAIMERS:

• (G)LMMs are hard - harder than you may think based on what you may have learned in your second statistics class, which probably focused on picking the appropriate sums of squares terms and degrees of freedom for the numerator and denominator of an \(F\) test. ‘Modern’ mixed model approaches, although more powerful (they can handle more complex designs, lack of balance, crossed random factors, some kinds of non-Normally distributed responses, etc.), also require a new set of conceptual tools. In order to use these tools you should have at least a general acquaintance with classical mixed-model experimental designs but you should also, probably, read something about modern mixed model approaches. @littell_sas_2006 and @pinheiro_mixed-effects_2000 are two places to start, although Pinheiro and Bates is probably more useful if you want to use R. Other useful references include @gelman_data_2006 (focused on Bayesian methods) and @zuur_mixed_2009. If you are going to use generalized linear mixed models, you should understand generalized linear models (@dobson_introduction_2008, @faraway_extending_2006, and @McCullaghNelder1989 are standard references; the last is the canonical reference, but also the most challenging).
• All of the issues that arise with regular linear or generalized-linear modeling (e.g.: inadequacy of p-values alone for thorough statistical analysis; need to understand how models are parameterized; need to understand the principle of marginality and how interactions can be treated; dangers of overfitting, which are not mitigated by stepwise procedures; the non-existence of free lunches) also apply, and can apply more severely, to mixed models.
• When SAS (or Stata, or Genstat/AS-REML or …) and R differ in their answers, R may not be wrong. Both SAS and R may be `right’ but proceeding in a different way/answering different questions/using a different philosophical approach (or both may be wrong …)
• The advice in this FAQ comes with absolutely no warranty of any sort.

References

linear mixed models

web/open

• UCLA IDRE statistical consulting
• @barr_learning_2020 Chapters 5-8

books (dead-tree/closed)

• pinheiro_mixed-effects_2000: LMM only.
• @zuur_mixed_2009: Focused on ecology.
• @gelman_data_2006: LMM and GLMM; Bayesian; examples from social science. Intermediate mathematics.
• (Rethinking)

Model definition

Model specification

The following formula extensions for specifying random-effects structures in R are used by

• lme4
• nlme (nested effects only, although crossed effects can be specified with more work)
• glmmADMB and glmmTMB

MCMCglmm uses a different specification, inherited from AS-REML.

(Modified from Robin Jeffries, UCLA:)

formula	meaning
(1|group)	random group intercept
(x|group) = (1+x|group)	random slope of x within group with correlated intercept
(0+x|group) = (-1+x|group)	random slope of x within group: no variation in intercept
(1|group) + (0+x|group)	uncorrelated random intercept and random slope within group
(1|site/block) = (1|site)+(1|site:block)	intercept varying among sites and among blocks within sites (nested random effects)
site+(1|site:block)	fixed effect of sites plus random variation in intercept among blocks within sites
(x|site/block) = (x|site)+(x|site:block) = (1 + x|site)+(1+x|site:block)	slope and intercept varying among sites and among blocks within sites
(x1|site)+(x2|block)	two different effects, varying at different levels
x*site+(x|site:block)	fixed effect variation of slope and intercept varying among sites and random variation of slope and intercept among blocks within sites
(1|group1)+(1|group2)	intercept varying among crossed random effects (e.g. site, year)

Or in a little more detail:

equation	formula
\(ß_0 + ß_{1}X_{i} + e_{si}\)	n/a (Not a mixed-effects model)
\((ß_0 + b_{S,0s}) + ß_{1}X_i + e_{si}\)	~ X + (1|Subject)
\((ß_0 + b_{S,0s}) + (ß_{1} + b_{S,1s}) X_i + e_{si}\)	~ X + (1 + X|Subject)
\((ß_0 + b_{S,0s} + b_{I,0i}) + (ß_{1} + b_{S,1s}) X_i + e_{si}\)	~ X + (1 + X|Subject) + (1|Item)
As above, but \(S_{0s}\), \(S_{1s}\) independent	~ X + (1|Subject) + (0 + X| Subject) + (1|Item)
\((ß_0 + b_{S,0s} + b_{I,0i}) + ß_{1}X_i + e_{si}\)	~ X + (1|Subject) + (1|Item)
\((ß_0 + b_{I,0i}) + (ß_{1} + b_{S,1s})X_i + e_{si}\)	~ X + (0 + X|Subject) + (1|Item)

Modified from: http://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet?lq=1 (Livius)

The magic development version of the equatiomatic package can handle mixed models (remotes::install_github("datalorax/equatiomatic")), e.g.

library(lme4)
library(equatiomatic)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
equatiomatic::extract_eq(fm1)

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{1j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \]

It doesn’t handle GLMMs (yet), but you could fit two fake models — one LMM like your GLMM but with a Gaussian response, and one GLM with the same family/link function as your GLMM but without the random effects — and put the pieces together.

More possibly useful links:

• Rense Nieuwenhuis’s blogpost/lesson on lme4 model specification
• CrossValidated’s lmer cheat sheet
• Kristoffer Magnusson’s Using R and lme/lmer to fit different two- and three-level longitudinal models

Should I treat factor xxx as fixed or random?

This is in general a far more difficult question than it seems on the surface. There are many competing philosophies and definitions. For example, from @gelman_analysis_2005:

Before discussing the technical issues, we briefly review what is meant by fixed and random effects. It turns out that different—in fact, incompatible—definitions are used in different contexts. [See also Kreft and de Leeuw (1998), Section 1.3.3, for a discussion of the multiplicity of definitions of fixed and random effects and coefficients, and Robinson (1998) for a historical overview.] Here we outline five definitions that we have seen: 1. Fixed effects are constant across individuals, and random effects vary. For example, in a growth study, a model with random intercepts αi and fixed slope β corresponds to parallel lines for different individuals i, or the model yit = αi + βt. Kreft and de Leeuw [(1998), page 12] thus distinguish between fixed and random coefficients. 2. Effects are fixed if they are interesting in themselves or random if there is interest in the underlying population. Searle, Casella and McCulloch [(1992), Section 1.4] explore this distinction in depth. 3. “When a sample exhausts the population, the corresponding variable is fixed; when the sample is a small (i.e., negligible) part of the population the corresponding variable is random” [Green and Tukey (1960)]. 4. “If an effect is assumed to be a realized value of a random variable, it is called a random effect” [LaMotte (1983)]. 5. Fixed effects are estimated using least squares (or, more generally, maximum likelihood) and random effects are estimated with shrinkage [“linear unbiased prediction” in the terminology of Robinson (1991)]. This definition is standard in the multilevel modeling literature [see, e.g., Snijders and Bosker (1999), Section 4.2] and in econometrics.

Another useful comment (via Kevin Wright) reinforcing the idea that “random vs. fixed” is not a simple, cut-and-dried decision: from @schabenberger_contemporary_2001, p. 627:

Before proceeding further with random field linear models we need to remind the reader of the adage that one modeler’s random effect is another modeler’s fixed effect.

@clark2015should address this question from a mostly econometric perspective, focusing mostly on practical variance/bias/RMSE criteria.

One point of particular relevance to ‘modern’ mixed model estimation (rather than ‘classical’ method-of-moments estimation) is that, for practical purposes, there must be a reasonable number of random-effects levels (e.g. blocks) – more than 5 or 6 at a minimum. This is not surprising if you consider that random effects estimation is trying to estimate an among-block variance. For example, from @Crawley2002 p. 670:

Are there enough levels of the factor in the data on which to base an estimate of the variance of the population of effects? No, means [you should probably treat the variable as] fixed effects.

Some researchers (who treat fixed vs random as a philosophical rather than a pragmatic decision) object to this approach.

Also see a very thoughtful chapter in @hodges_richly_2016.

Treating factors with small numbers of levels as random will in the best case lead to very small and/or imprecise estimates of random effects; in the worst case it will lead to various numerical difficulties such as lack of convergence, zero variance estimates, etc.. (A small simulation exercise shows that at least the estimates of the standard deviation are downwardly biased in this case; it’s not clear whether/how this bias would affect the point estimates of fixed effects or their estimated confidence intervals.) In the classical method-of-moments approach these problems may not arise (because the sums of squares are always well defined as long as there are at least two units), but the underlying problems of lack of power are there nevertheless.

Also see this thread on the r-sig-mixed-models mailing list.

Nested or crossed?

• Relatively few mixed effect modeling packages can handle crossed random effects, i.e. those where one level of a random effect can appear in conjunction with more than one level of another effect. (This definition is confusing, and I would happily accept a better one.) A classic example is crossed temporal and spatial effects. If there is random variation among temporal blocks (e.g. years) ‘’and’’ random variation among spatial blocks (e.g. sites), ‘’and’’ if there is a consistent year effect across sites and ‘’vice versa’’, then the random effects should be treated as crossed.
• lme4 does handled crossed effects, efficiently; if you need to deal with crossed REs in conjunction with some of the features that nlme offers (e.g. heteroscedasticity of residuals via weights/varStruct, correlation of residuals via correlation/corStruct, see p. 163ff of @pinheiro_mixed-effects_2000 (section 4.2.2: Google books link)
• I rarely find it useful to think of fixed effects as “nested” (although others disagree); if for example treatments A and B are only measured in block 1, and treatments C and D are only measured in block 2, one still assumes (because they are fixed effects) that each treatment would have the same effect if applied in the other block. (One might like to estimate treatment-by-block interactions, but in this case the experimental design doesn’t allow it; one would have to have multiple treatments measured within each block, although not necessarily all treatments in every block.) One would code this analysis as response~treatment+(1|block) in lme4. Also, in the case of fixed effects, crossed and nested specifications change the parameterization of the model, but not anything else (e.g. the number of parameters estimated, log-likelihood, model predictions are all identical). That is, in R’s model.matrix function (which implements a version of Wilkinson-Rogers notation) a*b and a/b (which expand to 1+a+b+a:b and 1+a+a:b respectively) give model matrices with the same number of columns.
• Whether you explicitly specify a random effect as nested or not depends (in part) on the way the levels of the random effects are coded. If the ‘lower-level’ random effect is coded with unique levels, then the two syntaxes (1|a/b) (or (1|a)+(1|a:b)) and (1|a)+(1|b) are equivalent. If the lower-level random effect has the same labels within each larger group (e.g. blocks 1, 2, 3, 4 within sites A, B, and C) then the explicit nesting (1|a/b) is required. It seems to be considered best practice to code the nested level uniquely (e.g. A1, A2, …, B1, B2, …) so that confusion between nested and crossed effects is less likely.

(When) can I include a predictor as both fixed and random?

See blog post by Thierry Onkelinx

Model extensions

Overdispersion

Testing for overdispersion/computing overdispersion factor

• with the usual caveats, plus a few extras – counting degrees of freedom, etc. – the usual procedure of calculating the sum of squared Pearson residuals and comparing it to the residual degrees of freedom should give at least a crude idea of overdispersion. The following attempt counts each variance or covariance parameter as one model degree of freedom and presents the sum of squared Pearson residuals, the ratio of (SSQ residuals/rdf), the residual df, and the \(p\)-value based on the (approximately!!) appropriate \(\chi^2\) distribution. Do PLEASE note the usual, and extra, caveats noted here: this is an APPROXIMATE estimate of an overdispersion parameter. Even in the GLM case, the expected deviance per point equaling 1 is only true as the distribution of individual deviates approaches normality, i.e. the usual \(\lambda>5\) rules of thumb for Poisson values and \(\textrm{min}(Np, N(1-p)) > 5\) for binomial values (e.g. see @venables_modern_2002, p. 208-209). (And that’s without the extra complexities due to GLMM, i.e. the “effective” residual df should be large enough to make the sums of squares converge on a \(\chi^2\) distribution …)
• Remember that (1) overdispersion is irrelevant for models that estimate a scale parameter (i.e. almost anything but Poisson or binomial: Gaussian, Gamma, negative binomial …) and (2) overdispersion is not estimable (and hence practically irrelevant) for Bernoulli models (= binary data = binomial with \(N=1\)).
• The recipes below may need adjustment for some of the more complex model types allowed by glmmTMB (e.g. zero-inflation/variable dispersion), where it’s less clear what to measure to estimate overdispersion.

The following function should work for a variety of model types (at least glmmADMB, glmmTMB, lme4, …).

overdisp_fun <- function(model) {
    rdf <- df.residual(model)
    rp <- residuals(model,type="pearson")
    Pearson.chisq <- sum(rp^2)
    prat <- Pearson.chisq/rdf
    pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
    c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

Example:

library(lme4)
library(glmmTMB)

set.seed(101)  
d <- data.frame(x=runif(1000),
                f=factor(sample(1:10,size=1000,replace=TRUE)))
suppressMessages(d$y <- simulate(~x+(1|f), family=poisson,
                          newdata=d,
                          newparams=list(theta=1,beta=c(0,2)))[[1]])
m1 <- glmer(y~x+(1|f),data=d,family=poisson)
overdisp_fun(m1)

       chisq        ratio          rdf            p 
1035.9966325    1.0391140  997.0000000    0.1902294 

m2 <- glmmTMB(y~x+(1|f),data=d,family="poisson")
overdisp_fun(m2)

       chisq        ratio          rdf            p 
1035.9961394    1.0391135  997.0000000    0.1902323 

The gof function in the aods3 provides similar functionality (it reports both deviance- and \(\chi^2\)-based estimates of overdispersion and tests).

Fitting models with overdispersion?

• quasilikelihood estimation: MASS::glmmPQL. Quasi- was deemed unreliable in lme4, and is no longer available. (Part of the problem was questionable numerical results in some cases; the other problem was that DB felt that he did not have a sufficiently good understanding of the theoretical framework that would explain what the algorithm was actually estimating in this case.) geepack::geelgm may be workable (haven’t tried it)

  If you really want quasi-likelihood analysis for glmer fits, you can do it yourself by adjusting the coefficient table - i.e., by multiplying the standard error by the square root of the dispersion factor ² and recomputing the \(Z\)- and \(p\)-values accordingly, as follows:

## extract summary table; you may also be able to do this via
##  broom::tidy or broom.mixed::tidy
quasi_table <- function(model,ctab=coef(summary(model)),
                           phi=overdisp_fun(model)["ratio"]) {
    qctab <- within(as.data.frame(ctab),
    {   `Std. Error` <- `Std. Error`*sqrt(phi)
        `z value` <- Estimate/`Std. Error`
        `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE)
    })
    return(qctab)
}
printCoefmat(quasi_table(m1),digits=3)

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.2277     0.2700    0.84      0.4    
x             2.0640     0.0528   39.11   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## to use this with glmmTMB, we need to separate out the
##  conditional component of the summary
printCoefmat(quasi_table(m2,
                         ctab=coef(summary(m2))[["cond"]]),
             digits=3)

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.2277     0.2700    0.84      0.4    
x             2.0640     0.0528   39.09   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Another version, this one tidyverse-centric:

library(broom.mixed)
library(dplyr)
tidy_quasi <- function(model, phi=overdisp_fun(model)["ratio"],
                       conf.level=0.95) {
    tt <- (tidy(model, effects="fixed")
        %>% mutate(std.error=std.error*sqrt(phi),
                   statistic=estimate/std.error,
                   p.value=2*pnorm(abs(statistic), lower.tail=FALSE))
    )
    return(tt)
}
tidy_quasi(m1)

# A tibble: 2 x 6
  effect term        estimate std.error statistic p.value
  <chr>  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 fixed  (Intercept)    0.228    0.270      0.843   0.399
2 fixed  x              2.06     0.0528    39.1     0    

tidy_quasi(m2)

# A tibble: 2 x 7
  effect component term        estimate std.error statistic p.value
  <chr>  <chr>     <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 fixed  cond      (Intercept)    0.228    0.270      0.843   0.399
2 fixed  cond      x              2.06     0.0528    39.1     0    

These functions make some simplifying assumptions: (1) this overdispersion computation is approximate

In this case using quasi-likelihood doesn’t make much difference, since the data we simulated in the first place were Poisson.) Keep in mind that once you switch to quasi-likelihood you will either have to eschew inferential methods such as the likelihood ratio test, profile confidence intervals, AIC, etc., or make more heroic assumptions to compute “quasi-” analogs of all of the above (such as QAIC).

• observation-level random effects (OLRE: this approach should work in most packages). If you want to a citation for this approach, try @elston_analysis_2001, who cite @lawson_disease_1999; apparently there is also an example in section 10.5 of @maindonald_data_2010, and (according to an R-sig-mixed-models post) this is also discussed by @rabehesketh_multilevel_2008. Also see @browne_variance_2005 for an example in the binomial context (i.e. logit-normal-binomial rather than lognormal-Poisson). Agresti’s excellent (2002) book @agresti_categorical_2002 also discusses this (section 13.5), referring back to @breslow_extrapoisson_1984 and @hinde_compound_1982. [Notes: (a) I haven’t checked all these references myself, (b) I can’t find the reference any more, but I have seen it stated that observation-level random effect estimation is probably dodgy for PQL approaches as used in Elston et al 2001]
• alternative distributions
  • Poisson-lognormal model for counts or binomial-logit-Normal model for proportions (see above, “observation-level random effects”)
  • negative binomial for counts or beta-binomial for proportions
    • lme4::glmer.nb() should fit a negative binomial, although it is somewhat slow and fragile compared to some of the other methods suggested here. lme4 cannot fit beta-binomial models (these cannot be formulated as a part of the exponential family of distributions)
    • glmmTMB will fit two parameterizations of the negative binomial: family="nbinom2" gives the classic parameterization with \(\sigma^2=\mu(1+\mu/k)\) (“NB2” in Hardin and Hilbe’s terminology) while family="nbinom1" gives a parameterization with \(\sigma^2=\phi \mu\), \(\phi>1\) (“NB1” to Hardin and Hilbe). The latter might also be called a “quasi-Poisson” parameterization because it matches the mean-variance relationship assumed by quasi-Poisson models, i.e. the variance is strictly proportional to the mean (although the proportionality constant must be >1, a limitation that does not apply to quasi-likelihood approaches). (glmmADMB will also fit these models, with family="nbinom" for NB2, but is deprecated in favour of glmmTMB.)
    • glmmTMB allows beta-binomial models ([@harrison_comparison_2015] suggests comparing beta-binomial with OLRE models to assess reliability)
    • the brms package has a negbinomial family (no beta-binomial, but it does have a wide range of other families)
• other packages/approaches (less widely used, or requiring a bit more effort)
  • gamlss.mx:gamlssNP
  • WinBUGS/JAGS (via R2WinBUGS/Rjags)
  • AD Model Builder (possibly via R2admb package) or TMB
  • gnlmm in the repeated package (off-CRAN)
  • ASREML

Negative binomial models in glmmTMB and lognormal-Poisson models in glmer (or MCMCglmm) are probably the best quick alternatives for overdispersed count data. If you need to explore alternatives (different variance-mean relationships, different distributions), then ADMB, TMB, WinBUGS, Stan, NIMBLE are the most flexible alternatives.

Underdispersion

Underdispersion (much less variability than expected) is a less common problem than overdispersion.

• mild overdispersion is sometimes ignored, since it tends in general to lead to conservative rather than anti-conservative results
• quasi-likelihood (and the quasi-hack listed above) can handle under- as well as underdispersion
• some other solutions exist, but are less widely implemented
  • for distributions with a small range (e.g. litter sizes of large mammals), one can treat responses as ordinal (e.g. using the ordinal package, or MCMCglmm or brms for Bayesian solutions)
  • the COM-Poisson distribution and generalized Poisson distributions, implemented in glmmTMB, can handle underdispersion (J. Hilbe recommends the latter in this CrossValidated answer). (VGAM has a generalized Poisson distribution, but doesn’t handle random effects.)

Gamma GLMMs

While one (well, OK I) would naively think that GLMMs with Gamma distributions would be just as easy (or hard) as any other sort of GLMMs, it seems that they are in fact harder to implement. Basic simulated examples of Gamma GLMMs can fail in lme4 despite analogous problems with Poisson, binomial, etc. distributions. Solutions: - the default inverse link seems particularly problematic; try other links (especially family=Gamma(link="log")) if that is possible/makes sense - consider whether a lognormal model (i.e. a regular LMM on logged data) would work/makes sense. - @lo_transform_2015 argue that the Gamma family with an identity link is superior to lognormal models for reaction-time data. I (BMB) don’t find their argument particularly convincing, but lots of people want to do this. Unfortunately this is technically challenging (see here), because it is likely that some “illegal” values (predicted responses \(\le 0\)) will occur while fitting the model, even if the final fitted model makes no impossible predictions. Thus something has to be done to make the model-fitting machinery tolerant of such values (i.e. returning NA for these model evaluations, or clamping illegal values to the constrained space with an appropriate smooth penalty function).

Gamma models can be fitted by a wide variety of platforms (lme4::glmer, MASS::glmmPQL, glmmADMB, glmmTMB, MixedModels.jl, MCMCglmm, brms … not sure about others.

Beta GLMMs

Proportion data where the denominator (e.g. maximum possible number of successes for a given observation) is not known can be modeled using a Beta distribution. @smithson_better_2006 is a good introduction for non-statisticians (not in the mixed-model case), and the betareg package [@cribari-neto_beta_2009] handles non-mixed Beta regressions. The glmmTMB and brms packages handle Beta mixed models (brms also handles zero-inflated and zero-one inflated models).

Zero-inflation

See e.g. @martin_zero_2005 or @warton_many_2005 (“many zeros does not mean zero inflation”) or @zuur_zero-truncated_2009 for general information on zero-inflation.

Count data

• MCMCglmm handles zero-truncated, zero-inflated, and zero-altered models, although specifying the models is a little bit tricky: see Sections 5.3 to 5.5 of the CourseNotes vignette
• glmmADMB handles
  • zero-inflated models (with a single zero-inflation parameter – i.e., the level of zero-inflation is assumed constant across the whole data set)
  • truncated Poisson and negative binomial distributions (which allows two-stage fitting of hurdle models)
• glmmTMB handles a variety of Z-I and Z-T models (allows covariates, and random effects, in the zero-alteration model)
• brms does too
• so does GLMMadaptive
• Gavin Simpson has a detailed writeup showing that mgcv::gam() can do simple mixed models (Poisson, not NB) with zero-inflation, and comparing mgcv with glmmTMB results
• gamlssNP in the gamlss.mx package should handle zero-inflation, and the gamlss.tr package should handle truncated (i.e. hurdle) models – but I haven’t tried them
• roll-your-own: ADMB/R2admb, WinBUGS/R2WinBUGS, TMB, Stan, …

Continuous data

Continuous data are a special case where the mixture model for zero-inflated data is less relevant, because observations that are exactly zero occur with probability (but not probability density) zero. There are two cases of interest:

Probability density of \(x\) zero or infinite

In this case zero is a problematic observation for the distribution; it’s either impossible or infinitely (locally) likely. Some examples:

• Gamma distribution: probability density at zero is infinite (if shape<1) or zero (if shape>1); it’s finite only for an exponential distribution (shape==1)
• Lognormal distribution: the probability density at zero is zero.
• Beta distribution: the probability densities at 0 and 1 are zero (if the corresponding shape parameter is >1) or infinite (if shape<1)

The best solution depends very much on the data-generating mechanism.

• If the bad (0/1) values are generated by rounding (e.g. proportions that are too close to the boundaries are reported as being on the boundaries), the simplest solution is to “squeeze” these in slightly, e.g. \(y \to (y +a)/2a\) for some sensible value of \(a\) [@smithson_better_2006]
• If you think that zero values are generated by a separate process, the simplest solution is to fit a Bernoulli model to the zero/non-zero data, then a conditional continuous model for the non-zero values; this is effectively a hurdle model.
• you might have censored data where all values below a certain limit (e.g. a detection limit) are recorded as zero; in this case you might be able to use survreg() and frailty() in the survival package for random-intercept models (as suggested on r-help by Thomas Lumley in 2003 or on StackOverflow by user 42- in 2014. The lmec package handles linear mixed models.
• The cplm package handles ‘Tweedie compound Poisson linear models’, which in a particular range of parameters allows for skewed continuous responses with a spike at zero

Probability density of \(x\) positive and finite

In this case (e.g. a spike of zeros in the center of an otherwise continuous distribution), the hurdle model probably makes the most sense.

Tests for zero-inflation

• you can use a likelihood ratio test between the regular and zero-inflated version of the model, but be aware of boundary issues (search “boundary” elsewhere on this page …) – the null value (no zero inflation) is on the boundary of the feasible space
• you can use AIC or variations, with the same caveats
• you can use Vuong’s test, which is often recommended for testing zero-inflation in GLMs, because under some circumstances the various model flavors under consideration (hurdle vs zero-inflated vs “vanilla”) are not nested. Vuong’s test is implemented (and referenced) in the pscl package, and should be feasible to implement for GLMMs, but I don’t know of an implementation. Someone should let me (BMB) know if they find one.
• two untested but reasonable approaches:
  • use a simulate() method if it exists to construct a simulated distribution of the proportion of zeros expected overall from your model, and compare it to the observed proportion of zeros in the data set
  • compute the probability of a zero for each observation. On the basis of (conditionally) independent Bernoulli trials, compute the expected number of zeros and the confidence intervals – compare it with the observed number.

Spatial and temporal correlation models, heteroscedasticity (“R-side” models)

In nlme these so-called R-side (R for “residual”) structures are accessible via the weights/VarStruct (heteroscedasticity) and correlation/corStruct (spatial or temporal correlation) arguments and data structures. This extension is a bit harder than it might seem. In LMMs it is a natural extension to allow the residual error terms to be components of a single multivariate normal draw; if that MVN distribution is uncorrelated and homoscedastic (i.e. proportional to an identity matrix) we get the classic model, but we can in principle allow it to be correlated and/or heteroscedastic.

It is not too hard to define marginal correlation structures that don’t make sense. One class of reasonably sensible models is to always assume an observation-level random effect (as MCMCglmm does for computational reasons) and to allow that random effect to be MVN on the link scale (so that the full model is lognormal-Poisson, logit-normal binomial, etc., depending on the link function and family).

For example, a relatively simple Poisson model with spatially correlated errors might look like this:

\[ \begin{split} \eta & \sim \textrm{MVN}(a + b x, \Sigma) \\ \Sigma_{ij} & = \sigma^2 \exp(-d_{ij}/s) \\ y_i & \sim \textrm{Poisson}(\lambda=\exp(\eta_i)) \end{split} \]

That is, the marginal distributions of the response values are Poisson-lognormal, but on the link (log) scale the latent Normal variables underlying the response are multivariate normal, with a variance-covariance matrix described by an exponential spatial correlation function with scale parameter \(s\).

How can one achieve this?

• These types of models are not implemented in lme4, for either LMMs or GLMMs; they are fairly low priority, and it is hard to see how they could be implemented for GLMMs (the equivalent for LMMs is tedious but should be straightforward to implement).
• For LMMs, you can use the spatial/temporal correlation structures that are built into (n)lme
• You can use the spatial/temporal correlation structures available for (n)lme, which include basic geostatistical (space) and ARMA-type (time) models.

  library(sos)
  findFn("corStruct")