# read data
leprosy <- read.table("leprosy.dat",header=T)

# before and after scores
lep.fit <- aov(after~before+drug, leprosy)

attach(leprosy)
mplot(before,after,pch=drug)
mlines(before, predict(lep.fit), group=drug, lty=1:3)

tmp <- data.frame(drug=levels(drug),before=rep(mean(before),3))
abline(v=tmp$before[1],lty=4)
points(tmp$before,predict(sc.fit,tmp),pch=18)
points(tapply(before,drug,mean),tapply(after,drug,mean),pch=5)

legend(5,20,c("A","B","C"),lty=1:3)
lsd.bar(lep.fit,sc,xpos=18,ypos=5)
detach()

# plot residuals
mplot(predict(lep.fit),resid(lep.fit),pch=leprosy$drug,
	xlab="predicted",ylab="residual")
abline(0,0,lty=2)

legend(0,8,c("A","B","C"),lty=1:3)
se.bar(18,-3,std.dev(lep.fit),cap="SD")

# diagnostic partial plots
lep.drug <- aov(after~drug,leprosy)

attach(leprosy)
mplot(before,resid(lep.drug),pch=drug,
	ylim=c(-12,12),xlab="(a) before",ylab="")
abline(0,0,lty=2)
tmp <- data.frame(after=resid(lep.drug),before=before)
lines(sort(before), predict(lm(after~before,tmp))[order(before)],lty=1)
se.bar(17,-5,std.dev(lep.drug),cap="SD")
detach()

lep.before <- aov(after~before,leprosy)

attach(leprosy)
mplot(predict(lep.before),resid(lep.before),pch=drug,
	ylim=c(-12,12),xlab="(b) predicted",ylab="residual")
abline(0,0,lty=2)
tmp <- data.frame(y=resid(lep.before),drug=drug)
mlines(predict(lep.before),predict(aov(y~drug,tmp)),group=drug)
legend(0,-6,c("A","B","C"),lty=1:3)
se.bar(14,-9,std.dev(lep.before),cap="SD")
detach()

# sequential tests
lep.fit <- aov(after~drug+before,leprosy)
summary(lep.fit)
summary(aov(after~before+drug,leprosy))

# adjusted tests
drop1(lep.fit,formula(lep.fit)

# different slopes by factor
attach(leprosy)
mplot(predict(lep.fit),resid(lep.fit),pch=drug,
	xlab="(a) predicted",ylab="residual")
abline(0,0,lty=2)
tmp <- data.frame(after=resid(lep.fit),before=predict(lep.fit),drug=drug)
mlines(tmp$before,predict(aov(after~before*drug,tmp)),group=drug)

legend(13,9,c("A","B","C"),lty=1:3)
se.bar(17,-3,std.dev(sc.fit),cap="SD")
detach()

lep.ifit <- aov(after~before*drug,leprosy)
summary(lep.ifit)

attach(leprosy)
mplot(before,after,pch=drug,
	xlab="(b) before",ylab="after")
sc.ifit <- aov(y~x*drug,sc)
mlines(before,predict(sc.ifit),group=drug)

tmp <- data.frame(drug=levels(drug), before=rep(mean(before),3))
abline(v=tmp$before[1],lty=4)
points(tmp$before, predict(sc.ifit,tmp), pch=18)
points(tapply(before,drug,mean), tapply(after,drug,mean), pch=5)

se.bar(18,5,std.dev(lep.ifit),cap="SD")
detach()