# This example shows the analyses for the RBD         
# using the wheat data example we looked at in class  

# Entering the data and defining the variables:  

##########
##
# Reading the data into R:

my.datafile <- tempfile()
cat(file=my.datafile, "
 A       1       31.0
 A       2       39.5
 A       3       30.5
 A       4       35.5
 A       5       37.0
 B       1       28.0
 B       2       34.0
 B       3       24.5
 B       4       31.5
 B       5       31.5
 C       1       25.5
 C       2       31.0
 C       3       25.0
 C       4       33.0
 C       5       29.5
", sep=" ")

options(scipen=999) # suppressing scientific notation

wheat <- read.table(my.datafile, header=FALSE, col.names=c("variety", "block", "yield")) 
  
# Note we could also save the data columns into a file and use a command such as:
# wheat <- read.table(file = "z:/stat_516/filename.txt", header=FALSE, col.names = c("variety", "block", "yield"))

attach(wheat)

# The data frame called wheat is now created, 
# with three variables, variety, block, and yield.
##
#########


################################################################** 
# What if we ignored the blocks and just did a one-way CRD ANOVA?  
                  
# We specify that variety is a (qualitative) factor with the factor() function: 
                                    
# Making "variety" a factor:

variety <- factor(variety)

# The lm statement specifies that yield is the response                  
# and variety is the factor 
# The ANOVA table is produced by the anova() function         

wheat.fit <- lm(yield ~ variety);
anova(wheat.fit)

# The effect of variety is not significant at the 5% level (p-value = .0590). 
# This is NOT the best way to analyze these data!!!!                          

############################################################################### 
# Here we treat the experiment as a Randomized Block Design (RBD).              
# The lm() and anova() functions will do a standard ANOVA for a RBD. 
                                 
# We specify that block and variety are factors with the factor() function:

variety <- factor(variety)
block <- factor(block) 

wheat.RBD.fit <- lm(yield ~ variety + block);
anova(wheat.RBD.fit)

# Note that there are significant effects among treatments     
# (F statistic = 27.34, P-value = .0003).                      
# This implies that the mean yield is significantly different  
# for the different varieties of wheat.                        

# There is also significant variation among blocks             
# (F statistic = 20.68, P-value = .0003).                      
# This may or may not be of interest.                          
   

###############################################################################** 
# Question:  Is Variety A superior to the others?                                 
# We can again answer this question with a statement about a contrast.  

# Variety A vs. Others:

contrasts(variety) <- matrix(c(1, -1/2, -1/2), nrow=3, ncol=1)

print("Variety A vs. Others:")
summary(lm(yield ~ variety + block))$coef["variety1",]

# This code shows line labeled "variety1" which gives
# the two-sided P-value of t-test about contrast.

# Yes, clearly there is strong evidence (t = 7.28) that variety A has a superior  
# mean yield to the others.                                                       

# Two-sided P-value is given as 0.0000855.
# What is the P-value for this one-sided test?