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