# R example for analyzing a BIBD  

# These are the taste test data from class and from chapter 28.  

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

tastetest.data <- read.table(file = "z:/My Documents/teaching/stat_705/tastetestdata.txt", 
header=FALSE, col.names = c('block', 'formulation', 'rating'))

# attaching the data frame:

attach(tastetest.data)

# Making "block" and "formulation" factors:

block <- factor(block)
formulation <- factor(formulation)

# We use R to analyze the data from the BIBD.  

# If blocks are considered random, we use a mixed model:
# formulations are fixed, blocks are random. 

# To get a mixed effects model fit, we can install the "lme4" package in R:

# To install the "lme4" package:
# Go to "Packages" menu and choose
# "Install package(s) from CRAN"
# scroll to the "lme4" package and choose it
# An Internet connection is needed to install it.

# Then you can proceed with the code...

# Load the "lme4" package:

library(lme4)   

# A mixed model can be fit with the lmer() function in the lme4 package:
# This syntax tells R that "formulations" has fixed levels 
# but "block" has random levels:

tastetest.fit <- lmer(rating ~ formulation + (1|block) )

summary(tastetest.fit)

# Note the MSE is 0.80228 (Find this in the summary output.):

MSE <- 0.80228

anova(tastetest.fit)

# We can do an F-test for significant "formulation" effects
# by dividing the MS(Formulation) by the MSE:

# The (approximate) error d.f. is not shown, but can be found by dividing
# SSE by MSE:

error.df <- sum((resid(tastetest.fit))^2)/MSE

F.star <- anova(tastetest.fit)["formulation","Mean Sq"] / MSE; 
print(paste("F*=", round(F.star,4)))
P.value <- pf(F.star, df1 = anova(tastetest.fit)["formulation","Df"], df2=error.df, lower.tail=F)
print(paste("P-value=", round(P.value,4)))

# The significant effect due to formulation is seen by the     
# F* = 16, P-value < .0001 in the output.  


############################################################

# If blocks are considered fixed, we can simply use the lm() and
# anova() functions to fit the model:

# When it is of interest to test for TREATMENT effects:
# The BLOCKING FACTOR should go FIRST in the model!
# This makes a difference when the design is INCOMPLETE.

tastetest.fit.fixed <- lm(rating ~ block + formulation)
anova(tastetest.fit.fixed)     

# The significant effect due to formulation is seen by the    
# F* = 15.66, P-value < .0001 in the output.  

# If we want to test for BLOCK effects:
# The TREATMENT factor should go first in the model:

tastetest.fit.fixed2 <- lm(rating ~ formulation + block)
anova(tastetest.fit.fixed2)                              

# There is also a significant block effect (F* = 5.56, P-value = .0015)  


# Some code for the
# Tukey procedure:

my.std.error <- coef(summary(tastetest.fit.fixed))["formulation2", "Std. Error"]

# Since there are 5 levels of the treatment factor:
Tukey.cutoff <- qtukey(0.95, 5, anova(tastetest.fit.fixed)["Residuals", "Df"])*my.std.error / sqrt(2)

# Viewing the fitted model coefficients:
coef(tastetest.fit.fixed)
# Picking out the ones associated with "formulations":
mycoefs <- c(0,coef(tastetest.fit.fixed)[11:14])

# Showing which differences are significant at the 0.05 family level:
abs(outer(mycoefs, mycoefs, "-")) > Tukey.cutoff

# From the Tukey simultaneous comparisons, we see:                 
# Treatments 1 and 5 are not significantly different,              
# Treatments 2 and 3 are not significantly different,              
# Treatments 2 and 4 are not significantly different,              
# Treatments 3 and 4 are not significantly different,              
# All other pairs of treatment means are significantly different.  
# This uses a 0.05 family significance level.