Lecture 5 linear mixed effects models
Download
1 / 28

Lecture 5 Linear Mixed Effects Models - PowerPoint PPT Presentation


  • 224 Views
  • Uploaded on

Advanced Research Skills. Lecture 5 Linear Mixed Effects Models. Olivier MISSA, [email protected] Outline. Explore options available when assumptions of classical linear models are untenable.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about ' Lecture 5 Linear Mixed Effects Models' - hanley


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Lecture 5 linear mixed effects models

Advanced Research Skills

Lecture 5 Linear Mixed Effects Models

Olivier MISSA, [email protected]


Outline
Outline

  • Explore options available when assumptions of classical linear models are untenable.

  • In this lecture: What can we do when observations (and thus residuals) are not strictly independent ?


Classical linear models
Classical Linear Models

  • Defined by three assumptions:

  • (1) the response variable is continuous.

  • (2) the residuals (ε) are normally distributed and ...

  • (3) ... independently and identically distributed.

  • Today, we will consider a range of options availablewhen we either know or suspect that our data are not strictly independent from each other.

  • (Departures from the other assumptions will be dealt with later)


Non independent residuals

A non-linear trend

Non-independent Residuals

  • In previous lectures:

  • We merely checked the independence of our residuals by inspecting the plot of residualsvs.fitted values.

  • Example from lecture 2:

which suggested that our linear model was probably misspecified


Non independent observations
Non-independent Observations

  • Data collection can often lead to non-independence among your observations.

  • A few examples:

  • Repeated (longitudinal) observations on the same "individuals"

  • (on different days, weeks, months, years)

  • Collecting data from a few locations (spatial structure)

  • (surveys conducted in schools/streets, fields/sites/islands)

  • Collecting data on related individuals

  • (father & sons, twins, species within the same genus/tribe/family).

  • If we treat all these observations as fully independent, we are likely to overestimate the number of degrees of freedom.

  • which may lead us to wrongly reject a null hypothesis (type I error)


Non independent observations1
Non-independent Observations

  • Two ways to cope with non-independent observations

  • When design is balanced ("equal sample size")

  • We can use factors to partition our observations in different "groups" and analyse them as an ANOVA or ANCOVA.

  • We already know how to do that (when factors are "crossed")

  • We just need to figure out how to cope with nested factors.

  • When design is unbalanced ("uneven sample size")

  • Mixed effect models are then called for.


Nested anova

control

irrigated

high

medium

Block

low

N

P

NP

Nested Anova

  • Example:

  • A designed field experiment on crop yield with three treatments :

  • irrigation (control, irrigated)

  • sowing density (low, medium, high)

  • fertilizer (N, P, NP)

Split plot design

each block has 18 different subplots


Nested anova1
Nested Anova

  • Example:

  • A designed field experiment on crop yield with three treatments

> yields <- read.table("splityield.txt", header=T)

> attach(yields)

> names(yields)

[1] "yield" "block" "irrigation" "density" "fertilizer"

> str(yields)

'data.frame': 72 obs. of 5 variables:

$ yield : int 90 95 107 92 89 92 81 92 93 80 ...

$ block : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1...

$ irrigation: Factor w/ 2 levels "control","irrigated": 1 1...

$ density : Factor w/ 3 levels "high","low","medium": 2 2...

$ fertilizer: Factor w/ 3 levels "N","NP","P": 1 3 2 1 3 2...


Nested anova2
Nested Anova

  • Example:

  • A designed field experiment on crop yield with three treatments

> model0 <- aov(yield ~ irrigation*density*fertilizer)

## non-nested version (incorrect !!)

> summary(model0)

Df Sum Sq Mean Sq F value Pr(>F)

irrigation 1 8277.6 8277.6 59.5746 2.813e-10 ***

density 2 1758.4 879.2 6.3276 0.0033972 **

fertilizer 2 1977.4 988.7 7.1160 0.0018070 **

irrigation:density 2 2747.0 1373.5 9.8853 0.0002197 ***

irrigation:fertilizer 2 953.4 476.7 3.4310 0.0395615 *

density:fertilizer 4 304.9 76.2 0.5486 0.7008151

irrigation:density:fertilizer 4 234.7 58.7 0.4223 0.7918283

Residuals 54 7503.0 138.9

Sum 71 23756.44


Nested anova3
Nested Anova

> model1 <- aov(yield ~ irrigation*density*fertilizer +

Error(block/irrigation/density) )

## Correct nested version, nesting from large to small

> summary(model1)

Error:block Df Sum Sq Mean Sq F value Pr(>F)

Residuals 3 194.44 64.815

Error:block:irrigation Df Sum Sq Mean Sq F value Pr(>F)

irrigation 1 8277.6 8277.6 17.59 0.0247 *

Residuals 3 1411.8 470.6

Error:block:irrigation:density Df Sum Sq Mean Sq F value Pr(>F)

density 2 1758.4 879.18 3.784 0.05318 .

irrigation:density 2 2747.0 1373.51 5.912 0.01633 *

Residuals 12 2787.9 232.33

Error:Within Df Sum Sq Mean Sq F value Pr(>F)

fertilizer 2 1977.4 988.72 11.449 0.000142 ***

irrigation:fertilizer 2 953.4 476.72 5.520 0.008108 **

density:fertilizer 4 304.9 76.22 0.883 0.484053

irrigation:density:fertilizer 4 234.7 58.68 0.679 0.610667

Residuals 36 3108.8 86.36

Res Sum 54 7503.00Gd Sum 71 23756.44


Nested anova4

control

irrigated

high

medium

Block

low

N

P

NP

Nested Anova

  • Comparison between nested and non-nested results

Non-nested Nested

Df F value Pr(>F) F value Pr(>F)

irrigation 1 59.5746 2.81e-10 17.5896 0.024725

density 2 6.3276 0.003397 3.7842 0.053181

fertilizer 2 7.1160 0.001807 11.4493 0.000142

irrig:dens 2 9.8853 0.000220 5.9119 0.016331

irrig:ferti 2 3.4310 0.039562 5.5204 0.008108 dens:ferti 4 0.5486 0.700815 0.8826 0.484053

irrig:dens:ferti 4 0.4223 0.791828 0.6795 0.610667


Recognizing nestedness is key
Recognizing Nestedness is key !

  • Being able to distinguish crossed factors (independent from each other) from nested factors is essential.

  • Nestedness occurs most often from spatial structure

  • Student surveys in different classes from different schools.

  • Samples from individual branches on sets of trees within a number of forest patches.

  • But can also occur from temporal structure

  • Samples taken from the same individuals every fortnight for 2 months on two successive years.


When the design is not balanced
When the design is not balanced

  • We need a different modelling framework:Mixed Effects Models.

  • So called because they mix together fixed effectsand random effects.

  • Until now, we have only used fixed effects in our models,

  • each effect having an estimated parameter

  • (intercept, slope, mean, ...).

  • But in certain circumstances, these parameters may not be very informative and one would be better off trying to "estimate" the underlying distribution they come from.

  • An example will help clarify the difference between these 2 approaches.


Mixed effects modelling
Mixed Effects Modelling

  • Example: railway rails tested for longitudinal stress.

  • 6 rails chosen at random and tested three times with ultrasound.

> library(nlme) ## package dedicated to mixed effects models

> data(Rail)

> names(Rail)

[1] "Rail" "travel"

> stripchart(Rail$travel ~ Rail$Rail, pch=16,

ylab="Ultrasonic Travel Time (nanosecs)",

xlab="Rail number", vertical =T,

col=rainbow(6) )

> abline(h=mean(Rail$travel),

col="Gray85", lty=2, lwd=2)

Classically in a linear model, we would be able to tell whether the rails differ significantly from each other.

But it doesn't help us make predictions about other rails.


Mixed effects modelling1
Mixed Effects Modelling

  • Random effects: interested in explaining the variance of the response.

  • Fixed effects: interested in explaining the response itself.

Fixed effects

Male & Female

Control & Treatment

Wet vs. Dry

Light vs. Shade

Random effects

Blocks

Individuals withRepeated measures

Genotypes

Sites


Mixed effects modelling2

Makes Rail a simple factor (not an ordered one)

Mixed Effects Modelling

  • Back to our example:

> Rail2 <- data.frame(travel=Rail$travel, Rail=factor(as.character(Rail$Rail)) )

> Rail.lm <- lm(travel ~ Rail, data=Rail2) ## LINEAR MODEL

> summary(Rail.lm)

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 54.000 2.321 23.262 2.37e-11 ***

Rail2 -22.333 3.283 -6.803 1.90e-05 ***

Rail3 30.667 3.283 9.341 7.44e-07 ***

Rail4 42.000 3.283 12.793 2.36e-08 ***

Rail5 -4.000 3.283 -1.218 0.246

Rail6 28.667 3.283 8.732 1.52e-06 ***

---

Residual standard error: 4.021 on 12 degrees of freedom

Multiple R-squared: 0.9796, Adjusted R-squared: 0.9711

F-statistic: 115.2 on 5 and 12 DF, p-value: 1.033e-09


Mixed effects modelling3

Standard Deviation associated with rails

Standard Deviation of residuals

Grand Average

Mixed Effects Modelling

> anova(Rail.lm)

Analysis of Variance Table

Response: travel

Df Sum Sq Mean Sq F value Pr(>F)

Rail 5 9310.5 1862.1 115.18 1.033e-09 ***

Residuals 12 194.0 16.2

## now as a MIXED EFFECT MODEL

> Rail.lme <- lme(travel ~ 1, data=Rail, random= ~1|Rail)

> summary(Rail.lme)

Linear mixed-effects model fit by REML

Data: Rail

AIC BIC logLik

128.177 130.6766 -61.0885

Random effects:

Formula: ~1 | Rail

(Intercept) Residual

StdDev: 24.80547 4.020779

Fixed effects: travel ~ 1

Value Std.Error DF t-value p-value

(Intercept) 66.5 10.17104 12 6.538173 0


Mixed effects modelling4

Standard Deviation associated with rails

Standard Deviation of residuals

Grand Average

Mixed Effects Modelling

> anova(Rail.lm)

Analysis of Variance Table

Response: travel

Df Sum Sq Mean Sq F value Pr(>F)

Rail 5 9310.5 1862.1 115.18 1.033e-09 ***

Residuals 12 194.0 16.2

## now as a MIXED EFFECT MODEL

> Rail.lme <- lme(travel ~ 1, data=Rail, random= ~1|Rail)

> summary(Rail.lme)

Linear mixed-effects model fit by REML

Data: Rail

AIC BIC logLik

128.177 130.6766 -61.0885

Random effects:

Formula: ~1 | Rail

(Intercept) Residual

StdDev: 24.80547 4.020779

Fixed effects: travel ~ 1

Value Std.Error DF t-value p-value

(Intercept) 66.5 10.17104 12 6.538173 0


Mixed effects modelling5
Mixed Effects Modelling

Is the Random effect significant ?

> Rail.lme$call ## random effect model

lme.formula(fixed = travel ~ 1, data = Rail, random = ~1 | Rail)

> AIC(Rail.lme) ## exact AIC version

[1] 128.177

> Rail.lm0 <- lm(travel ~ 1, data=Rail2) ## NULL linear model

> AIC(Rail.lm0)

[1] 167.9265

> Rail.lm$call ## model with Rail as a fixed effect factor

lm(formula = travel ~ Rail, data = Rail)

> AIC(Rail.lm)

[1] 107.8765

Comparing models which only differ in their random effects is easy (with AIC).

Comparing models which differ in their fixed effects is a little harder. Can only be done using "maximum likelihood" (not the default method in lme).


Mixed effects modelling6
Mixed Effects Modelling

Applied on the split plot study of crop yield

> yields <- read.table("splityield.txt", header=T)

> yield.lme <- lme(yield ~ irrigation*density*fertilizer, data=yields,

random= ~1|block/irrigation/density)

> summary(yield.lme)

Linear mixed-effects model fit by REML

Data: yields

AIC BIC logLik

481.6212 525.3789 -218.8106

Random effects:

Formula: ~1 | block

(Intercept)

StdDev: 0.0006600339

Formula: ~1 | irrigation %in% block

(Intercept)

StdDev: 1.982461

Formula: ~1 | density %in% irrigation %in% block

(Intercept) Residual

StdDev: 6.975554 9.292805 ...


Mixed effects modelling7
Mixed Effects Modelling

Fixed effects: yield ~ irrigation * density * fertilizer

Value Std.Error DF t-value p-value

(Intercept) 80.50 5.893741 36 13.658558 0.0000

irrigirrig 31.75 8.335008 3 3.809234 0.0318

dnslow 5.50 8.216282 12 0.669403 0.5159

dnsmed 14.75 8.216282 12 1.795216 0.0978

fertiNP 5.50 6.571005 36 0.837010 0.4081

fertiP 4.50 6.571005 36 0.684827 0.4978

irrigirrig:dnslow -39.00 11.619577 12 -3.356404 0.0057

irrigirrig:dnsmed -22.25 11.619577 12 -1.914872 0.0796

irrigirrig:fertiNP 13.00 9.292805 36 1.398932 0.1704

irrigirrig:fertiP 5.50 9.292805 36 0.591856 0.5576

dnslow:fertiNP 3.25 9.292805 36 0.349733 0.7286

dnsmed:fertiNP -6.75 9.292805 36 -0.726368 0.4723

dnslow:fertiP -5.25 9.292805 36 -0.564953 0.5756

dnsmed:fertiP -5.50 9.292805 36 -0.591856 0.5576

irrigirrig:dnslow:fertiNP 7.75 13.142011 36 0.589712 0.5591

irrigirrig:dnsmed:fertiNP 3.75 13.142011 36 0.285344 0.7770

irrigirrig:dnslow:fertiP 20.00 13.142011 36 1.521837 0.1368

irrigirrig:dnsmed:fertiP 4.00 13.142011 36 0.304367 0.7626


Mixed effects modelling8
Mixed Effects Modelling

> anova(yield.lme)

numDF denDF F-value p-value

(Intercept) 1 36 2674.6301 <.0001

irrigation 1 3 30.9207 0.0115

density 2 12 3.7842 0.0532

fertilizer 2 36 11.4493 0.0001

irrigation:density 2 12 5.9119 0.0163

irrigation:fertilizer 2 36 5.5204 0.0081

density:fertilizer 4 36 0.8826 0.4841

irrigation:density:fertilizer 4 36 0.6795 0.6107

## We should probably remove the three-way interaction

## But if we are fiddling with the fixed effects, we ought ## to fit the model through Maximum Likelihood and base our ## decisions on its AIC values and Likelihood Ratio Tests

> yield.lme.ml <- update(yield.lme, ~. ,method="ML")

> AIC(yield.lme.ml)

[1] 573.5108


Mixed effects modelling9
Mixed Effects Modelling

> yield.lme.ml2 <- update(yield.lme.ml, ~.

- irrigation:density:fertilizer)

> yield.lme.ml2$method

[1] "ML" ## just checking that update() kept using "ML"

> AIC(yield.lme.ml2)

[1] 569.0046 ## an improvement

> anova(yield.lme.ml2)

numDF denDF F-value p-value

(Intercept) 1 40 2872.7394 <.0001

irrigation 1 3 33.2110 0.0104

density 2 12 4.0645 0.0449

fertilizer 2 40 11.4341 0.0001

irrigation:density 2 12 6.3499 0.0132

irrigation:fertilizer 2 40 5.5131 0.0077

density:fertilizer 4 40 0.8815 0.4837

> yield.lme.ml3 <- update(yield.lme.ml2, ~.

- density:fertilizer)

> AIC(yield.mle.lm3)

[1] 565.1933


Mixed effects modelling10
Mixed Effects Modelling

> anova(yield.lme.ml,yield.lme.ml2)

Model df AIC BIC logLik Test L.Ratio p-value

yield.lme.ml 1 22 573.5108 623.5974 -264.7554

yield.lme.ml2 2 18 569.0046 609.9845 -266.5023 1vs2 3.49379 0.4788

> anova(yield.lme.ml2, yield.lme.ml3)

Model df AIC BIC logLik Test L.Ratio p-value

yield.lme.ml2 1 18 569.0046 609.9845 -266.5023

yield.lme.ml3 2 14 565.1933 597.0667 -268.5967 1vs2 4.18877 0.3811

> anova(yield.lme.ml3)

numDF denDF F-value p-value

(Intercept) 1 44 3070.8771 <.0001

irrigation 1 3 35.5016 0.0095

density 2 12 4.3448 0.0381

fertilizer 2 44 11.2013 0.0001

irrigation:density 2 12 6.7878 0.0107

irrigation:fertilizer 2 44 5.4008 0.0080

> yield.lme.ml4 <- update(yield.lme.ml3, ~.

–irrigation:density)

> AIC(yield.mle.ml4)

[1] 572.9022


Mixed effects modelling11

ml4 is one simplification too far

column

matrix

Mixed Effects Modelling

> anova(yield.lme.m3,yield.lme.ml4)

Model df AIC BIC logLik Test L.Ratio p-value

yield.lme.ml3 1 14 565.1933 597.0667 -268.5967

yield.lme.ml4 2 12 572.9022 600.2221 -274.4511 1vs2 11.7088 0.0029

> anova(yield.lme.ml4)

numDF denDF F-value p-value

(Intercept) 1 44 2138.9678 <.0001

irrigation 1 3 24.7281 0.0156

density 2 14 2.6264 0.1075

fertilizer 2 44 11.5626 0.0001

irrigation:fertilizer 2 44 5.5750 0.0069

## here comes Model Checking

> shapiro.test(yield.lme.ml3$residuals[,"fixed"]) # Best Model

Shapiro-Wilk normality test

data: yield.lme.ml3$residuals[, "fixed"]

W = 0.9797, p-value = 0.2958


Mixed effects modelling12

including all random effects

excluding all random effects

Mixed Effects Modelling

> res <- yield.lme.ml3$resid[,"fixed"]

> st.res <- res/sd(res)

> qqnorm(st.res, pch=16, main="")

> qqline(st.res, col="red", lwd=2)

> res <- yield.lme.ml3$resid[,4]

> st.res <- res/sd(res)

> qqnorm(st.res, pch=16, main="")

> qqline(st.res, col="red", lwd=2)


Mixed effects modelling13
Mixed Effects Modelling

> plot(yield.lme.ml3)

## by default Residuals vs Fitted values

> plot(yield.lme.ml3, yield ~ fitted(.) )

## Observed vs Fitted values


Mixed effects modelling14
Mixed Effects Modelling

> qqnorm(yield.lme.ml3, ~resid(.) | block)

## qqplot but broken down by block


ad