Mixed Models
A Student
2022-03-30
CHAPTER 10 OUTPUT
This shows the outputs from Chapter 10 using R. The textbook is “A Modern Approach to Regression with R” by Simon J. Sheather (2008). The R code provided with the book has been updated.
library(nlme)
Orthodont <- read_excel("MARData.xlsx", sheet = "Orthodont")
FOrthodont <- Orthodont[Orthodont$Sex == "Female",]
#Figure 10.1 on page 333
plot(FOrthodont)
#Output from R on page 334
DistFAge8 <- FOrthodont$distance[FOrthodont$age == 8]
DistFAge10 <- FOrthodont$distance[FOrthodont$age == 10]
DistFAge12 <- FOrthodont$distance[FOrthodont$age == 12]
DistFAge14 <- FOrthodont$distance[FOrthodont$age == 14]
T <- cbind(DistFAge8,DistFAge10,DistFAge12,DistFAge14)
c<-cor(T)
round(c,3)## DistFAge8 DistFAge10 DistFAge12 DistFAge14
## DistFAge8 1.000 0.830 0.862 0.841
## DistFAge10 0.830 1.000 0.895 0.879
## DistFAge12 0.862 0.895 1.000 0.948
## DistFAge14 0.841 0.879 0.948 1.000#Figure 10.2 on page 335
pairs(~DistFAge14+DistFAge12+DistFAge10+DistFAge8,lower.panel = NULL)
#Output from R on page 337
mFRI <- lme(distance~age,data = FOrthodont,random = ~1|Subject,method = "REML")
summary(mFRI)## Linear mixed-effects model fit by REML
## Data: FOrthodont
## AIC BIC logLik
## 149.2183 156.169 -70.60916
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 2.06847 0.7800331
##
## Fixed effects: distance ~ age
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 0.8587419 32 20.230440 0
## age 0.479545 0.0525898 32 9.118598 0
## Correlation:
## (Intr)
## age -0.674
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2736479 -0.7090164 0.1728237 0.4122128 1.6325181
##
## Number of Observations: 44
## Number of Groups: 11#Figure 10.3 on page 338
plot(mFRI, form = distance ~ fitted(.) | Subject,layout = c(3,4), between = list(y = c(0, 0, 0, 0.5)),
aspect = 1.0, abline = c(0,1))
#Fixed intercepts in Table 10.2 on page 338
mFFI <- lm(distance~factor(Subject)+age - 1,data = FOrthodont)
summary(mFFI)##
## Call:
## lm(formula = distance ~ factor(Subject) + age - 1, data = FOrthodont)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.81136 -0.50000 0.08409 0.36477 1.14545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## factor(Subject)F01 16.10000 0.69768 23.076 < 2e-16 ***
## factor(Subject)F02 17.72500 0.69768 25.406 < 2e-16 ***
## factor(Subject)F03 18.47500 0.69768 26.481 < 2e-16 ***
## factor(Subject)F04 19.60000 0.69768 28.093 < 2e-16 ***
## factor(Subject)F05 17.35000 0.69768 24.868 < 2e-16 ***
## factor(Subject)F06 15.85000 0.69768 22.718 < 2e-16 ***
## factor(Subject)F07 17.72500 0.69768 25.406 < 2e-16 ***
## factor(Subject)F08 18.10000 0.69768 25.943 < 2e-16 ***
## factor(Subject)F09 15.85000 0.69768 22.718 < 2e-16 ***
## factor(Subject)F10 13.22500 0.69768 18.956 < 2e-16 ***
## factor(Subject)F11 21.10000 0.69768 30.243 < 2e-16 ***
## age 0.47955 0.05259 9.119 2.06e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.78 on 32 degrees of freedom
## Multiple R-squared: 0.9991, Adjusted R-squared: 0.9988
## F-statistic: 3122 on 12 and 32 DF, p-value: < 2.2e-16#Random intercepts in Table 10.2 on page 338
coef(mFRI)## (Intercept) age
## F01 16.14369 0.4795455
## F02 17.71291 0.4795455
## F03 18.43716 0.4795455
## F04 19.52353 0.4795455
## F05 17.35078 0.4795455
## F06 15.90228 0.4795455
## F07 17.71291 0.4795455
## F08 18.07503 0.4795455
## F09 15.90228 0.4795455
## F10 13.36740 0.4795455
## F11 20.97204 0.4795455#Figure 10.4 on page 339
MOrthodont <- Orthodont[Orthodont$Sex == "Male",]
plot(MOrthodont)
#Output from R on pages 340 and 341
DistMAge8 <- MOrthodont$distance[MOrthodont$age == 8]
DistMAge10 <- MOrthodont$distance[MOrthodont$age == 10]
DistMAge12 <- MOrthodont$distance[MOrthodont$age == 12]
DistMAge14 <- MOrthodont$distance[MOrthodont$age == 14]
T <- cbind(DistMAge8,DistMAge10,DistMAge12,DistMAge14)
c<-cor(T)
round(c,3)## DistMAge8 DistMAge10 DistMAge12 DistMAge14
## DistMAge8 1.000 0.437 0.558 0.315
## DistMAge10 0.437 1.000 0.387 0.631
## DistMAge12 0.558 0.387 1.000 0.586
## DistMAge14 0.315 0.631 0.586 1.000#Figure 10.5 on page 340
pairs(~DistMAge14+DistMAge12+DistMAge10+DistMAge8,lower.panel = NULL)
#Output from R on page 341
mMRI <- lme(distance~age,data = MOrthodont,random = ~1|Subject,method = "REML")
summary(mMRI)## Linear mixed-effects model fit by REML
## Data: MOrthodont
## AIC BIC logLik
## 281.448 289.9566 -136.724
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.625019 1.67822
##
## Fixed effects: distance ~ age
## Value Std.Error DF t-value p-value
## (Intercept) 16.340625 1.1287202 47 14.477126 0
## age 0.784375 0.0938154 47 8.360838 0
## Correlation:
## (Intr)
## age -0.914
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.008054554 -0.640688586 0.007833248 0.534480581 3.052946887
##
## Number of Observations: 64
## Number of Groups: 16#Figure 10.6 on page 342
plot(mMRI, form = distance ~ fitted(.) | Subject,layout = c(4,4), between = list(y = c(0, 0, 0, 0.5)),
aspect = 1.0, abline = c(0,1))
#Output from R on pages 343 and 344
m10.5 <- lme(distance~age*Sex,data = Orthodont,random = ~1|Subject,method = "REML",
weights = varIdent(form = ~1|Sex))
summary(m10.5)## Linear mixed-effects model fit by REML
## Data: Orthodont
## AIC BIC logLik
## 429.2205 447.7312 -207.6102
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.847571 0.7813007
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Sex
## Parameter estimates:
## Female Male
## 1.000000 2.137235
## Fixed effects: distance ~ age * Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 0.8123609 79 21.385479 0.0000
## age 0.479545 0.0526753 79 9.103804 0.0000
## SexMale -1.032102 1.4039842 25 -0.735124 0.4691
## age:SexMale 0.304830 0.1071828 79 2.844016 0.0057
## Correlation:
## (Intr) age SexMal
## age -0.713
## SexMale -0.579 0.413
## age:SexMale 0.351 -0.491 -0.840
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.00556473 -0.63419466 0.01890458 0.55016878 3.06446973
##
## Number of Observations: 108
## Number of Groups: 27#Output from R on pages 344 and 345
m10.6 <- lme(distance~age*Sex,data = Orthodont,random = ~1|Subject,method = "REML")
summary(m10.6)## Linear mixed-effects model fit by REML
## Data: Orthodont
## AIC BIC logLik
## 445.7572 461.6236 -216.8786
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.816214 1.386382
##
## Fixed effects: distance ~ age * Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 1.1835071 79 14.679023 0.0000
## age 0.479545 0.0934698 79 5.130483 0.0000
## SexMale -1.032102 1.5374208 25 -0.671321 0.5082
## age:SexMale 0.304830 0.1214209 79 2.510520 0.0141
## Correlation:
## (Intr) age SexMal
## age -0.869
## SexMale -0.770 0.669
## age:SexMale 0.669 -0.770 -0.869
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.59804400 -0.45461690 0.01578365 0.50244658 3.68620792
##
## Number of Observations: 108
## Number of Groups: 27#Output from R on page 345
anova(m10.6,m10.5)## Model df AIC BIC logLik Test L.Ratio p-value
## m10.6 1 6 445.7572 461.6236 -216.8786
## m10.5 2 7 429.2205 447.7312 -207.6102 1 vs 2 18.53677 <.0001#Figure 10.7 on page 347
qqnorm(m10.5,~ranef(.),id = 0.05,cex = 0.7)
#Output from R on page 347
MRes <- resid(m10.5, level = 0)
MRAge8 <- MRes[Orthodont$age == 8]
MRAge10 <- MRes[Orthodont$age == 10]
MRAge12 <- MRes[Orthodont$age == 12]
MRAge14 <- MRes[Orthodont$age == 14]
T <- cbind(MRAge8,MRAge10,MRAge12,MRAge14)
c<-cor(T)
round(c,3)## MRAge8 MRAge10 MRAge12 MRAge14
## MRAge8 1.000 0.560 0.660 0.522
## MRAge10 0.560 1.000 0.560 0.718
## MRAge12 0.660 0.560 1.000 0.728
## MRAge14 0.522 0.718 0.728 1.000#Figure 10.8 on page 348
pairs(~MRAge14+MRAge12+MRAge10+MRAge8,lower.panel = NULL)
#Output from R on page 348
CRes <- resid(m10.5, level = 1)
CRAge8 <- CRes[Orthodont$age == 8]
CRAge10 <- CRes[Orthodont$age == 10]
CRAge12 <- CRes[Orthodont$age == 12]
CRAge14 <- CRes[Orthodont$age == 14]
T <- cbind(CRAge8,CRAge10,CRAge12,CRAge14)
c<-cor(T)
round(c,3)## CRAge8 CRAge10 CRAge12 CRAge14
## CRAge8 1.000 -0.307 -0.120 -0.620
## CRAge10 -0.307 1.000 -0.544 -0.005
## CRAge12 -0.120 -0.544 1.000 -0.083
## CRAge14 -0.620 -0.005 -0.083 1.000#Figure 10.9 on page 349
pairs(~CRAge14+CRAge12+CRAge10+CRAge8,lower.panel = NULL)
#Figure 10.10 on page 350
plot(m10.6,form = resid(., type = "response",level = 1) ~ fitted(.,level = 1) | Sex,abline = 0)
#Figure 10.11 on page 351
plot(m10.5,form = resid(., type = "response",level = 1) ~ fitted(.,level = 1) | Sex,abline = 0)
#Figure 10.12 on page 352
#Choleski Residuals for m10.5
attach(Orthodont)
m10.5.a <- lm(distance~age*Sex, data = Orthodont)
m10.5.b <- lm(distance~(Subject-1), data = Orthodont)
m10.5.X <- model.matrix(m10.5.a)
m10.5.Z <- model.matrix(m10.5.b)
m10.5.G <- diag(rep(getVarCov(m10.5),ncol(m10.5.Z)))
m10.5.R <- diag(attr(m10.5[[15]],"std")^2)
m10.5.V <- m10.5.Z %*% m10.5.G %*% t(m10.5.Z) + m10.5.R
m10.5.tCHOLinv <- solve(t(chol(m10.5.V)))
# Now to premultiply terms by m10.5.tCHOLinv to get Choleski residuals, etc.
dist.CHOL.m10.5<- m10.5.tCHOLinv %*% Orthodont$distance
x1.CHOL.m10.5 <- m10.5.tCHOLinv %*% m10.5.X[,1]
x2.CHOL.m10.5 <- m10.5.tCHOLinv %*% m10.5.X[,2]
x3.CHOL.m10.5 <- m10.5.tCHOLinv %*% m10.5.X[,3]
x4.CHOL.m10.5 <- m10.5.tCHOLinv %*% m10.5.X[,4]
m10.5.CHOL <- lm(dist.CHOL.m10.5 ~ x1.CHOL.m10.5 + x2.CHOL.m10.5 +
x3.CHOL.m10.5 + x4.CHOL.m10.5 -1)
#summary(m10.5.CHOL)
#Choleski Residuals for m10.6
m10.6.a <- lm(distance~age*Sex, data = Orthodont)
m10.6.b <- lm(distance~(Subject-1), data = Orthodont)
m10.6.X <- model.matrix(m10.6.a)
m10.6.Z <- model.matrix(m10.6.b)
m10.6.G <- diag(rep(getVarCov(m10.6),ncol(m10.6.Z)))
m10.6.R <- diag(attr(m10.6[[15]],"std")^2)
m10.6.V <- m10.6.Z %*% m10.6.G %*% t(m10.6.Z) + m10.6.R
m10.6.tCHOLinv <- solve(t(chol(m10.6.V)))
# Now to premultiply terms by m10.6.tCHOLinv to get Choleski residuals, etc.
dist.CHOL.m10.6<- m10.6.tCHOLinv %*% Orthodont$distance
x1.CHOL.m10.6 <- m10.6.tCHOLinv %*% m10.6.X[,1]
x2.CHOL.m10.6 <- m10.6.tCHOLinv %*% m10.6.X[,2]
x3.CHOL.m10.6 <- m10.6.tCHOLinv %*% m10.6.X[,3]
x4.CHOL.m10.6 <- m10.6.tCHOLinv %*% m10.6.X[,4]
m10.6.CHOL <- lm(dist.CHOL.m10.6 ~ x1.CHOL.m10.6 + x2.CHOL.m10.6 +
x3.CHOL.m10.6 + x4.CHOL.m10.6 -1)
#summary(m10.6.CHOL)
CholeskyResid10.5 <- m10.5.CHOL$residuals
CholeskyResid10.6 <- m10.6.CHOL$residuals
#
# par(mfrow = c(1,2))
# plot(CholeskyResid10.6~Sex,ylab = "Cholesky residuals from model (10.6)",ylim = c(-3,4))
# plot(CholeskyResid10.5~Sex,ylab = "Cholesky residuals from model (10.5)",ylim = c(-3,4))
#p-value for Levene's test on page 353
library(car)
leveneTest(CholeskyResid10.6,Sex)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 6.535 0.01199 *
## 106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1detach(Orthodont)library(foreign)
library(nlme)
pigweights <- read_excel("MARData.xlsx", sheet = "pigweights")
pigweights$Animal <- pigweights$pigid
pigweights$Time <- pigweights$weeknumber
attach(pigweights)
pigweights <- groupedData(weight ~ Time | Animal,data = pigweights,
labels = list(x = "Time",y = "Weight"),
units = list(x = "(weeks)", y = "(kg)"))
#Figure 10.13 on page 354
par(mfrow = c(1,1))
plot(Time,weight)
lines(Time,weight,lty = 2)
#Output from R on page 354
T1 <- weight[Time == 1]
T2 <- weight[Time == 2]
T3 <- weight[Time == 3]
T4 <- weight[Time == 4]
T5 <- weight[Time == 5]
T6 <- weight[Time == 6]
T7 <- weight[Time == 7]
T8 <- weight[Time == 8]
T9 <- weight[Time == 9]
T <- cbind(T1,T2,T3,T4,T5,T6,T7,T8,T9)
c<-cor(T)
round(c,2)## T1 T2 T3 T4 T5 T6 T7 T8 T9
## T1 1.00 0.92 0.80 0.80 0.75 0.71 0.66 0.63 0.56
## T2 0.92 1.00 0.91 0.91 0.88 0.84 0.78 0.71 0.66
## T3 0.80 0.91 1.00 0.96 0.93 0.91 0.84 0.82 0.77
## T4 0.80 0.91 0.96 1.00 0.96 0.93 0.87 0.83 0.79
## T5 0.75 0.88 0.93 0.96 1.00 0.92 0.85 0.81 0.79
## T6 0.71 0.84 0.91 0.93 0.92 1.00 0.96 0.93 0.89
## T7 0.66 0.78 0.84 0.87 0.85 0.96 1.00 0.96 0.92
## T8 0.63 0.71 0.82 0.83 0.81 0.93 0.96 1.00 0.97
## T9 0.56 0.66 0.77 0.79 0.79 0.89 0.92 0.97 1.00#Figure 10.14 on page 355
pairs(~T9+T8+T7+T6+T5+T4+T3+T2+T1,lower.panel = NULL)
#Output from R on page 356 and 357
m10p13 <- gls(weight~Time,data = pigweights,
correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "REML")
summary(m10p13)## Generalized least squares fit by REML
## Model: weight ~ Time
## Data: pigweights
## AIC BIC logLik
## 1635.943 1826.941 -770.9717
##
## Correlation Structure: General
## Formula: ~1 | Animal
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6 7 8
## 2 0.894
## 3 0.722 0.895
## 4 0.767 0.910 0.946
## 5 0.742 0.879 0.889 0.956
## 6 0.694 0.836 0.876 0.930 0.923
## 7 0.650 0.774 0.806 0.863 0.857 0.963
## 8 0.602 0.720 0.812 0.835 0.810 0.927 0.953
## 9 0.546 0.669 0.753 0.789 0.788 0.891 0.917 0.968
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time
## Parameter estimates:
## 1 2 3 4 5 6 7 8
## 1.000000 1.133914 1.514697 1.524492 1.825743 1.798020 2.005045 2.212076
## 9
## 2.561563
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 19.072858 0.3292663 57.92533 0
## Time 6.174365 0.0791561 78.00235 0
##
## Correlation:
## (Intr)
## Time 0.009
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.53084894 -0.50616954 0.06082772 0.72617499 2.77318390
##
## Residual standard error: 2.476839
## Degrees of freedom: 432 total; 430 residual#Choleski Residuals for m10.13
m10p13.V <- getVarCov(m10p13)
m10p13.tCHOLinv <- diag(1,length(unique(pigweights$Animal))) %x% solve(t(chol(m10p13.V)))
ystar <- m10p13.tCHOLinv %*% pigweights$weight
Intstar <- m10p13.tCHOLinv %*% rep(1,length(pigweights$Time))
Timestar <- m10p13.tCHOLinv %*% pigweights$Time
m10p13star <- lm(ystar~Intstar+Timestar-1)
summary(m10p13star)##
## Call:
## lm(formula = ystar ~ Intstar + Timestar - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.06004 -0.59531 -0.03663 0.67704 3.04806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Intstar 19.07286 0.32927 57.92 <2e-16 ***
## Timestar 6.17437 0.07916 78.00 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1 on 430 degrees of freedom
## Multiple R-squared: 0.9561, Adjusted R-squared: 0.9559
## F-statistic: 4678 on 2 and 430 DF, p-value: < 2.2e-16Resm10p13star <- m10p13star$residuals
#Figure 10.15 on page 358
par(mfrow = c(1,1))
plot(Timestar,Resm10p13star,ylab = "Cholesky Residuals",xlab = "x*",ylim = c(-3.5,3.5))
lines(lowess(Timestar,Resm10p13star,iter = 1,f = 1/3),lty = 2)
#F-statistic p-value on page 357
mres <- lm(Resm10p13star~Timestar+I(Timestar^2)+I(Timestar^3)+I(Timestar^4)+I(Timestar^5))
summary(mres)##
## Call:
## lm(formula = Resm10p13star ~ Timestar + I(Timestar^2) + I(Timestar^3) +
## I(Timestar^4) + I(Timestar^5))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.01093 -0.57400 -0.04689 0.66439 2.75290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1224 0.6014 3.529 0.000463 ***
## Timestar -12.2288 3.4030 -3.594 0.000364 ***
## I(Timestar^2) -1.2916 1.7644 -0.732 0.464525
## I(Timestar^3) 87.4975 23.8325 3.671 0.000272 ***
## I(Timestar^4) -134.0208 37.8657 -3.539 0.000445 ***
## I(Timestar^5) 57.4926 17.0288 3.376 0.000802 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9771 on 426 degrees of freedom
## Multiple R-squared: 0.05406, Adjusted R-squared: 0.04296
## F-statistic: 4.869 on 5 and 426 DF, p-value: 0.0002414#Output from R on pages 360 and 361: REML fit of model (10.15)
m10p15 <- gls(weight~Time+TimeM2Plus+TimeM3Plus+TimeM4Plus+TimeM5Plus+TimeM6Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "REML")
summary(m10p15)## Generalized least squares fit by REML
## Model: weight ~ Time + TimeM2Plus + TimeM3Plus + TimeM4Plus + TimeM5Plus + TimeM6Plus + TimeM7Plus + TimeM8Plus
## Data: pigweights
## AIC BIC logLik
## 1613.634 1832.192 -752.817
##
## Correlation Structure: General
## Formula: ~1 | Animal
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6 7 8
## 2 0.916
## 3 0.802 0.912
## 4 0.796 0.908 0.958
## 5 0.749 0.881 0.928 0.962
## 6 0.705 0.835 0.906 0.933 0.922
## 7 0.655 0.776 0.843 0.868 0.855 0.963
## 8 0.625 0.713 0.817 0.829 0.810 0.928 0.959
## 9 0.558 0.664 0.769 0.786 0.786 0.889 0.917 0.969
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time
## Parameter estimates:
## 1 2 3 4 5 6 7 8
## 1.000000 1.130223 1.435540 1.512623 1.836839 1.802342 2.014341 2.197064
## 9
## 2.566117
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 18.260417 0.3801324 48.03699 0.0000
## Time 6.760417 0.1623942 41.62967 0.0000
## TimeM2Plus 0.322917 0.2294243 1.40751 0.1600
## TimeM3Plus -1.552083 0.2925824 -5.30477 0.0000
## TimeM4Plus 0.229167 0.2432420 0.94213 0.3467
## TimeM5Plus 0.531250 0.3931492 1.35127 0.1773
## TimeM6Plus -0.281250 0.2646040 -1.06291 0.2884
## TimeM7Plus 0.833333 0.2974985 2.80114 0.0053
## TimeM8Plus -0.927083 0.2757118 -3.36251 0.0008
##
## Correlation:
## (Intr) Time TmM2Pl TmM3Pl TmM4Pl TmM5Pl TmM6Pl TmM7Pl
## Time -0.355
## TimeM2Plus 0.295 -0.414
## TimeM3Plus 0.004 -0.176 -0.695
## TimeM4Plus 0.051 0.265 0.141 -0.559
## TimeM5Plus -0.131 -0.247 0.081 0.023 -0.613
## TimeM6Plus 0.148 0.111 -0.065 0.045 0.326 -0.706
## TimeM7Plus 0.081 -0.174 0.312 -0.197 0.042 -0.066 -0.377
## TimeM8Plus -0.173 0.414 -0.353 0.147 0.147 -0.239 0.200 -0.482
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.57036756 -0.63836354 -0.04507358 0.59998227 2.72299338
##
## Residual standard error: 2.468838
## Degrees of freedom: 432 total; 423 residual#Output from R on page 361: Comparing ML fits of models (10.13) and (10.15)
m10p15.ML <- gls(weight~Time+TimeM2Plus+TimeM3Plus+TimeM4Plus+TimeM5Plus+TimeM6Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "ML")
m10p13.ML <- gls(weight~Time,data = pigweights,
correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "ML")
anova(m10p13.ML,m10p15.ML)## Model df AIC BIC logLik Test L.Ratio p-value
## m10p13.ML 1 47 1632.303 1823.520 -769.1517
## m10p15.ML 2 54 1600.992 1820.687 -746.4958 1 vs 2 45.31183 <.0001#Output from R on page 363: Comparing ML fits of models (10.15) and (10.16)
m10p16.ML <- gls(weight~Time+TimeM3Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "ML")
anova(m10p16.ML,m10p15.ML)## Model df AIC BIC logLik Test L.Ratio p-value
## m10p16.ML 1 50 1600.125 1803.546 -750.0624
## m10p15.ML 2 54 1600.992 1820.687 -746.4958 1 vs 2 7.133239 0.129#Output from R on page 363: REML fit of model (10.16)
m10p16 <- gls(weight~Time+TimeM3Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "REML")
summary(m10p16)## Generalized least squares fit by REML
## Model: weight ~ Time + TimeM3Plus + TimeM7Plus + TimeM8Plus
## Data: pigweights
## AIC BIC logLik
## 1608.513 1811.352 -754.2563
##
## Correlation Structure: General
## Formula: ~1 | Animal
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6 7 8
## 2 0.916
## 3 0.800 0.909
## 4 0.796 0.909 0.955
## 5 0.749 0.881 0.924 0.963
## 6 0.703 0.832 0.906 0.929 0.918
## 7 0.651 0.771 0.844 0.863 0.849 0.964
## 8 0.620 0.706 0.816 0.821 0.802 0.927 0.959
## 9 0.553 0.657 0.769 0.779 0.778 0.889 0.917 0.970
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time
## Parameter estimates:
## 1 2 3 4 5 6 7 8
## 1.000000 1.130602 1.434683 1.512722 1.835001 1.803820 2.019456 2.206958
## 9
## 2.575040
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 18.256909 0.3550365 51.42262 0.0000
## Time 6.820640 0.1377674 49.50839 0.0000
## TimeM3Plus -0.981702 0.1382102 -7.10297 0.0000
## TimeM7Plus 0.828703 0.2106230 3.93453 0.0001
## TimeM8Plus -0.767969 0.2497387 -3.07509 0.0022
##
## Correlation:
## (Intr) Time TmM3Pl TmM7Pl
## Time -0.298
## TimeM3Plus 0.364 -0.787
## TimeM7Plus 0.005 -0.123 -0.091
## TimeM8Plus -0.113 0.267 -0.082 -0.542
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.62768833 -0.62997203 -0.03143121 0.60053129 2.79484335
##
## Residual standard error: 2.467269
## Degrees of freedom: 432 total; 427 residual#Output from R on page 363: Comparing ML fits of models (10.16) and (10.17)
m10p17.ML <- gls(weight~Time+TimeM3Plus+TimeM5Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corSymm(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "ML")
anova(m10p17.ML,m10p16.ML)## Model df AIC BIC logLik Test L.Ratio p-value
## m10p17.ML 1 51 1599.149 1806.639 -748.5744
## m10p16.ML 2 50 1600.125 1803.546 -750.0624 1 vs 2 2.976025 0.0845#Output from R on page 365: REML fit of model (10.16) with AR errors
m10p16.AR1 <- gls(weight~Time+TimeM3Plus+TimeM7Plus+TimeM8Plus,
data = pigweights,correlation = corAR1(form = ~1|Animal),weights = varIdent(form = ~1|Time),method = "REML")
summary(m10p16.AR1)## Generalized least squares fit by REML
## Model: weight ~ Time + TimeM3Plus + TimeM7Plus + TimeM8Plus
## Data: pigweights
## AIC BIC logLik
## 1612.402 1673.254 -791.2011
##
## Correlation Structure: AR(1)
## Formula: ~1 | Animal
## Parameter estimate(s):
## Phi
## 0.9459712
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time
## Parameter estimates:
## 1 2 3 4 5 6 7 8
## 1.000000 1.129122 1.370541 1.378917 1.644524 1.563853 1.678025 1.788617
## 9
## 2.080040
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 18.124615 0.3585964 50.54321 0.0000
## Time 6.901604 0.1213124 56.89117 0.0000
## TimeM3Plus -1.011724 0.1608420 -6.29017 0.0000
## TimeM7Plus 0.953293 0.2509486 3.79876 0.0002
## TimeM8Plus -0.933510 0.3453570 -2.70303 0.0071
##
## Correlation:
## (Intr) Time TmM3Pl TmM7Pl
## Time -0.085
## TimeM3Plus 0.030 -0.788
## TimeM7Plus 0.032 0.032 -0.260
## TimeM8Plus 0.082 0.084 -0.004 -0.592
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.65953926 -0.62183173 -0.07058175 0.56220717 2.92971160
##
## Residual standard error: 2.768652
## Degrees of freedom: 432 total; 427 residual#Output from R on page 365: Comparing ML fits of models (10.16) with different errors
anova(m10p16.AR1,m10p16)## Model df AIC BIC logLik Test L.Ratio p-value
## m10p16.AR1 1 15 1612.402 1673.254 -791.2011
## m10p16 2 50 1608.513 1811.352 -754.2563 1 vs 2 73.88972 1e-04detach(pigweights)