# R source file for the solutions to problems given in Chapter 4
# Time-stamp: Time-stamp: <2008-02-28 00:15:03 chl>

##############################################################################
WCGS <- read.table("wcgs.txt",header=TRUE,sep="\t")
WCGS$arcus_f <- as.factor(WCGS$arcus)
WCGS.lm1 <- lm(chol~arcus_f,data=WCGS)
summary(WCGS.lm1)


##############################################################################
levels(WCGS$arcus_f)
WCGS$arcus_f2 <- relevel(WCGS$arcus_f,"1")
levels(WCGS$arcus_f2)
WCGS.lm2 <- lm(chol~arcus_f2,data=WCGS)
summary(WCGS.lm2)
with(WCGS, tapply(chol,arcus_f,mean,na.rm=T))
with(WCGS, tapply(chol,arcus_f2,mean,na.rm=T))
table(WCGS$arcus_f2)


##############################################################################
WCGS$arcus2 <- WCGS$arcus
WCGS$arcus2[WCGS$arcus2==0] <- 2
WCGS.lm3 <- lm(chol~arcus2,data=WCGS)
summary(WCGS.lm3)


##############################################################################
WCGS$arcus3 <- WCGS$arcus
WCGS$arcus3[WCGS$arcus==0] <- -1
WCGS.lm5 <- lm(chol~arcus3,data=WCGS)
summary(WCGS.lm5)


##############################################################################
yhat <- fitted(WCGS.lm3)
length(yhat)
length(WCGS$chol[!is.na(WCGS$arcus)])


##############################################################################
print(idx <- attr(WCGS.lm3$model,"na.action"))


##############################################################################
cor(yhat,WCGS$chol[-as.numeric(idx)])^2
summary(WCGS.lm3)$r.squared


##############################################################################
HERS <- read.table("hersdata.txt",header=TRUE,na.strings=".")
HERS$physactf <- as.factor(HERS$physact)
physact.desc <- c("1=Much less active","2=Somewhat less active","3=About as active","4=Somewhat more active","5=Much more active")
HERS.lm <- lm(glucose~physactf,data=HERS,subset=diabetes==0)
summary(HERS.lm)
confint(HERS.lm)
par(mfrow=c(1,2))
plot(glucose~physactf,data=HERS)
with(HERS[HERS$diabetes==0,], stripchart(glucose~physactf,vertical=T,cex=.8,pch=1))
abline(HERS.lm1 <- lm(glucose~physact,data=HERS,subset=diabetes==0),col="red")


##############################################################################
anova(HERS.lm)


##############################################################################
head(model.matrix(HERS.lm))


##############################################################################
HERS.lm2 <- lm(glucose~-1+physactf,data=HERS,subset=diabetes==0)
summary(HERS.lm2)


##############################################################################
head(model.matrix(HERS.lm2))


##############################################################################
ni <- 20
x <- gl(5,ni)
y <- c(rnorm(ni,10,2),replicate(3,rnorm(ni,12,2)),rnorm(ni,6,2.5))
xy <- data.frame(x=x,y=y)
summary(xy)
boxplot(y~x,xlab="Group",ylab="Response")
lines(1:5,tapply(y,x,mean),type="b",col="red",lwd=2)
abline(h=mean(y[x==1]),col="lightgray",lwd=2)


##############################################################################
xy.lm <- lm(y~x,xy)
library(contrast)
contrast(xy.lm,list(x='1'),list(x='2'))


##############################################################################
C <- c(1,1,-1,-1)
aov(cbind(x11,x12,x21,x22) %*% C ~ group, data)


##############################################################################
options(contrasts = c("contr.helmert", "contr.poly"))


##############################################################################
help.search("contrast")


##############################################################################
summary(xy.lm)


##############################################################################
n <- 100
x <- gl(2,n/2,n,labels=0:1)
y <- c(rnorm(n/2,25,3),rnorm(n/2,30,2.5))
tapply(y,x,mean)
tapply(y,x,var)
plot(x,y)


##############################################################################
y.lm <- lm(y~x)
summary(y.lm)


##############################################################################
t.test(y~x,var.equal=TRUE)


##############################################################################
x <- gl(3,n/3,n,labels=0:2)
y <- as.numeric(x) + rnorm(n)
by(y,x,summary)
plot(x,y)
stripchart(y~x,vertical=T)
points(1:3,tapply(y,x,mean),pch="x",cex=2,col="red")
abline(lm(y~as.numeric(x)),col="red")


##############################################################################
summary(lm(y~x))


##############################################################################
mean(y[x==1])-mean(y[x==0])


##############################################################################
summary(aov(y~x))


##############################################################################
summary.lm(aov(y~x))


##############################################################################
model.tables(aov(y~x))