- 224 Views
- Uploaded on

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

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

- 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

- 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

- 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 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.

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 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 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 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

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 !

- 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

- 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

- 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 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

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

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

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 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 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 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 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 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 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

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

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 Modelling

> plot(yield.lme.ml3)

## by default Residuals vs Fitted values

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

## Observed vs Fitted values

Mixed Effects Modelling

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

## qqplot but broken down by block

Download Presentation

Connecting to Server..