Mixed model lab #1

Licensed under the Creative Commons attribution-noncommercial license. Please share & remix noncommercially, mentioning its origin.

Preliminaries

Load packages:

## modeling packages
library("lme4")     ## basic (G)LMMs
library("nlme")     ## more LMMs
library("glmmTMB")
library("afex")     ## helper functions
library("emmeans")

## graphics
library("ggplot2"); theme_set(theme_bw())
## squash panels together
zero_margin <- theme(panel.spacing=grid::unit(0,"lines")) 
library("lattice")
library("dotwhisker")  ## coefficient plots

## data manipulation
library("broom.mixed")

Linear mixed model: starling example

Data from Toby Marthews, r-sig-mixed-models mailing list: you can find it here.

• subject: individual bird
• roostsitu (tree, nest-box, inside, other): location of nest box
• mnth (Nov, Jan): month
• stmass: mass of bird

Graphics: spaghetti plots

## cosmetic tweaks
load("data/starling.RData")   ## loads "dataf"
ggplot(dataf,aes(x=mnth,y=stmass))+
  geom_point()+
  geom_line(aes(group=subject))+  ## connect subjects with a line
  facet_grid(.~roostsitu)+        ## 1 row of panels by roost
  zero_margin                         ## squash together

You could also try

ggplot(dataf,aes(mnth,stmass,colour=roostsitu))+
    geom_point()+
    geom_line(aes(group=subject))

Overkill for this data set, but sometimes it can be useful to put every individual in its own facet:

## reorder individuals by magnitude of difference
dataf <- transform(dataf, subject=reorder(subject,stmass,FUN=diff))
ggplot(dataf,aes(mnth,stmass,colour=roostsitu,group=subject))+
    geom_point()+geom_line()+facet_wrap(~subject)+
    zero_margin

It’s pretty obvious that the starting (November) mass varies among roost situations (tree/nest-box/etc.), and that mass increases from November to January, but we might like to know if the slopes differ among situations. That means our fixed effects would be ~roostsitu*mnth, with our attention focused on the roostsitu:mnth (interaction) term. For random effects, we can allow both intercept (obviously) and slope (maybe) to vary among individuals, via (1+mnth|subject) or equivalent … in this case, because measurements are only taken in two months, we can also write the random term as (1|subject/mnth).

However, it turns out that we can’t actually estimate random slopes for this model, because every individual is only measured twice. That means that the variance in the slopes would be completely confounded with the residual variance.

If we fit this with lme4:

lmer1 <- lmer(stmass~mnth*roostsitu+(1|subject/mnth),data=dataf)

## Error: number of levels of each grouping factor must be < number of observations (problems: mnth:subject)

We get an error. We can actually do this (lmerControl(check.nobs.vs.nlev="ignore"); see ?lmerControl for details), and it will actually work, but it will pick an arbitrary division between the residual variance and the month:subject variance.

In principle we can do this model with glmmTMB, by dropping the residual variance term (for technical reasons, this isn’t possible in lme4):

glmmTMB1 <- glmmTMB(stmass~mnth*roostsitu+(1|subject/mnth),
                    data = dataf,
                    ## REML = TRUE,  ## figure out why this gives a npd Hessian ...
                    dispformula = ~0)

But now let’s forget about including the (unidentifiable) random slopes.

lmer2 <- lmer(stmass~mnth*roostsitu+(1|subject),data=dataf)

We can now get estimates, although the subject-level random effects are very uncertain: see confint(lmer2)

Walking through the output:

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: stmass ~ mnth * roostsitu + (1 | subject)
##    Data: dataf

• This gives us the formula, reminds us that the model was fitted by restricted maximum likelihood, and gives us the “REML criterion” (equivalent to -2 log likelihood for a model fitted by maximum likelihood).

## REML criterion at convergence: 429.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max

This is a quick reminder of the scaled residuals; these should be approximately unskewed (median \(\approx 0\)) and the extremes should be somewhere around \(\pm 2\) (bigger in a big data set )

## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept)  0.3444  0.5869

• This tells us about the random effect - it tells us the variances and standard deviations of the random effect (the standard deviations are not an estimate of the uncertainty of the estimate – the estimate itself is a variance, or standard deviation), and the variance and standard deviation of the residual variance. We can see that the standard deviation of the intercept among subjects is much less than the standard deviation of the residual error …

##  Residual             17.3500  4.1653

• The listed number of observations and groups is very useful for double-checking that the random effects grouping specification is OK. Are the numbers of groups what we expected from the experimental design?

## 
## Fixed effects:
##                           Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)                 83.600      1.330 71.973  62.847  < 2e-16 ***
## mnthJan                      7.200      1.863 36.000   3.865 0.000446 ***
## roostsitunest-box           -4.200      1.881 71.973  -2.233 0.028685 *  
## roostsituinside             -5.000      1.881 71.973  -2.658 0.009682 ** 
## roostsituother              -8.200      1.881 71.973  -4.359 4.27e-05 ***
## mnthJan:roostsitunest-box    3.600      2.634 36.000   1.367 0.180244

• the standard fixed-effect output. lmer does not tell us the denominator degrees of freedom for the test (although we can get a rough idea of importance/significance fro the \(t\) statistics; e.g. \(t>2.5\) will be significant at \(p<0.05\) for 6 or more degrees of freedom). We’ll come back to this in the inference section.

Diagnostic plots

performance::check_model() gives a relatively complete set of model diagnostics.

performance::check_model(lmer2)

• linearity is automatically satisfied when all predictors are categorical and all interactions are included
• similarly, in a balanced factorial design all observations have the same leverage
• I usually ignore the collinearity diagnostics
• residuals are slightly light-tailed, not a cause for concern
• random effect estimates (BLUPs/conditional modes) look Normal

DHARMa (diagnostics via simulated residuals: usually unnecessary for LMMs):

sr2 <- DHARMa::simulateResiduals(lmer2)
plot(sr2)

lme4 also has some built-in plotting capabilities (skipped during live-coding session)

Basic diagnostic plots: fitted vs. residuals, coloured according to the roostsitu variable (col=dataf$roostsitu), point type according to month (pch=).

plot(lmer2,col=dataf$roostsitu,pch=dataf$mnth,id=0.05)

The id argument specifies to mark points that are beyond the specified Normal confidence limits (i.e. \(\alpha\)-values) (these are not outliers in any formal statistical sense but might be points you want to look at if they are otherwise interesting). By default they are labeled by random-effect group; use idLabels to change the default labels, and idLabels=~.obs to label by observation number.

You can get a scale-location plot by specifying

## FIXME: why doesn't this match performance scale-resid plot?
ss <- residuals(lmer2)/(1-hatvalues(lmer2))
xyplot(sqrt(abs(ss)) ~ fitted(lmer2),type=c("p","smooth"))

Adding the smoothed line is helpful because uneven sampling can influence your perception of the pattern.

Q-Q plot (a straight line indicates normality)

qqmath(lmer2,col=dataf$roostsitu,pch=dataf$mnth)

We can also do this with broom.mixed::augment() and stat_qq() (although there are also some some limitations - colour argument doesn’t actually work…)

ggplot(augment(lmer2),
       aes(sample=.resid/sd(.resid)))+  ## scale to variance=1
    stat_qq(aes(group=1,colour=roostsitu))+
    geom_abline(intercept=0,slope=1)

## Warning: The following aesthetics were dropped during statistical
## transformation: colour.
## ℹ This can happen when ggplot fails to infer the correct grouping
##   structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a
##   numerical variable into a factor?

There are some deviations here, but not enough that I would worry very much. In particular, the distribution is slightly thin-tailed (the smallest residuals are largest than expected, and the largest residuals are smaller than expected), which would make the results slightly conservative (a fat-tailed distribution would make them anticonservative).

Boxplots of residuals subdivided by roostsitu (you have to put the grouping variable on the left side of the formula here):

plot(lmer2,roostsitu~resid(.))

Might be easier with augment:

aa <- augment(lmer2)
ggplot(aa,aes(roostsitu,.resid))+
    geom_boxplot()+coord_flip()

Plot random effects to look for outliers:

dotplot(ranef(lmer2))

## $subject

Or with tidy + ggplot:

tt <- tidy(lmer2,effects="ran_vals")
tt <- transform(tt,level=reorder(level,estimate))
ggplot(tt,aes(level,estimate))+
    geom_pointrange(aes(ymin=estimate-1.96*std.error,
                        ymax=estimate+1.96*std.error))+
    coord_flip()

In general, try plotting residuals and \(\sqrt{\textrm{abs}(r)}\) both as a function of the fitted values and as a function of particularly important predictors.

• if you have your own example data set handy that you can write a simple model for, run these diagnostics and interpret them
• the hatvalues() function can also be useful (although all hat values are identical here)
• also see the influence.ME package. Try plot(influence(lmer2,group="subject")) and interpret the results …

Try it yourself: tundra carbon example

• Download the file tundra_agg.rda (from the data directory given above). The GS.NEE records the net ecosystem exchange (in grams of carbon/m^2/year) in a given year and site; cYear is a centered version of the year; n records the number of observations taken in a given site/year (the NEE value given is an average of these observations).
• Draw some plots (with ggplot2 if you can, otherwise with base R graphics) to inspect the data
• We are interested in the change in NEE over time (i.e., a random-slopes model). Run an lmer fit (note this data set does not have the issue of the starling data set), draw some preliminary conclusions, and inspect the plot diagnostics.