# This file posted as weight2.ssc
#
# This code is applied to the weighting
# data from Littel, et. al. 1991.
temp <- read.table("c:/courses/st511/snew/weight2.dat",
col.names=c("Subj","Program","S1",
"S2","S3","S4","S5","S6","S7"))
temp
# Create a data file where responses at different
# time points are on different lines
attach(temp)
weight <- data.frame(rbind(
cbind(Subj,Program, 2,S1),
cbind(Subj,Program, 4,S2),
cbind(Subj,Program, 6,S3),
cbind(Subj,Program, 8,S4),
cbind(Subj,Program,10,S5),
cbind(Subj,Program,12,S6),
cbind(Subj,Program,Time=14,
Strength=S7)))
detach( )
# Create factors
weight$Subj <- as.factor(weight$Subj)
weight$Program <- as.factor(weight$Program)
weight$Timef <- as.factor(weight$Time)
# Sort the data set by suject
i <- order(weight$Subj,weight$Time)
weight <- weight[i,]
# Delete the list called i
rm(i)
# Compute sample means
means <- tapply(weight$Strength,
list(weight$Time,weight$Program),mean)
means
# Make a profile plot of the means
# Unix users should insert the motif( )
# command
# motif()
x.axis <- unique(weight$Time)
par(fin=c(6.0,6.0),pch=18,mkh=.1,mex=1.5,
cex=1.2,lwd=3)
matplot(c(2,14), c(79,85.7), type="n",
xlab="Time(Days)", ylab="Strength",
main= "Observed Strength Means")
matlines(x.axis,means,type='l',lty=c(1,3,7))
matpoints(x.axis,means, pch=c(16,17,15))
legend(2.1,85.69,legend=c("RI program",
'WI Program','Controls'),lty=c(1,3,7),bty='n')
# Use the lme function. This application
# assumes that each subject has a different
# identification value
weight.lme <- lme(Strength ~ Program*Timef,
random= ~ 1|Subj, data=weight,
method=c("REML"))
summary(weight.lme)
anova(weight.lme)
names(weight.lme)
# Use the gls( ) function to fit a
# model where the errors have a
# compound symmetry covariance structure
# within subjects. Random effects are
# not used to induce correlation.
weight.glscs <- gls(Strength ~ Program*Timef,
data=weight,
correlation = corCompSymm(form=~1|Subj),
method=c("REML"))
summary(weight.glscs)
anova(weight.glscs)
# Try an auto regressive covariance
# structures across time within
# subjects
weight.glsar <- gls(Strength ~ Program*Timef,
data=weight,
correlation = corAR1(form=~1|Subj),
method=c("REML"))
summary(weight.glsar)
anova(weight.glsar)
# Use an arbitray covariance matrix for
# observations at different time
# points within subjects
weight.gls <- gls(Strength ~ Program*Timef,
data=weight,
correlation = corSymm(form=~1|Subj),
weight = varIdent(form = ~ 1|Timef),
method=c("REML"))
summary(weight.gls)
anova(weight.gls)
# Compare the fit of various covariance
# structures.
anova(weight.gls, weight.glscs)
anova(weight.gls, weight.glsar)
# Treat time as a continuous variable and
# fit quadratic trends in strength
# over time
weight.time <- gls(Strength ~ Program+ Time+
Program*Time+Time^2+Program*Time^2,
data=weight,
correlation = corAR1(form=~1|Subj),
method=c("REML"))
summary(weight.time)
anova(weight.time)
# To compare the continuous time model to the
# model where we fit a different mean at each
# time point, we must compare likelihood values
# instead of REML likelihood values.
weight.glsarmle <- gls(Strength ~ Program*Timef,
data=weight,
correlation = corAR1(form=~1|Subj),
method=c("ML"))
weight.timemle <- gls(Strength ~ Program+ Time+
Program*Time+Time^2+Program*Time^2,
data=weight,
correlation = corAR1(form=~1|Subj),
method=c("ML"))
anova(weight.glsarmle, weight.timemle)
# Do not fit different quadratic trends
# for different programs
weight.timemle <- gls(Strength ~ Program + Time+
Program*Time+Time^2, data=weight,
correlation = corAR1(form=~1|Subj),
method=c("ML"))
anova(weight.glsarmle, weight.timemle)
# Fit a model with random regression coefficients
# for individual subjects
weight.timer <- lme(Strength ~ Program + Time+
Program*Time+Time^2,
random = ~ Time + Time^2 | Subj,
data=weight,
correlation = corAR1(form=~1|Subj),
method=c("REML"))
summary(weight.timer)
anova(weight.timer)
weight.timeru <- lme(Strength ~ Program + Time+
Program*Time+Time^2,
random = ~ Time + Time^2 | Subj,
data=weight,
method=c("REML"))
anova(weight.timer,weight.timeru)