# scab1.R
scabdata <- read.table(file="scab.dat",header=T) # read the data from an
# external file and grab the
# variable names from the
# 1st line of that file.
# The data are placed into
# a "data frame", which I've
# called scabdata
head(scabdata) #look at the first few rows of scabdata
#side-by-side boxplots. Notice the ordering of the treatments
plot(scabindx~trt,data=scabdata)
title(main="Boxplots of scab index for each treatment")
is.factor(scabdata$trt)
is.ordered(scabdata$trt)
scabdata$trtord <- ordered(scabdata$trt,
levels=c("0","F3","F6","F12","S3","S6","S12"))
#side-by-side boxplots (reordered treatments)
plot(scabindx~trtord,data=scabdata)
title(main="Boxplots of scab index for each treatment")
# Set "treatment" contrasts for ordered factors (for this application)
# instead of the default, which is orthogonal polynomial contrasts for
# ordered factors (we'll talk about what this means later)
op <- options(contrasts=c("contr.treatment", "contr.treatment"))
options()$contrasts # now you can see that they've been changed
op # and op holds the original value of contrasts so that we can reset back to
# the original more easily later.
m1 <- aov(scabindx~trtord,scabdata)
m1 # this "prints" the fitted model
summary(m1) # this summarizes the fitted model.
coefficients(m1) # this gives the actual fitted model parameters
summary.lm(m1) # this summarizes the fitted model treating it as generic lm
# not necessarily an aov model. This is an easy way to get some
# additional useful output.
# now estimate the mean for each treatment:
lsmeans(m1, specs = ~ trtord)
# Contrasts in the model paramters can be tested using the linearHypothesis
# function from the car package. This is useful, but sometimes a more
# convenient alternative is to define the corresponding contrast in the
# lsmeans and use the lsmeans() function to test the contrast.
( c1 <- c(6,-1,-1,-1,-1,-1,-1) )
( c2 <- c(0,1,1,1,-1,-1,-1) )
c3 <- c(0,2,-1,-1,2,-1,-1)
c4 <- c(0,0,-1,1,0,-1,1)
c5 <- c2*c3
c6 <- c2*c4
lsmeans(m1, specs = lsm ~ trtord,
contr=list(lsm=list( ctrl.vs.trted=c1/6, fall.vs.sp=c2/3 ) ) )
# The above contrasts can be tested with the linearHypothesis function.
# It is convenient to switch to a cell-means model first, though, to
# make the specification of the contrasts easier
m1a <- lm(scabindx~trtord-1,scabdata) # the "-1" removes the intercept
# now print out the coefficients from this fitted model, which are the
# estimated treatment means
m1a
linearHypothesis(m1a,hypothesis.matrix=c1) # ctrl vs treated
linearHypothesis(m1a,hypothesis.matrix=c2) # fall vs spring
linearHypothesis(m1a,hypothesis.matrix=rbind(c3,c4)) # dose
linearHypothesis(m1a,hypothesis.matrix=rbind(c5,c6)) # dose*season
options(op) # reset contrasts to their original values