ec <- read.table("http://users.stat.ufl.edu/~winner/data/egyptcttn.dat", header=F, col.names=c("variety","lum","lngrade","grade", "giza69","giza67","giza70","giza68")) attach(ec) head(ec); tail(ec) variety1 <- rep(1:5,each=4) par(mfrow=c(2,1)) plot(lum ~ grade, col=variety1) plot(lum ~ lngrade, col=variety1) par(mfrow=c(1,1)) ### Direct coding of "new" dummy bariables from "variety" g69 <- ifelse(variety == "Giza69" , 1, 0) g67 <- ifelse(variety == "Giza67" , 1, 0) g70 <- ifelse(variety == "Giza70" , 1, 0) g68 <- ifelse(variety == "Giza68" , 1, 0) g69gr <- lngrade*giza69 g67gr <- lngrade*giza67 g70gr <- lngrade*giza70 g68gr <- lngrade*giza68 cbind(g69gr,g67gr,g70gr,g68gr) (TSS <- sum((lum-mean(lum))^2)) ec.mod1 <- lm(lum ~ lngrade) summary(ec.mod1) anova(ec.mod1) confint(ec.mod1) coef(ec.mod1) X.grid <- seq(0,5,0.1) yhat1 <- coef(ec.mod1)[1] + coef(ec.mod1)[2]*X.grid plot(lum ~ lngrade, pch=16, xlim=c(0,5), ylim=c(82,92)) lines(X.grid,yhat1, col="blue", lwd=2) ec.mod2 <- lm(lum ~ lngrade + giza69 + giza67 + giza70 + giza68) summary(ec.mod2) anova(ec.mod2) confint(ec.mod2) coef(ec.mod2) X.grid <- seq(0,5,0.1) yhat1 <- coef(ec.mod2)[1] + coef(ec.mod2)[2]*X.grid yhat2 <- yhat1 + coef(ec.mod2)[3] yhat3 <- yhat1 + coef(ec.mod2)[4] yhat4 <- yhat1 + coef(ec.mod2)[5] yhat5 <- yhat1 + coef(ec.mod2)[6] plot(lum ~ lngrade, pch=16) lines(X.grid,yhat1, col="blue", lwd=2) lines(X.grid,yhat2, col="red", lwd=2) lines(X.grid,yhat3, col="green4", lwd=2) lines(X.grid,yhat4, col="purple", lwd=2) lines(X.grid,yhat5, col="brown", lwd=2) anova(ec.mod1, ec.mod2) ec.mod3 <- lm(lum ~ lngrade + giza69 + giza67 + giza70 + giza68 + g69gr + g67gr + g70gr + g68gr) summary(ec.mod3) anova(ec.mod3) confint(ec.mod3) X.grid <- seq(0,5,0.1) yhat1 <- coef(ec.mod3)[1] + coef(ec.mod3)[2]*X.grid yhat2 <- yhat1 + coef(ec.mod3)[3] + coef(ec.mod3)[7]*X.grid yhat3 <- yhat1 + coef(ec.mod3)[4] + coef(ec.mod3)[8]*X.grid yhat4 <- yhat1 + coef(ec.mod3)[5] + coef(ec.mod3)[9]*X.grid yhat5 <- yhat1 + coef(ec.mod3)[6] + coef(ec.mod3)[10]*X.grid plot(lum ~ lngrade, pch=16) lines(X.grid,yhat1, col="blue", lwd=2) lines(X.grid,yhat2, col="red", lwd=2) lines(X.grid,yhat3, col="green4", lwd=2) lines(X.grid,yhat4, col="purple", lwd=2) lines(X.grid,yhat5, col="brown", lwd=2) anova(ec.mod1,ec.mod2) anova(ec.mod2,ec.mod3) anova(ec.mod1,ec.mod3)