##############################################
###  Lec3.r   STA 4211  Script file with:
###      I.  Code to imitate for Hw07 Question 1
###     II.  Code for lecture notes pp. 180 - 185
###
###   Q#1 requires use of the lm() function in R
###           and setting up dummy variables for factor effects model.
###          This process is illustrated here for the Castle Bakery dataset.
###           Needs dataset   CH19TA07.txt
###
##############################################
####  Castle Bakery data 3 x 2 balanced factorial design, w/ n=2
##############################################
###  Read the dataset into R.  Change the integer vectors of Factor A
###    and B levels into factor objects with appropriate labels of the levels.
###    Print the table of cell means.
castle <- read.table("Datasets/CH19TA07.txt",
                     col.names=c("Y","height","width","id"))
castle$height <- factor(castle$height,
                           labels=c("bottom","middle","top"))  #Factor A
castle$width <- factor(castle$width,
                       labels=c("regular","wide"))  #Factor B
cell.means <- tapply(castle$Y,list(castle$height,castle$width),
                     mean); cell.means
##  By hand, complete the table with row and column averages,
##    row and column main effects
##    and table of interaction effects.
##
##
## Next:  Find these effect estimates using "brute-force" linear
##    models approach.  
## Dummy variables for factor-effects formulation
##  
A1 <- rep(c(1,0,-1),c(4,4,4)); A1 # Height, Factor A, dummy var. 1
A2 <- rep(c(0,1,-1),c(4,4,4)); A2 # Height, Factor A, dummy var. 2
B1 <- rep(c(1,1,-1,-1),3); B1  # Width, Factor B, dummy var. 1 (and only)
A1B1 <- A1*B1  # Interaction, dummy 1, for cell (1,1)
A2B1 <- A2*B1  # interaction, dummy 2, for cell (2,1)
##  The lm() function for regression automatically includes the intercept
fit.lm <- lm(Y ~ A1 + A2 + B1 + A1B1 + A2B1, data=castle)
summary(fit.lm)
##
##
##  Get sums of squares and df for general linear F tests
##           first for main effects of A and B
##           then for the AB interaction
fit.full <- fit.lm
SSE.full <- sum(residuals(fit)^2); df.full <- fit$df.residual
SSE.full; df.full
Y <- castle$Y
fit.noA <- lm(Y ~ B1 + A1B1 + A2B1)
SSE.noA <- sum(residuals(fit.noA)^2); df.noA <- fit.noA$df.residual
SSE.noA; df.noA
F.A <- ((SSE.noA - SSE.full)/(df.noA - df.full))/(SSE.full/df.full); F.A
##
######################################
## The Randomized Complete Blocks Design, Chapter 21
######################################
## For a RCB design with three treatments and five blocks, use
##  R to assign treatments to units.
######################################
set.seed(16)
for (i in 1:5)
  {
    print(sample(1:3,replace=FALSE))
  }
##
##
##  Create the dataset for the Risk Premium example
##
confidence <- c(1,2,7,6,12,5,8,9,13,14,8,14,16,18,17)
method <- factor(rep(c("utility","worry","comp"),c(5,5,5)),
                 levels=c("utility","worry","comp"),
                 labels=c("utility","worry","comp"))
block <- factor(rep(1:5,3))
riskprem <- data.frame(confidence=confidence,method=method,block=block)
## Fit the additive model appropriate for randomized block design.
##  Obtain the ANOVA table; plot the data and the residuals.
riskprem.aov <- aov(confidence ~ block + method, data=riskprem)
summary(riskprem.aov)
anova(riskprem.aov)
attach(riskprem)
plot.design(confidence ~ block + method, data=riskprem)
interaction.plot(method, block, response=confidence)
plot(fitted(riskprem.aov),resid(riskprem.aov))
##
## Tukey's all pairwise comparisons
tapply(confidence,list(method),FUN=mean)
##  util 5.6   worry 9.8   comp 14.6
plot(TukeyHSD(riskprem.aov,which="method"))
##
## Relative efficiency of RCB versus CR design
nb <- 5; r <- 3; msbl <-42.833; msbl.tr <- 2.983
varhat.oneway <- ((nb - 1)*msbl + nb*(r-1)*msbl.tr)/(nb*r - 1)
varhat.oneway
ehat <- varhat.oneway/msbl.tr; ehat
##
## Don't forget to clean up, e.g.:
detach(riskprem)
######################################