# R Example with Mixed Effect Model  
# Training Data Example  
# These are NOT the same data as in the training example in the book!  

# Save the data file into a directory and 
# use the full path name:

training.data <- read.table(file = "z:/My Documents/teaching/stat_705/trainingdata.txt", 
header=FALSE, col.names = c('method', 'instructor', 'improvement'))

# attaching the data frame:

attach(training.data)

# Making "method" and "instructor" factors:

method <- factor(method)
instructor <- factor(instructor)
 
# Fitting the model and viewing the ANOVA table:

training.fit <- lm(improvement ~ method + instructor + method:instructor)
anova(training.fit)

# The lm() function has no equivalent to the RANDOM statement,
# but we must keep in mind that "instructor" has random levels.  

# For the test about the fixed-effect factor (method), we must tell R  
# that the correct denominator is MSAB (method*instructor). 
# The following code will do the correct test:

MSA <- anova(training.fit)["method", "Mean Sq"] 
MSAB <- anova(training.fit)["method:instructor", "Mean Sq"] 
df.A <- anova(training.fit)["method", "Df"]
df.AB <- anova(training.fit)["method:instructor", "Df"]

F.star <- MSA/MSAB; P.value <- pf(F.star, df.A, df.AB, lower.tail=F)

print(paste("F* =", round(F.star,3), "P.value =", round(P.value,4)))
  
# Note the correct results of this test (F*=135.7, P-value < .0001) are given   
# here.            
       
##########################

# Some code for the
# Tukey procedure:

# Specifying the family significance level:
alpha <- 0.05

df.A <- 1+anova(training.fit)["method", "Df"]
df.B <- 1+anova(training.fit)["instructor", "Df"]
n <- 1+(anova(training.fit)["Residuals", "Df"]/(df.A*df.B))

pairwise.diffs <- TukeyHSD(aov(training.fit),conf.level=0)$method[,"diff"]

my.std.error <- sqrt(2*MSAB/(df.B*n))

my.T <- qtukey(1-alpha, df.A, (df.A-1)*(df.B-1)) / sqrt(2)

# Tukey CIs for pairwise treatment mean differences among the methods:

lwr <- pairwise.diffs-my.T*my.std.error; upr <- pairwise.diffs+my.T*my.std.error

Tukey.CIs <- cbind(pairwise.diffs, lwr, upr)         
print(Tukey.CIs)         

# Since there was a significant effect due to method, we can use Tukey's procedure  
# to determine which methods specifically differ. At family significance level 0.05, 
# method 1 differs from each of methods 2, 3, and 4.  In addition,  methods 2 and 4 
# are significantly different.                                                                                                  

###############################################################################*