# R program for RCBD with Hardness Testing Experiment

      
 tip1 <- c(9.3,  9.4,  9.6, 10.0)
 tip2 <- c(9.4,  9.3,  9.8,  9.9)
 tip3 <- c(9.2,  9.4,  9.5,  9.7)
 tip4 <- c(9.7,  9.6, 10.0, 10.2)
 
 data <- data.frame(rbind(tip1, tip2, tip3, tip4))
 
 data

 apply(data, 1, mean) # Treatment averages
                     # "1" refers to the first of the pair (row,column)


 apply(data,2,mean) 

# Arrange responses by row
 y <- c(tip1, tip2, tip3, tip4)
 y

# Generate block index for response values
 blocks <- as.factor(rep(1:4, times = 4))
 blocks

# Generate treatment level for response values
 treatments <- as.factor(rep(1:4, each = 4))
 treatments

 par(mfrow=c(1,2))  # A 1 by 2 panel of plots
 boxplot(y~blocks, xlab="blocks")
 boxplot(y~treatments, xlab="treatments")


# Get the ANOVA table of RCBD
 g <- lm(y~treatments + blocks)
 anova(g)

# to check normality assumption and view residuals
 windows()
 par(mfrow=c(2,2))
 resid <- g$residuals
 qqnorm(resid)
 qqline(resid)

 plot(c(treatments),resid)  # Use c(treaments) so 'treatments' is once agina viewed as a numberiacl variable; 
                            # otherwise you get a boxplot
 abline(h=0)  # What does this do?
 plot(c(blocks),resid)
 abline(h=0)
 fits <- g$fitted.values
 plot(fits,resid)
 abline(h=0)


 bartlett.test(y,treatments)

 bartlett.test(y,blocks)

 

# Levene's test for equality of variances within blocks:
# Make the data into a matrix, rather than a frame, then 
# compute the abs. value of differences between the data 
# and the block (row) medians

 data.matrix <- as.matrix(data)
 abs.diff <- matrix(nrow=4, ncol=4)
 for(i in 1:4) {for (j in 1:4) abs.diff[i,j] <- abs(data.matrix[i,j] - median(data.matrix[ ,j]))}

 abs.diff
  g.blocks <- lm(c(t(abs.diff))~blocks)  
       # Transpose the matrix 'abs.diff', 
       # then string out the columns as one long column
  anova(g.blocks) 

# Friedman's test depends on the ranks within each block:

 data.ranks <- apply(data,2,rank)
 data.ranks

# Note that tip 4 has the largest rank each time, indicating that 
# its mean effect is greatest.

  friedman.test(y,treatments,blocks)