# ---------------------------------------------------------- # Exercise 15.1 # ---------------------------------------------------------- data <- scan(what = list(diet=0,subject=0,time=0,glucose=0), multi.line=F) 1 1 15 28 1 1 30 34 1 1 45 32 1 2 15 15 1 2 30 29 1 2 45 27 1 3 15 12 1 3 30 33 1 3 45 28 1 4 15 21 1 4 30 44 1 4 45 39 2 5 15 22 2 5 30 18 2 5 45 12 2 6 15 23 2 6 30 22 2 6 45 10 2 7 15 18 2 7 30 16 2 7 45 9 2 8 15 25 2 8 30 24 2 8 45 15 3 9 15 31 3 9 30 30 3 9 45 39 3 10 15 28 3 10 30 27 3 10 45 36 3 11 15 24 3 11 30 26 3 11 45 36 3 12 15 21 3 12 30 26 3 12 45 32 data prob15.1 <- data.frame(data) # using the within function to convert variables to factor prob15.1 <- within(prob15.1, { diet <- factor(diet) subject <- factor(subject) time <- factor(time) } ) # get an interaction plot (profile plot) with(prob15.1,interaction.plot(diet,time,glucose)) # first check the ANOVA table anova(lm(glucose ~ diet +subject:diet +time +diet:time, data = prob15.1)) # ----------------------------------------------------------------------- # Split Plot Analysis - gives mostly the same information as Proc Mixed # ----------------------------------------------------------------------- library(lme4) library(lmerTest) prob15.1.lme1 <- lmer(glucose ~ diet +(1|subject:diet) +time +diet:time, data = prob15.1) anova(prob15.1.lme1) lsmeans(prob15.1.lme1) difflsmeans(prob15.1.lme1) # residual plots par(mfrow=c(2,1)) plot(fitted(prob15.1.lme1),resid(prob15.1.lme1),xlab="Predicted Values", ylab="Residuals",main="Residual by Predicted plot") qqnorm(resid(prob15.1.lme1)) qqline(resid(prob15.1.lme1)) par(mfrow=c(1,1)) # or with afex library(afex) prob15.1.afex1 <- mixed(glucose ~ diet +(1|subject:diet) +time +diet:time, data = prob15.1) # ----------------------------------------------------------------------- # Mauchly tests and adjusted p values # ----------------------------------------------------------------------- # using the aov_ez function from the afex package prob15.1.afex2 <- aov_ez(id="subject",dv="glucose",data=prob15.1,between=c("diet"),within=c("time")) summary(prob15.1.afex2) nice(prob15.1.afex2) # ---------------------------------------------------------- # Exercise 15.5 # ---------------------------------------------------------- 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) # an easier way to get the sphericity test and epsilon-adjustments in the afex package prob155$vb <- 10*as.numeric(prob155$variety) +as.numeric(prob155$block) prob15.5.afexRPTMS2 <- aov_car(yield ~ variety +block + Error(vb/(dist)),data=prob155) summary(prob15.5.afexRPTMS2) # ----------------------------------------------------------------------------------- # this code gives the epsilon-adjusted repeated measures results for Problem 15.5 # using the car library in a more complicated way 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 # ---------------------------------------------------------- # Exercise 11.7 # ---------------------------------------------------------- data <- scan(what = list(rep=0,day=0,a=0,b=0,c=0,d=0,y=0), multi.line=F) 1 1 0 0 0 0 40 1 1 1 1 0 0 33 1 1 1 0 1 0 31 1 1 0 1 1 0 38 1 1 1 0 0 1 22 1 1 0 1 0 1 37 1 1 0 0 1 1 49 1 1 1 1 1 1 30 1 2 1 0 0 0 24 1 2 0 1 0 0 31 1 2 0 0 1 0 27 1 2 1 1 1 0 23 1 2 0 0 0 1 48 1 2 1 1 0 1 35 1 2 1 0 1 1 29 1 2 0 1 1 1 37 2 3 0 0 0 0 43 2 3 1 1 0 0 30 2 3 1 0 1 0 30 2 3 0 1 1 0 32 2 3 1 0 0 1 26 2 3 0 1 0 1 33 2 3 0 0 1 1 40 2 3 1 1 1 1 31 2 4 1 0 0 0 28 2 4 0 1 0 0 35 2 4 0 0 1 0 28 2 4 1 1 1 0 20 2 4 0 0 0 1 44 2 4 1 1 0 1 36 2 4 1 0 1 1 25 2 4 0 1 1 1 34 data prob11.7 <- data.frame(data) # using the within function to convert variables to factor prob11.7 <- within(prob11.7, { rep <- factor(rep) day <- factor(day) a <- factor(a) b <- factor(b) c <- factor(c) d <- factor(d) } ) # see the ANOVA table anova(lm(y ~ rep + day:rep + a*b*c*d, data = prob11.7 )) # it is easier to see the results using the afex package library(afex) prob11.7.afex1 <- mixed(y ~ rep + (1|day:rep) +(a+b+c+d)^3, data = prob11.7) summary(prob11.7.afex1) nice(prob11.7.afex1) lsmeans(prob11.7.afex1, ~ a) lsmeans(prob11.7.afex1, ~ b) lsmip(prob11.7.afex1, a ~ b|c) lsmip(prob11.7.afex1, a ~ b|d)