####################################################### ### ANCOVA example (glue data from the ppt lectures) ####################################################### glue <- read.table("Data/glue.txt",header=TRUE) # convert formulation codes to factor variables for proper statistical handling glue\$Formulation <- as.factor(glue\$Formulation) # fit linear models: full, thickness only, formulation only full.lm <- lm(Strength ~ Formulation + Thickness, data=glue) thick.lm <- lm(Strength ~ Thickness, data=glue) formu.lm <- lm(Strength ~ Formulation, data=glue) # test for sig of formulation anova(thick.lm,full.lm) # test for sig of thickness anova(formu.lm,full.lm) # plot the data & add the fitted lines for each formulation plot(glue\$Thickness,glue\$Strength,type="n",xlab="Thickness",ylab="Strength") text(glue\$Thickness,glue\$Strength,as.character(glue\$Formulation)) abline(58.92788,-0.95445,lty=1) # form 1 abline(58.92788+0.63466,-0.95445,lty=2) # form 2 abline(58.92788+0.87644,-0.95445,lty=3) # form 3 abline(58.92788+0.00911,-0.95445,lty=4) # form 4 legend(14.5,51,c("Formulation 1","Formulation 2","Formulation 3","Formulation 4"),lty=c(1,2,3,4)) # note however that due to no differences in formulations, it suffices to plot # the following line, obtained from thick.lm summary(thick.lm) abline(59.9294,-1.0044,col="blue") ####################################################### ### ANCOVA example (golf course speed data from ppt lectures) ####################################################### golf <- read.table("Data/Ott5thEdDataCh16/ch16-case_study.csv",header=TRUE,sep=",") region <- factor(golf\$R) cult <- factor(golf\$C) humid <- golf\$HUM speed <- golf\$SPEED # plot speed vs. humid #pdf(file="Plots/golf1.pdf",pointsize=3,width=6,height=5) plot(humid,speed,type="n",xlab="Humidity",ylab="Speed",main="Speed vs. Humidity by Cultivar (1=C1, 2=C2, 3=C3)") text(humid,speed,as.character(cult)) #dev.off() # Fit Models 1,2,3 model1 <- lm(speed ~ humid + region + cult + humid*cult) model2 <- lm(speed ~ humid + region + cult) model3 <- lm(speed ~ humid + region) # compare models 1 to 2 (test significance of different cult slopes) anova(model2,model1) # compare models 2 to 3 (test significance of cult) anova(model3,model2) # final model anova(model2) summary(model2) # plot speed vs. humid (add fitted cult lines) #pdf(file="Plots/golf2.pdf",pointsize=3,width=6,height=5) plot(humid,speed,type="n",xlab="Humidity",ylab="Speed",main="Speed vs. Humidity With Fitted Cultivar Lines") text(humid,speed,as.character(cult)) abline(8.421762,-0.022765,lty=1) # cult 1 abline(8.421762+0.917971,-0.022765,lty=2) # cult 2 abline(8.421762+1.885567,-0.022765,lty=3) # cult 3 legend(80,9.8,c("Cultivar 1","Cultivar 2","Cultivar 3"),lty=c(1,2,3)) #dev.off() # average region (block) effect mregion <- sum(c(-0.072989,-0.084832,-0.186642,0.340397,0.434006,0.433041,0.252458))/8 # adjusted (for region and humid) means of C1, C2, C3 muC1 <- 8.421762-0.022765*mean(humid)+mregion muC2 <- 8.421762+0.917971-0.022765*mean(humid)+mregion muC3 <- 8.421762+1.885567-0.022765*mean(humid)+mregion # print these means c(muC1,muC2,muC3) # Check model fit #pdf(file="Plots/golf3.pdf",pointsize=3,width=6,height=5) par(mfrow=c(2,1)) qqnorm(model2\$res) plot(model2\$fitted,model2\$res); abline(0,0) par(mfrow=c(1,1)) #dev.off() # all pairwise Tukey comparisons among cult means (needs "multcomp" package) library(multcomp) summary(glht(model2, linfct=mcp(cult="Tukey")))