library(polynom) Ddata <- read.table(file="vitD.dat",header=T) # read the data in head(Ddata) is.factor(Ddata\$dose) is.ordered(Ddata\$dose) m1 <- aov(D~dose,Ddata) summary(m1) # this summarizes the fitted model. # R uses a different parameterization than SAS does in the version of the # model fit above in which there is an intercept. This means that # contrasts in the treatment means are not straightforward to specify # for the model above. It is easier to specify them for the cell means # version of the model, which can be specified just by removing the intercept. # Therefore, we refit model m1 using a cell means parameterization. Call this # model m1a m1a <- aov(D~dose-1,Ddata) summary(m1a) # this summarizes the fitted model. # Note that models m1 and m1a are equivalent models, but when the cell-means # parameterization of m1a is used, the F test that the summary function # produces is no longer the test of equal means. Rather it is a test that # all the treatment means are equal to 0. This is a hypothesis that is not # generally of any interest, so the F test given by the above summary(m1a) # code should be ignored. Nevertheless, model m1a is useful because its # parameter estimates are the estimated treatment means, and therefore # contrasts in the treatment means are contrasts in the parameters of model # m1a. Below orthogonal polynomial contrasts are specied and tested. # first treat dose as equally spaced ( hypmat1 <- matrix( c(-3,-1,1,3, 1,-1,-1,1, -1,3,-3,1), nrow=3,byrow=T ) ) linearHypothesis(m1a,hypothesis.matrix=hypmat1[1,]) # linear effect of dose linearHypothesis(m1a,hypothesis.matrix=hypmat1[2:3,])#nonlin effects of dose # Now treat dose as unequally spaced using the true dose levels doses <- c(2000,1000,400,0) # these are the true dose levels # The next two lines generate the orthogonal polynomial contrast coefficients # given the true spacing of dose junk <- poly.orth(doses,3,norm=T) (polyContrasts <- sapply(junk, predict, doses)) ( hypmat2 <- t(polyContrasts[,2:4]) ) # The t() function in the previous line transposes a matrix # (switches rows and columns) linearHypothesis(m1a,hypothesis.matrix=hypmat2[1,]) # linear effect of dose linearHypothesis(m1a,hypothesis.matrix=hypmat2[2:3,])#nonlin effects of dose