eiv <- read.table("eiv.dat",header=T)
eiv$Clin <- eiv$Clevel
eiv$Clevel <- ordered(eiv$Clevel)

library(pda)
mplot(eiv$Cactual,eiv$Resp,group=eiv$Clevel) # plot of data

eiv.aov <- aov(Resp ~ Clevel, eiv)	# analysis of variance
print(eiv.aov)
lsmean(eiv.aov)			# least squares means
plot(eiv.aov)				# plot of ANOVA fit
coef(eiv.aov)/sqrt(2)			# compare linear coef with SAS

eiv.reg <- lm(Resp ~ Clin, eiv)	# regression on discrete levels
print(eiv.reg)
summary.aov(eiv.reg)
plot(eiv.reg)				# default plot of reg fit

					# test for linearity
eiv.lin <- lm(Resp ~ Clin+Clevel, eiv, singular.ok=T)
print(eiv.lin)
summary.aov(eiv.lin)

eiv.act <- lm(Resp ~ Cactual, eiv)	# regression on actual carbon
print(eiv.act)
summary.aov(eiv.act)
plot(eiv.act)				# default plot of fit

					# error in variables
eiv.ervar <- lm(Resp ~ Clevel+Cactual, eiv, singular.ok=T)
print(eiv.ervar)
summary.aov(eiv.ervar)
plot(eiv.ervar)			# plot of fit

					# plot of data and all fits
mplot(eiv$Cactual,eiv$Resp,group=eiv$Clevel)
mlines(eiv$Cactual,predict(eiv.aov),group=eiv$Clevel,lty=1)
mlines(eiv$Cactual,predict(eiv.reg),group=eiv$Clevel,lty=2)
mlines(eiv$Cactual,predict(eiv.act),group=eiv$Clevel,lty=3)
mlines(eiv$Cactual,predict(eiv.ervar),group=eiv$Clevel,lty=4)
legend(3.5,2250,legend=c("aov","reg","act","erv"),lty=1:4)

				# THIS LAST RUN IS FOR YOUR INFORMATION ONLY
				# error in variables depending on group
eiv.ervar2 <- lm(Resp ~ Clevel+Cactual+Cactual*Clevel, eiv, singular.ok=T)
summary.aov(eiv.ervar2)