# ---------------------------------------------------------------
 # Problem 15.4
 # ---------------------------------------------------------------


#  first read the data directly from Dr. Oehlerts website

Prb15.4 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/pr15.4",header=TRUE)

Prb15.4$judge   <- as.factor(Prb15.4$judge)
Prb15.4$A  <- as.factor(Prb15.4$A)
Prb15.4$B  <- as.factor(Prb15.4$B)
Prb15.4$C  <- as.factor(Prb15.4$C)
Prb15.4$D  <- as.factor(Prb15.4$D)

Prb15.4$logy <- log(Prb15.4$y) 


 # ------------- use Box-Cox transformation --------------------------------
library(MASS)
boxcox(y ~ judge + A*B*C*D, data=Prb15.4,lambda = seq(-2.00, 2.00, length = 50))
boxcox(y ~ judge + A*B*C*D, data=Prb15.4,lambda = seq(-0.50, 0.50, length = 50))
 # -------------------------------------------------------------------------

library(emmeans)

Prb15.4.lm1 <- lm(logy ~ judge + A*B*C*D, data = Prb15.4)

anova(Prb15.4.lm1)

lsmeans(Prb15.4.lm1, ~A)
lsmeans(Prb15.4.lm1, ~B)
lsmeans(Prb15.4.lm1, ~D)


par(mfrow=c(2,2))  
plot(Prb15.4.lm1)
par(mfrow=c(1,1)) 



 # ---------------------------------------------------------------
 # Problem 16.8
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Pr16.8 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/pr16.8",header=TRUE)

Pr16.8

Pr16.8$time2 <- 3*(Pr16.8$timing -1)     # convert 1, 2, 3  to 0, 3, 6

Pr16.8$sign  <- as.factor(Pr16.8$sign)
Pr16.8$timing  <- as.factor(Pr16.8$timing)
Pr16.8$interchange  <- as.factor(Pr16.8$interchange)


 # view sums of squares with model assuming fixed effects

Pr16.8.lm1  <- lm(speed ~ sign +sign:interchange +timing +sign:timing, data = Pr16.8)

anova(Pr16.8.lm1)


library(afex)

library(emmeans)

#-------------------------------------------------------------------------#
#  Mauchley test and downweighted degrees of freedom                      #
#-------------------------------------------------------------------------#

Pr16.8.afex1 <- aov_ez(id="interchange",dv="speed",data=Pr16.8,
                                         between=c("sign"),within=c("timing"))

 # the summary function produces Type III anova results assuming sphericity,
 # and includes the Mauchly tests for sphericity

summary(Pr16.8.afex1)

 # the nice function produces Type III anova results correcting for 
 # possible lack of sphericity by downweighting split plot degrees of freedom

nice(Pr16.8.afex1)




#-------------------------------------------------------------------------#
#  Random coefficient models                                              #
#-------------------------------------------------------------------------#

Pr16.8.afexrc1  <- mixed(speed ~ sign +time2 +sign:time2 +(1 + time2|interchange), 
                           data = Pr16.8)


summary(Pr16.8.afexrc1)
nice(Pr16.8.afexrc1)

qqnorm(summary(Pr16.8.afexrc1)$residuals)
qqline(summary(Pr16.8.afexrc1)$residuals)


#-------------------------------------------------------------------------#
#  Modeling covariance structure                                          #
#-------------------------------------------------------------------------#


library(nlme)

 # first we use the groupedData function to create a data structure approriate for
 # use with the repeated measures analysis

Pr16.8.longplot <- groupedData(speed ~ as.numeric(time2) | interchange, outer = ~ sign,
                                data = Pr16.8)

plot(Pr16.8.longplot,outer = ~ sign) 


 # now a structure for analysis

Pr16.8.long <- groupedData(speed ~ as.numeric(sign)*as.numeric(timing) | interchange, 
                                data = Pr16.8)


 # compound symmetric covariance structure

Pr16.8.cs <- gls(speed ~ sign * timing, data = Pr16.8.long, 
                          corr = corCompSymm(, form= ~ 1 | interchange) )
summary(Pr16.8.cs)
anova(Pr16.8.cs)


 # order 1 autoregressive covariance structure

Pr16.8.ar1 <- gls(speed ~ sign * timing, data = Pr16.8.long, 
                          corr = corAR1(, form= ~ 1 | interchange) )
summary(Pr16.8.ar1)
anova(Pr16.8.ar1)


 # order 1 heterogeneous autoregressive covariance structure

Pr16.8.arh1 <- gls(speed ~ sign * timing, data = Pr16.8.long, 
 corr = corAR1(form = ~ 1 | interchange),  weights = varIdent(form = ~ 1 | timing) )

summary(Pr16.8.arh1)
anova(Pr16.8.arh1)


 # an unstructured covariance matrix

Pr16.8.un <- gls(speed ~ sign * timing, data = Pr16.8.long, 
 corr=corSymm(form = ~ 1 | interchange),  weights = varIdent(form = ~ 1 | timing) )

summary(Pr16.8.un)
anova(Pr16.8.un)