library(MASS)
help(anorexia)

library(lattice)
attach(anorexia)

contrasts(Treat) <- contr.sum(length(levels(Treat)))
levels(Treat)
xyplot(Postwt~Prewt|Treat,data=anorexia,
	panel=function(x,y){panel.lmline(x,y,type="l")
	panel.abline(b=1,a=0,lty=3)
	panel.xyplot(x,y,col=2,pch=16)})


Prewt.centered <- Prewt - mean(Prewt)
xyplot(Postwt~Prewt.centered|Treat,data=anorexia,
	panel=function(x,y){panel.lmline(x,y,type="l")
	panel.abline(b=1,a=0,lty=3)
	panel.xyplot(x,y,col=2,pch=16)})

contrasts(Treat) <- contr.sum(length(levels(Treat)))
lm.obj <- lm(Postwt~Treat+Prewt.centered)
anova(lm.obj)
summary(lm.obj)
levels(Treat)

beta <- lm.obj$coef
beta[1]+beta[2]
beta[1]+beta[3]
beta[1]-beta[2]-beta[3]


xyplot(lm.obj$residuals~lm.obj$fitted.values|Treat,data=anorexia,xlab="Fitted",ylab="Residuals",
	panel=function(x,y){panel.lmline(x,y,type="l")
	panel.xyplot(x,y,col=1,pch=16)})

############################### interaction model #################
contrasts(Treat) <- contr.sum(length(levels(Treat)))

lm.obj <- lm(Postwt~Treat*Prewt.centered)
anova(lm.obj)
summary(lm.obj)
levels(Treat)


xyplot(Postwt~Prewt|Treat,data=anorexia,
	panel=function(x,y){panel.lmline(x,y,type="l")
	panel.abline(b=1,a=0,lty=3)
	panel.xyplot(x,y,col=2,pch=16)})


xyplot(lm.obj$residuals~lm.obj$fitted.values|Treat,data=anorexia,xlab="Fitted",ylab="Residuals",
	panel=function(x,y){panel.lmline(x,y,type="l")
	panel.abline(b=0,a=0,lty=3)
	panel.xyplot(x,y,col=1,pch=16)})

lm.obj <- lm(Postwt~Treat*Prewt)
anova(lm.obj)
summary(lm.obj)
model.matrix(lm.obj)

I <- Treat=="FT"
y <- Postwt[I]
x <- Prewt[I]

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

data.new <- data.frame(x=seq(from=70,to=100,by=1))
pred <- predict(lm.obj,newdata=data.new,interval = "confidence",level = 0.95,type="response")
pred

plot(x=c(70,100),y=c(67,120),xlab="Pre weight",ylab="Post weight",type="n")
points(x,y,pch=16)

lines(seq(from=70,to=100,by=1),pred[,1],col=1)
lines(seq(from=70,to=100,by=1),pred[,2],col=2)
lines(seq(from=70,to=100,by=1),pred[,3],col=2)


pred <- predict(lm.obj,newdata=data.new,interval = "prediction",level = 0.95,type="response")
pred
lines(seq(from=70,to=100,by=1),pred[,2],col=3)
lines(seq(from=70,to=100,by=1),pred[,3],col=3)