Solutions to exercise 3 Rasmus Waagepetersen October 16, 2019 - - PowerPoint PPT Presentation

Solutions to exercise 3

Rasmus Waagepetersen October 16, 2019

Exercise 3 (rats)

We plot response versus log-age (and age) rats$logage=log(rats$age) library(lattice) xyplot(response~logage,group=rat,type="l",data=rats) Vs log-age Vs age

logage response 70 75 80 85 0.0 0.5 1.0 1.5 2.0 age response 70 75 80 85 50 60 70 80 90 100 110

There seems to be a linear dependence between response and log age but with rat specific intercepts.

We fit an ordinary linear regression with interaction between log age and treatment (this means that we have different slopes for different treatments). We conduct an F-test for no interaction.

> fit=lm(response~logage*treat,data=rats) > summary(fit) ... Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 71.82847 0.48178 149.089 <2e-16 *** logage 5.56738 0.41897 13.288 <2e-16 *** treathig

• 1.32020

0.65171

• 2.026

0.0439 * treatlow

• 0.02737

0.64526

• 0.042

0.9662 logage:treathig 0.27832 0.54319 0.512 0.6088 logage:treatlow 0.58649 0.54161 1.083 0.2799 ... Residual standard error: 2.154 on 246 degrees of freedom ... > fit2=lm(response~logage+treat,data=rats) > anova(fit,fit2) Analysis of Variance Table Model 1: response ~ logage * treat Model 2: response ~ logage + treat Res.Df RSS Df Sum of Sq F Pr(>F) 1 246 1141.0 2 248 1146.6 -2

• 5.5678 0.6002 0.5495

According to this analysis, we accept equal slopes across treatments.

If we continue to do an F-test for treatment we obtain a p-value of 5.433e-6 (thus there seems to be a treatment effect):

> fit3=lm(response~logage,data=rats) > anova(fit2,fit3) Analysis of Variance Table Model 1: response ~ logage + treat Model 2: response ~ logage Res.Df RSS Df Sum of Sq F Pr(>F) 1 248 1146.6 2 250 1271.7 -2

• 125.08 13.527 2.657e-06 ***

However, all this assumed equal intercepts among rats - which does not seem correct according to the initial plot.

We now try a mixed-model analysis with random rat specific intercepts:

> library(lme4) > library(pbkrtest) > fitm=lmer(response~logage*treat+(1|rat),data=rats) > summary(fitm) ... Random effects: Groups Name Variance Std.Dev. rat (Intercept) 3.512 1.874 Residual 1.378 1.174 Number of obs: 252, groups: rat, 50 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 71.6755 0.5523 64.7452 129.787 <2e-16 *** logage 5.8878 0.2509 209.9808 23.465 <2e-16 *** treathig

• 1.1508

0.7548 63.8359

• 1.525

0.132 treatlow 0.1223 0.7454 63.9991 0.164 0.870 logage:treathig

• 0.0799

0.3224 208.8490

• 0.248

0.805 logage:treatlow 0.2314 0.3206 208.7366 0.722 0.471

The residual variance 2.1542 = 4.64 from the previous model is split into variance 3.5 of random rat specific intercepts and (new) residual variance 1.4. There seems to be quite a strong between rats variance (the proportion of the total variance is 0.72=3.512/(3.512+1.378)).

We again compute an F-test for no interaction and subsequently an F-test for treatment:

> fitm2=lmer(response~logage+treat+(1|rat),data=rats) > KRmodcomp(fitm,fitm2) F-test with Kenward-Roger approximation; computing time: 0.52 sec. large : response ~ logage * treat + (1 | rat) small : response ~ logage + treat + (1 | rat) stat ndf ddf F.scaling p.value Ftest 0.6349 2.0000 208.2765 1 0.531 > fitm3=lmer(response~logage+(1|rat),data=rats) > KRmodcomp(fitm2,fitm3) F-test with Kenward-Roger approximation; computing time: 0.22 sec. large : response ~ logage + treat + (1 | rat) small : response ~ logage + (1 | rat) stat ndf ddf F.scaling p.value Ftest 3.0538 2.0000 46.8389 1 0.05667

The interaction is not significant. Moreover, according to the mixed model analysis, neither is treatment ! (at the 5% level - p-value is 0.057)

Model check

We check the initial mixed model using residuals and BLUPS of random effects.

> ratsres=resid(fitm) > ratsfitted=fitted(fitm) > plot(ratsfitted,ratsres) > boxplot(ratsres~rats$rat) > qqnorm(ratsres) > qqline(ratsres) > us=ranef(fitm)$rat > qqnorm(us[[1]]) > qqline(us[[1]])

• resid. vs. fitted

boxplots

• 70

75 80 85 −2 −1 1 2 ratsfitted rnorm(length(ratsfitted))

• 1

6 9 12 17 24 28 33 38 43 49 53 59 64 −2 −1 1 2 3 rats$rat ratsres

qq-plot for resid. qq-plot for BLUPS (random intcpts.)

• ●
• −3

−2 −1 1 2 3 −2 −1 1 2 3 Normal Q−Q Plot Theoretical Quantiles Sample Quantiles

• ●
• −2

−1 1 2 −4 −2 2 4 Normal Q−Q Plot Theoretical Quantiles Sample Quantiles

Not obvious pattern in plot of residuals versus fitted/rat number. Beautiful qq-plot for residuals. Maybe one outlying rat according to qq-plot for predictions (BLUPs) of random intercepts. Overall, the model appears satisfactory.