data <- scan(what = list(variety=0,block=0,dist=0,yield=0), multi.line=F)
 1 1 1 416.7  1 1 2 376.1  1 1 3 328.9  1 1 4 178.1 
 1 2 1 490.2  1 2 2 513.7  1 2 3 438.4  1 2 4 348.1 
 1 3 1 341.2  1 3 2 452.0  1 3 3 541.5  1 3 4 458.8 
 2 1 1 644.7  2 1 2 555.4  2 1 3 587.8  2 1 4 413.7 
 2 2 1 526.8  2 2 2 481.4  2 2 3 490.3  2 2 4 468.1  
 2 3 1 540.6  2 3 2 504.3  2 3 3 495.9  2 3 4 523.3 
 3 1 1 388.9  3 1 2 491.8  3 1 3 355.0  3 1 4 222.2 
 3 2 1 298.8  3 2 2 407.3  3 2 3 500.0  3 2 4 320.3 
 3 3 1 386.7  3 3 2 388.4  3 3 3 492.4  3 3 4 438.2 
 4 1 1 512.0  4 1 2 598.9  4 1 3 442.1  4 1 4 186.0 
 4 2 1 484.8  4 2 2 542.5  4 2 3 463.1  4 2 4 383.2 
 4 3 1 368.5  4 3 2 547.8  4 3 3 702.9  4 3 4 445.3 

data

prob155   <- data.frame(data)
rm(data)

prob155$variety         <- as.factor(prob155$variety)
prob155$block         <- as.factor(prob155$block)
prob155$dist   <- as.factor(prob155$dist)


 # check ANOVA table

anova(lm(yield ~ variety +block +variety:block +dist +variety:dist, data=prob155))


 # Using the afex package to get results like Proc Mixed

library(afex)

prob155.afex1  <- mixed(yield ~ variety +(1|block) +(1|variety:block) 
                                        +dist +variety:dist, data=prob155)

summary(prob155.afex1)

nice(prob155.afex1)

lsmeans(prob155.afex1, ~variety)

lsmeans(prob155.afex1, ~dist)

lsmip(prob155.afex1, variety ~ dist)

lsmip(prob155.afex1, dist ~ variety)


 # this code gives the epsilon-adjusted repeated measures results for Problem 15.5
 # using the car library

data <- scan(what = list(variety=0,block=0,yield1=0,yield2=0,yield3=0,yield4=0), multi.line=F)
 1 1  416.7  376.1  328.9  178.1 
 1 2  490.2  513.7  438.4  348.1 
 1 3  341.2  452.0  541.5  458.8 
 2 1  644.7  555.4  587.8  413.7 
 2 2  526.8  481.4  490.3  468.1 
 2 3  540.6  504.3  495.9  523.3 
 3 1  388.9  491.8  355.0  222.2 
 3 2  298.8  407.3  500.0  320.3 
 3 3  386.7  388.4  492.4  438.2 
 4 1  512.0  598.9  442.1  186.0 
 4 2  484.8  542.5  463.1  383.2 
 4 3  368.5  547.8  702.9  445.3 

data

prob15.5.multi <- data.frame(data)

prob15.5.multi  <- within(prob15.5.multi,
                    {avgyield <- (yield1+yield2+yield3+yield4)/4
                     variety  <- as.factor(variety)
                     block    <- as.factor(block)   } )


 # check the split plot sum of squares

anova(lm(avgyield ~ variety + block, data = prob15.5.multi))


 # use the car library for the multivariate approach to repeated measures,
 # to obtain the epsilon adjustment for general designs

library(car)

 # define the levels of the repeated measure

distance <- factor(c(1,2,3,4),levels=c("1","2","3","4"))

 # first fit a model to the whole plot structure

prob15.5.mmodel <- lm(cbind(yield1,yield2,yield3,yield4) ~ variety + block, prob15.5.multi)

 # create a data frame for the repeated measure

idata15.5  <- data.frame(distance)

 # now fit the repeated measures part of the model

prob15.5.rptms  <- Anova(prob15.5.mmodel, idata= idata15.5, idesign= ~distance, type=3)

 # gives the results with epsilon adjusted repeated measures analysis

summary(prob15.5.rptms, multivariate=FALSE)


 # all of the results match those of SAS Proc GLM except for the F statistic for distance

 # a very nice tutorial on this approach can be found at:
 # http://rtutorialseries.blogspot.com/2011/02/r-tutorial-series-one-way-repeated.html