Mixed effects models

Slide 1

Mixed Effects Modeling

Slide 2: What are Fixed and Random Effects?

• Fixed Effects:
• Parameters associated with the entire population or with certain repeatable levels of experimental units.
• The value for the vessel factor for the SS Titantic.
• Random Effects:
• Parameters associated with individual experimental units drawn at random from a population.
• The value of the vessel factor for the 23rd vessel in a database.

Slide 3: Mixed Effects Models

• A mixed effects model is a model that has both fixed effects and random effects.
• Often a random effect is assigned to observations sharing a common classification factor.
• Goal is to model both variation within GROUPS and variation among GROUPS

Slide 4: A Simple Example

• We select six streams at random and determine the density of fish at each stream three times.
• The questions:
• What is the density in a typical stream?
• What is the variation in density among streams?
• What is the variation in density estimates within a stream?

Slide 5: Three Options for Structured Data

• Your options are:
• 1. Complete pooling: Pool density estimates from all streams for analysis
• 2. No pooling: Include a factor for each stream. This will estimate a separate parameter for each stock.
• 3. Partial pooling: Model the factor as a random effect with some variance around the mean among streams. This separates variation within streams (j = 3 observations) and variation among streams (i = 6 streams)

Slide 6: Complete pooling: Fit of the Standard Linear Model(complete pooling of observations)

This analysis attributed all the error to within-stream variation. There is between-stream variation in density but (1) we have no way to comment on it, and (2) independence violated (are 1000 measurements from 2 streams just as informative about fish density in typical stream as 10 measurements from 200 streams)

• yi,j = β + εi,j
• i = 1,2,..,6; j = 1,2,3
• β̂ = 78.63; σ̂ = 16.03

Slide 7: No pooling: Fit a Model with a “Stream-Effect”

The residual pattern has been removed with fixed effects, but we have no overall mean and no way to comment on between-stream variance, i.e. we cannot say anything about the population of streams.

• yi,j = βi + εi,j
• i = 1,2,..,6; j = 1,2,3
• σ̂ = 6.1

Slide 8: Partial pooling: A Mixed Effects Modelor really random effect model in this simplified case

• We construct a random effects model as follows:
• β is the mean density across the population,
• bi is a random variable representing the deviation from the population mean:
• εi,j is a random variable representing the deviation for observation j from mean density for stream i, i.e.

Within stream variation relative to stream mean

Variation in stream mean

yi,j = β + bi + εi,j

bi ~ N(0; σ2b)

εi,j ~ N(0; σ2)

Slide 9: Finally: A Mixed Effects Model(Some notes-I)

• The residuals are assumed to be independent, normally distributed random variables with constant variance.
• The random effects are assumed to be normally distributed.
• The number of parameters in this model is three (β, σ and σb) because stream is just another level of random variation (with fixed effects model, each stream gets a parameter).

Slide 10: Finally: A Mixed Effects Model

The results from fitting the mixed effects model are:

What can we learn from this?

β = 78.63; σb = 15.79; σ = 6.1

Slide 11: Time for some R

• Fitting the fixed effects models:
• lm1 <- lm(Density~1,data=Streams)
• lm2 <- lm(Density~Stream-1,data=Streams)
• Fitting the mixed effects model:
• lm3 <- lme(Density ~ 1,data=Streams, random = ~ 1 | Stream)
• lme() is from the nlme package
• The call to LME first gives the fixed effects (the intercept in this case), followed by the data set, and finally the random component of the model.

Group

What varies by group

Slide 12: Packages

• The function lmer (in package lme4) also implements linear mixed models (including generalized linear models with glmer).
• lmer doesn’t have all the functionality of lme – however, both should be considered in any application.

Slide 13: The LME Function

• The full call is:
• lme(fixed, data, random, correlation, weights, subset, method, na.action, control)
• fixed – formula for the fixed effects;
• data – must be a data frame;
• random – the random effects (random effects can be nested);
• correlation – correlation structure of the residuals (optional);
• For example “correlation=corAR1()” for autoregressive term
• weights – variance structure of the residuals (optional to account for unequal variance);
• subset – which rows of the data set to include (vector of boolean)
• method – ML or REML
• control – optional control specifications (some of these are pretty useful).

Slide 14: Linear Mixed Models (Notes)

• The linear mixed model can be fit using Maximum Likelihood (ML) or Restricted Maximum Likelihood (REML). If you have many experimental units, this hardly matters, but if you have few replicates, this can make a major difference (use REML).
• REML is default
• The distinction is due to the fact that the MLE for the variance parameter is biased (it has 1/n rather than 1/n-1)
• This is only relevant for Normal distribution (lmm not glmm)
• REML uses an unbiased estimator for the variance
• This is nitpicky, but can matter with small sample size

Slide 15: Linear Mixed Models (Understanding the output)

• Linear mixed-effects model fit by REML
• Data: Streams
• AIC BIC logLik
• 133.7973 136.2969 -63.89865
• Random effects:
• Formula: ~ 1 | Stream
• (Intercept) Residual
• StdDev: 15.78877 6.083881
• Fixed effects: Density ~ 1
• Value Std.Error DF t-value p-value
• (Intercept) 78.62872 6.60332 12 11.90745 0
• Standardized Within-Group Residuals:
• Min Q1 Med Q3 Max
• -1.670431 -0.8239512 0.1256738 0.5295573 1.298799
• Number of Observations: 18
• Number of Groups: 6

• AIC and BIC
• for model selection

The among-stream SD

The within-stream SD

The fixed effects

Slide 16: Fixed and Random Effects Together

• Typically we will have cases in which there are some known sources of variability (the random effects) and some variables that may have different impacts.
• These mixed models will have both fixed and random effects.

Slide 17: A guide to model selection with fixed and random effects (Zuur et al)

• Set the fixed effects based on “as many explanatory variables as possible”.
• Select a random effects structure using AIC (using REML) if variance structure unclear.
• Fix the random effects structure and compare fixed effects options (using ML estimation). [Likelihood tests do not work with REML]
• Compute the final model using REML because ML has biased variance for small sample sizes.

Slide 18: Rule of thumb

• Compare random effects structures with REML.
• Compare fixed effects structures with ML.
• Report final results with REML.

Slide 19: Example:Fixed and Random Effects Together(Length-weight regressions)

• The Problem:
• We select 10 individual fish and measure the relationships between their lengths and weights as they grow. The relationship between length and weight is:
• We have repeated measurements from the same individual “Subject”.
• Should we use a fixed effect for “Subject” or a random effect?
• Modeling decision – do we care about estimate for each individual fish? Not at all. So let’s let intercept be random
• lme(Weight~Length,data=LenW,random = ~1 |Subject,method="REML")

W = aLb; ln Wi,j = ln ai + bi ln Li,j + εi,j

Note that weight and length have been log transformed in the data to make this power function linear

Slide 20: Where does the variation in weight arise from?

Slide 21: Fixed and Random Effects Together(Length-weight regressions)

Slide 22: Fixed and Random Effects Together(Length-weight regressions)

• Linear mixed-effects model fit by REML
• Data: LenW
• AIC BIC logLik
• 6.72662 17.06649 0.63669
• Random effects:
• Formula: ~ 1 | Subject
• (Intercept) Residual
• StdDev: 0.2113517 0.205754
• Fixed effects: Weight ~ Length
• Value Std.Error DF t-value p-value
• (Intercept) -0.026674 0.09560888 89 -0.2790 0.7809
• Length 3.017151 0.02958709 89 101.9753 <.0001
• Standardized Within-Group Residuals:
• Min Q1 Med Q3 Max
• -2.419116 -0.5799786 -0.05204277 0.5719096 2.320844
• Number of Observations: 100
• Number of Groups: 10

Slide 23: Tricks and Traps

• The standard output does not provide:
• the “best estimates” of the group-level random effects - use coef(lmout) to find them.
• the confidence intervals for the variance parameters – use intervals(lmout, 0.95)

Slide 24: Mixed Effects Models

Ben Bolker’s FAQ on GLMMs can help us:

https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html

Slide 25: Another Mixed Model Example

• Dataset: Measure fish trophic position = Dependent variable.
• Individuals from 3 fish species in 6 lakes sampled over a range of body lengths.
• 10 individuals per species.
• Question: Does trophic position increase with length?

Slide 26: How we should not analyze this dataset

• 1. No pooling. Separate regressions for every species in every lake
• estimate a slope and intercept parameter for each regression (2 parameters x 3 species X 6 lakes = 36 parameter estimates) and the sample size for each analysis would be 10.
• We’re also not scientifically interested in the lakes.
• 2. Complete pooling
• Only need to estimate 2 parameters.
• Pseudoreplication? Measurements within same lake or species correlated.
• Noise – Surely some is due to interlake or interspecies variation

Slide 27: Data exploration

Slide 28: Data exploration

Slide 29: Let’s look for structure in residuals

• Complete pooling
• Trophic_Position ~ Length

Lots of structure in residuals grouped by species, but not as bad for lake

Slide 30: Ignore structure in residuals?

• Biased standard errors and confidence intervals (underestimated).
• Inflated type I error for tests of significance.
• Wrong effect estimates if slope varies strongly between groups and you force a single slope.
• Poor predictive performance.

An lm() assumes independent, identically distributed residuals. If fish from the same species or lake are more similar, residuals will be correlated and lm() SEs and p-values will be wrong

Slide 31: Allow intercept to vary by species

lmer(Z_TP~Z_Length + (1|Fish_Species))

Slide 32: Allow intercept to vary by species

• lmer(Z_TP~Z_Length + (1|Fish_Species))
• You might also use
• lm(Z_TP~Z_Length + Fish_Species)
• What’s the difference? Which actually makes more sense here?

Slide 33: Should we use fixed or random effects for groups?

• Ask yourself what you want to infer:
• Do you want estimates/contrasts for the specific species or lakes in your dataset → Fixed effects.
• Do you want to generalize to a larger population of species/lakes (i.e., treat species/lake as a random sample from a population and estimate between-group variance)? → Random effects, if you can estimate them reliably.

Slide 34: Allow intercept to vary by lake

• lmer(Z_TP~Z_Length + (1|Lake))
• Can also do both species and length
• lmer(Z_TP~Z_Length + (1|Fish_Species)+(1|Lake))

Slide 35: But what if the effect of length on TP varies by species or lake?All our prior models were “random intercept” models. Now also “random slope”

lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species)+(1+Z_Length|Lake))

Slide 36: Borrowing strength

• “Shrinkage” toward the mean
• Species with few data rely more on average parameter values
• Here we have balanced design, so not issue
• use of random effects with unbalanced designs + shrinkage can lead to incorrect conclusions if, say, rare species are different from common species.

Slide 37: Which model to use for prediction?

Model 8 (AICc=36.3) and Model 2 (AICc=39.6) much lower than others

REML to compare random effects

Slide 38: Crossed vs Nested random effects

• Note that fish species appear in multiple lakes. This is a crossed random effect, but effects can also be nested
• Consider measurements on student within classes within schools
• If the classes are “remedial”, “average”, and “gifted”, then they could be crossed among schools.
• But if they are just different teachers, then this would be a nested design, R uses different syntax for each.

• Nested
• Use (1|School/Class)

• Crossed
• (1|School) + (1|Class)

Slide 39: Labeling matters!Remove ambiguity with distinct labelling

• Nested Specification
• Use (1|School/Class)

• Crossed Specification
• (1|School) + (1|Class)

• SAME RESULTS
• Both nested due to labelling

Slide 40: Glmer ExampleProportions are Rarely Proportions

• Imagine 4 out of 10 observed animals survive, or contract disease, or give birth
• Does this have the same information content as 40 of 100 animals?
• Now imagine 4 of 10 in one study area, 60 of 100 in another, 2 of 9 in another, and so so. Within study area is a binomial process with probability partially determined by covariates.
• Additional variation not explained by covariates can be modeled with random intercepts

Slide 41: White Stork Example

Oliver Gimenez’s Example

Slide 42

White stork: Is breeding success affected by weather conditions?

• nbchicks = c(151,105,73,107,113,87,77,108,118,122,112,120,122,89,69,71,
• 53,41,53,31,35,14,18)
• nbpairs = c(173,164,103,113,122,112,98,121,132,136,133,137,145,117,90,
• 80,67,54,58,39,42,23,23)
• temp = c(15.1,13.3,15.3,13.3,14.6,15.6,13.1,13.1,15.0,11.7,15.3,14.4,
• 14.4,12.7,11.7,11.9,15.9,13.4,14.0,13.9,12.9,15.1,13.0)
• rainfall = c(67,52,88,61,32,36,72,43,92,32,86,28,57,55,66,26,28,96,48,
• 90,86,78,87)
• year = seq(1,length(nbchicks))

286

Slide 43

Standard GLM on proportions

• nbchicksi ∼ Binomial(nbpairsi, pi)
• logit(pi) = β0 + β1 tempi + β2 rainfalli

287

Slide 44: What if we have substantial additional interannual variation not explained by weather?

Standard GLM on proportions in R

• stork_glm <- glm(cbind(nbchicks,nbpairs-nbchicks) ~ temp + rainfall,
• family = binomial)
• tidy(stork_glm)
• #> # A tibble: 3 x 5
• #> term estimate std.error statistic p.value
• #> <chr> <dbl> <dbl> <dbl> <dbl>
• #> 1 (Intercept) 1.62 0.629 2.58 0.00999
• #> 2 temp 0.0239 0.0465 0.515 0.607
• #> 3 rainfall -0.00656 0.00260 -2.52 0.0116

288

Slide 45

GLMM on proportions with year as a random effect

• nbchicksi ∼ Binomial(nbpairsi, pi)
• logit(pi) = β0i + β1 tempi + β2 rainfalli
• β0i ∼ Normal(μy, σ2y)

289

Slide 46

GLMM on proportions with year as a random effect in R

• m <- glmer(cbind(nbchicks,nbpairs-nbchicks) ~ temp + rainfall
• + (1 | year), family = binomial)

290

Slide 47

GLMM on proportions with year as a random effect in R

• #> Generalized linear mixed model fit by maximum likelihood (Laplace
• #> Approximation) [glmerMod]
• #> Family: binomial ( logit )
• #> Formula: cbind(nbchicks, nbpairs - nbchicks) ~ temp + rainfall + (1 |
• #> year)
• #> AIC BIC logLik deviance df.resid
• #> 159.7418 164.2838 -75.8709 151.7418 19
• #> Random effects:
• #> Groups Name Std.Dev.
• #> year (Intercept) 0.478
• #> Number of obs: 23, groups: year, 23
• #> Fixed Effects:
• #> (Intercept) temp rainfall
• #> 2.478277 -0.026616 -0.008226

291

Slide 48

Slide 49: Example 2 – Nonlinear mixed model for length and age

• We shall assume that:
• Growth follows a von Bertalanffy growth equation

ℓi,j = ℓ∞,i(1 − exp(−κ(ai,j − t0))) + εi,j

Slide 50: Example – length and age

• We shall assume that:
• Growth follows a von Bertalanffy growth equation
• t0 and κ are the same for all animals in the population while the asymptotic size differs among individuals.
• There is normally distributed error.
• R can be used to fit non-linear mixed effects models (function nlme).

ℓi,j = ℓ∞,i(1 − exp(−κ(ai,j − t0))) + εi,j

Slide 51: Example – length and age

The response variable

The random effect

The fixed effects

The covariate

ℓi,j = ℓ∞,i(1 − exp(−κ(ai,j − t0))) + εi,j

Residual variation

Slide 52: Fitting Length and Age Data(Background)

• We have 10 subjects (animals).
• Each animal is measured 10 times between age 1 and age 10.
• Growth follows a von Bertalanffy growth equation where the asymptotic size is subject-specific and the growth rate and age at zero size are the same across subjects.

Slide 53: Fitting Length and Age Data(The basic data)

Let L-inf (saturation) vary by individual

Slide 54: Fitting Length and Age Data(Standard non-linear model-I)

• lm1 <- nls(formula=Length~Linf*(1-exp(-1*Kappa*(Age-Tzero))),
• data=AgeLen, start=c(Linf=100,Kappa=0.2,Tzero=0))

• There is evidence in the
• residuals for model-
• mispecification. This is
• not heteroscedastic errors
• rather it is correlation
• among observations for the
• same animal.

Slide 55: Fitting Length and Age Data(Standard non-linear model-II)

boxplot(split(residuals(lm1),AgeLen$Subject),ylab="Residual",

xlab="Subject",csi=0.2)

Slide 56: The NLME Function-I

• The NLME function has several options:
• nlme(model, data, fixed, random, groups, start, correlation, weights, subset, method,na.action,naPattern, control,verbose)
• model: “response ~ fn(fixed effects)”
• fixed: “fixed = c(fix1~var1,fix2~var2)”.
• There are several ways to specify the fixed effects structure.

• random: “random=c(rnd1 ~ var1, rnd1 ~ 1)”. Leaving out a “random” specification implies all fixed effects have an associated random effect.
• start: “start = c(fix1=x1,fix2=x2,..)” – the initial values for the fixed parameters of the model.

Slide 57: The NLME Function-II

• random: “random=c(rnd1 ~ var1, rnd1 ~ 1)”. Leaving out a “random” specification implies all fixed effects have an associated random effect.
• start: “start = c(fix1=x1,fix2=x2,..)” – the initial values for the fixed parameters of the model.

• nlme(model, data, fixed, random, groups,
• start, correlation, weights, subset, method, control)

Slide 58: Fitting Length and Age Data(Non-linear mixed model-I)

• lm2 <- nlme(model=Length~Linf*(1-exp(-1*Kappa*(Age-Tzero))), data=AgeLen,
• random=Linf~1, fixed=Linf+Kappa+Tzero~1,start=c(Linf=100,Kappa=0.2,Tzero=0))
• Print(lm1$modelStruct)

• Random effect:
• Linf by subject

• Fixed effects:
• Linf population mean
• Kappa
• t0

Slide 59: Fitting Length and Age Data(Non-linear mixed model-II)

plot(lm2,Subject ~ resid(.),abline=0)

• Try “intervals (lm2)” to see
• the confidence intervals

Slide 60: Fitting Length and Age Data(Non-linear mixed model-III)

• Nonlinear mixed-effects model fit by maximum likelihood
• Model: Length ~ Linf * (1 - exp(-1 * Kappa * (Age - Tzero)))
• Data: AgeLen
• AIC BIC logLik
• -53.64647 -40.62062 31.82323
• Random effects:
• Formula: Linf ~ 1 | Subject
• Linf Residual
• StdDev: 8.28837 0.1060163
• Fixed effects: Linf + Kappa + Tzero ~ 1
• Value Std.Error DF t-value p-value
• Linf 97.02240 2.6629267 88 36.4345 0.0000
• Kappa 0.19960 0.0005324 88 374.8752 0.0000
• Tzero -0.00453 0.0040570 88 -1.1176 0.2668
• Correlation:
• Linf Kappa
• Kappa -0.034
• Tzero -0.027 0.876

Slide 61: Fitting Length and Age Data(Non-linear mixed model-IV)

• plot(LinfEst+lm2$coefficients$random$Subject,
• pch=16,csi=0.3)

• plot(lm2, form = resid(., type = "p")
• ~ fitted(.) | Subject,
• abline = 0,pch=16)

Slide 62: Fitting Length and Age Data(Non-linear mixed model-V)

plot(augPred(lm2),csi=0.3,pch=16,lwd=2)

• Plots the data and
• the model
• predictions