# R commands for Chapter 8

 

# Section 8.1

 

# Read the fecal fat data into a data frame (needs only to be done once if the workspace is saved)

fecfat=read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/v11/fecfat.txt",header=T)

 

 

# One-way ANOVA for the fecal fat example (cf Table 8.2, page 255):

fit.fecfat1=lm(fecfat~factor(pilltype), data=fecfat)

anova(fit.fecfat1)

 

# Two-way ANOVA for the fecal fat example (cf Table 8.3, page 256):

fit.fecfat2=lm(fecfat~factor(subject)+factor(pilltype), data=fecfat)

anova(fit.fecfat2)

 

 

# Section 8.3

 

# Read the data on bith weight and birth order into a data frame (needs only to be done once if the workspace is saved)

babies=read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/v11/gababies.txt",header=T)

 

# Summary statistic for first and last born babies (cf Table 8.4, page 263):

summary(babies[ ,c(10,8,5)])

 

# Regression of difference in birthweight on centered initial age (cf Table 8.5, page 263):

fit.diffweight=lm(delwght~cinitage,data=babies,subset=(birthord==1))

summary(fit.diffweight)

 

# Repeated measures regression of birthweight on birth order (only first and fifth) and centered initial age (cf Table 8.6, page 264):

#These commands use the package "gee", which first need to be installed (cf. the introduction to R given at the first lecture)

library(gee)??????? # load the gee-package

fit.repweight=gee(bweight~I(birthord==5)+cinitage+I(birthord==5):cinitage,

???????????????????????????????? id=momid,data=babies,subset=((birthord==1)|(birthord==5)),corstr="exchangeable")

summary(fit.repweight)

 

# Regression of final birthweight on centered initial age, adjusting for first birthweight (cf Table 8.7, page 265):

fit.lastweight=lm(lastwght~cinitage+initwght,data=babies,subset=(birthord==5))

summary(fit.lastweight)

 

 

# Section 8.4

 

# Plot of birthweight versus birth order (cf Fig 8.1, page 267):

plot(babies$birthord,babies$bweight,ylim=c(0,5000),xlab="Birth Order",ylab="Birthweight (g)")

lines(lowess(babies$birthord,babies$bweight),lwd=2,col='red')

 

# Boxplot of birthweight versus birth order (cf Fig 8.2, page 268):

boxplot(bweight~birthord,ylim=c(0,5000), data=babies,xlab="Birth Order",ylab="Birthweight (g)")

 

# Correlation of birthweights for different? birth orders (cf Table 8.8, page 268):

birthweight=matrix(0,nrow=200,ncol=5)

for (i in 1:5) birthweight[,i]=babies$bweight[babies$birthord==i]

birthweight=data.frame(birthweight)

names(birthweight)=c("first","second","third","fourth","fifth")

cor(birthweight)

 

# Correlation of birthweights for different? birth orders (cf Fig 8.3, page 269):

plot(birthweight)

 

# Generalized estimating equations model with and without robust standard errors (cf Tables 8.9 and 8.10, page 272):

fit.gee=gee(bweight~birthord+initage, id=momid, data=babies, corstr="exchangeable")

summary(fit.gee)

 

 

# Section 8.5

 

# Random effects linear regression model for birthweight (cf Table 8.14, page 277):

#These commands use the package "nlme", which first need to be installed (cf. the introduction to R given at the first lecture)

library(nlme)??????? # load the nlme-package

fit.random=lme(bweight~birthord+initage, random=~1|momid, data=babies, method="ML")

summary(fit.random)