# read data (level is score as a factor)
damage <- read.table("damage.dat",header=T)
damage$level <- factor(damage$score)
damage$logit <- log(damage$loss/(1-damage$loss))

# loss score against yield
pp.loss <- lm(yield~loss,damage)
attach(damage)
mplot(loss,yield,group=as.character(score),
   xlim=c(0,1),xaxt="n",xlab="(a) loss (%)",ylab="")
axis(1,at=c(0,0.5,1))
tmp <- order(loss)
lines(sort(loss),predict(pp.loss)[order(loss)],lty=2)
se.bar(0.6,93,std.dev(pp.loss),cap="SD")
detach("damage")

pp.level <- lm(yield~level,damage)
pp.score <- lm(yield~score,damage)
pp.mean <- tapply(damage$yield,damage$score,mean)
attach(damage)
plot(score,yield,
   xlim=c(1,5),xaxt="n",xlab="(b) score",ylab="yield")
axis(1,at=1:5)
lines(sort(score),predict(pp.score)[order(score)],lty=2)
for (i in 1:4)
	lines(i+c(-.4,.4),rep(pp.mean[i],2))
se.bar(3.5,93,std.dev(pp.score),cap="SD")
detach("damage")

# test for linearity
pp.lin <- aov(yield~score+level,damage)
summary(pp.lin)

# score error in variables (\myexpage{plerr}; \myfigpage{plerr})
pp.both <- aov(yield~level+logit,damage)
summary(pp.both)
attach(damage)
mplot(logit,yield,group=as.character(score),
	xlab="(b) logit(loss)",ylab="")
tmp <- order(logit)
mlines(logit,predict(pp.both),group=score,lty=1)
se.bar(0.6,93,std.dev(pp.both),cap="SD")
detach("damage")

pp.logit <- aov(yield~logit,damage)
summary(pp.logit)
attach(damage)
mplot(logit,yield,group=as.character(score),
	xlab="(b) logit(loss)",ylab="yield")
lines(sort(logit),predict(pp.logit)[order(logit)],lty=2)
se.bar(0.6,93,std.dev(pp.logit),cap="SD")
detach("damage")