# This example shows the analysis for the RBD with sampling  
# using the rubber 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, "                                                                                       
 1 1 A 72 
 2 1 A 79 
 3 1 A 61 
 4 1 A 71 
 5 1 B 133 
 6 1 B 129 
 7 1 B 123 
 8 1 B 156 
 9 1 C 37 
 10 1 C 36 
 11 1 C 26 
 12 1 C 24 
 13 1 D 63 
 14 1 D 49 
 15 1 D 63 
 16 1 D 43 
 17 1 E 35 
 18 1 E 26 
 19 1 E 24 
 20 1 E 61 
 21 1 F 31 
 22 1 F 32 
 23 1 F 28 
 24 1 F 26 
 25 1 G 43 
 26 1 G 40 
 27 1 G 35 
 28 1 G 38 
 29 2 A 61 
 30 2 A 49 
 31 2 A 57 
 32 2 A 61 
 33 2 B 129 
 34 2 B 125 
 35 2 B 136 
 36 2 B 127 
 37 2 C 20 
 38 2 C 14 
 39 2 C 30 
 40 2 C 27 
 41 2 D 51 
 42 2 D 52 
 43 2 D 62 
 44 2 D 52 
 45 2 E 27 
 46 2 E 27 
 47 2 E 26 
 48 2 E 26 
 49 2 F 22 
 50 2 F 20 
 51 2 F 29 
 52 2 F 28.0 
 53 2 G 32.0 
 54 2 G 29.0 
 55 2 G 45.0 
 56 2 G 40.0 
 57 3 A 70.0 
 58 3 A 62.0 
 59 3 A 62.0 
 60 3 A 76.0 
 61 3 B 121.0 
 62 3 B 125.0 
 63 3 B 109.0 
 64 3 B 128.0 
 65 3 C 33.0 
 66 3 C 33.0 
 67 3 C 27.0 
 68 3 C 29.5 
 69 3 D 58.0 
 70 3 D 64.0 
 71 3 D 56.0 
 72 3 D 55.0 
 73 3 E 28.0 
 74 3 E 28.0 
 75 3 E 27.0 
 76 3 E 29.0 
 77 3 F 27.0 
 78 3 F 30.5 
 79 3 F 27.0 
 80 3 F 27.0 
 81 3 G 44.0 
 82 3 G 44.0 
 83 3 G 45.0 
 84 3 G 49.0 
 85 4 A 36.0 
 86 4 A 39.0 
 87 4 A 41.0 
 88 4 A 45.0 
 89 4 B 57.0 
 90 4 B 58.0 
 91 4 B 59.0 
 92 4 B 67.0 
 93 4 C 27.0 
 94 4 C 24.0 
 95 4 C 22.0 
 96 4 C 25.0 
 97 4 D 38.0 
 98 4 D 38.0 
 99 4 D 37.0 
 100 4 D 38.0 
 101 4 E 22.0 
 102 4 E 23.0 
 103 4 E 20 
 104 4 E 20 
 105 4 F 22 
 106 4 F 23 
 107 4 F 22 
 108 4 F 22 
 109 4 G 31 
 110 4 G 31 
 111 4 G 28 
 112 4 G 30 
 113 5 A 58 
 114 5 A 57 
 115 5 A 58 
 116 5 A 53 
 117 5 B 122 
 118 5 B 98 
 119 5 B 107 
 120 5 B 110 
 121 5 C 34 
 122 5 C 27 
 123 5 C 26 
 124 5 C 26 
 125 5 D 53 
 126 5 D 47 
 127 5 D 48 
 128 5 D 47 
 129 5 E 25 
 130 5 E 25 
 131 5 E 21 
 132 5 E 19 
 133 5 F 26 
 134 5 F 25 
 135 5 F 22 
 136 5 F 18 
 137 5 G 43 
 138 5 G 35 
 139 5 G 43 
 140 5 G 36 
 141 6 A 52 
 142 6 A 56 
 143 6 A 52 
 144 6 A 50 
 145 6 B 109 
 146 6 B 120 
 147 6 B 112 
 148 6 B 107 
 149 6 C 30 
 150 6 C 31 
 151 6 C 31 
 152 6 C 28 
 153 6 D 50 
 154 6 D 50.0 
 155 6 D 50.0 
 156 6 D 51.0 
 157 6 E 25.0 
 158 6 E 25.0 
 159 6 E 26.0 
 160 6 E 26.0 
 161 6 F 24.0 
 162 6 F 26.0 
 163 6 F 25.0 
 164 6 F 26.0 
 165 6 G 38.0 
 166 6 G 41.0 
 167 6 G 40.0 
 168 6 G 43.0 
 169 7 A 40.7 
 170 7 A 45.9 
 171 7 A 43.1 
 172 7 A 37.3 
 173 7 B 80.0 
 174 7 B 71.9 
 175 7 B 75.8 
 176 7 B 63.7 
 177 7 C 26.5 
 178 7 C 27.1 
 179 7 C 26.6 
 180 7 C 25.6 
 181 7 D 38.8 
 182 7 D 39.4 
 183 7 D 40.7 
 184 7 D 38.0 
 185 7 E 23.0 
 186 7 E 22.9 
 187 7 E 22.5 
 188 7 E 35.7 
 189 7 F 22.2 
 190 7 F 23.9 
 191 7 F 22.6 
 192 7 F 25.5 
 193 7 G 29.4 
 194 7 G 31.6 
 195 7 G 29.6 
 196 7 G 29.3 
 197 8 A 68.1 
 198 8 A 69.8 
 199 8 A 65.9 
 200 8 A 62.1 
 201 8 B 135.0 
 202 8 B 151.0 
 203 8 B 143.0 
 204 8 B 142.0 
 205 8 C 38.1 
 206 8 C 37.4 
 207 8 C 37.9 
 208 8 C 37.1 
 209 8 D 64.5 
 210 8 D 65.7 
 211 8 D 64.0 
 212 8 D 62.5 
 213 8 E 32.1 
 214 8 E 35.2 
 215 8 E 33.0 
 216 8 E 34.9 
 217 8 F 32.7 
 218 8 F 32.4 
 219 8 F 30.3 
 220 8 F 35.6 
 221 8 G 50.2 
 222 8 G 50.4 
 223 8 G 42.5 
 224 8 G 45.0 
 225 9 A 46.0 
 226 9 A 47.0 
 227 9 A 46.0 
 228 9 A 45.0 
 229 9 B 69.0 
 230 9 B 69.0 
 231 9 B 73.0 
 232 9 B 70.0 
 233 9 C 26.0 
 234 9 C 26.0 
 235 9 C 25.0 
 236 9 C 25.0 
 237 9 D 40.0 
 238 9 D 38.0 
 239 9 D 39.0 
 240 9 D 39.0 
 241 9 E 24.0 
 242 9 E 24.0 
 243 9 E 24.0 
 244 9 E 25.0 
 245 9 F 23.0 
 246 9 F 24.0 
 247 9 F 24.0 
 248 9 F 23.0 
 249 9 G 32.0 
 250 9 G 31.0 
 251 9 G 32.0 
 252 9 G 30.0 
 253 10 A 77.0 
 254 10 A 74.0 
 255 10 A 77.0 
 256 10 A 72 
 257 10 B 132 
 258 10 B 129 
 259 10 B 141 
 260 10 B 137 
 261 10 C 45 
 262 10 C 41 
 263 10 C 39 
 264 10 C 38 
 265 10 D 71 
 266 10 D 69 
 267 10 D 66 
 268 10 D 68 
 269 10 E 36 
 270 10 E 33 
 271 10 E 35 
 272 10 E 25 
 273 10 F 38 
 274 10 F 36 
 275 10 F 38 
 276 10 F 38 
 277 10 G 56 
 278 10 G 48 
 279 10 G 48 
 280 10 G 50 
 281 11 A 76 
 282 11 A 55 
 283 11 A 60 
 284 11 A 58 
 285 11 B 118 
 286 11 B 109 
 287 11 B 115 
 288 11 B 106 
 289 11 C 27 
 290 11 C 32 
 291 11 C 26 
 292 11 C 26 
 293 11 D 52 
 294 11 D 45 
 295 11 D 48 
 296 11 D 54 
 297 11 E 22 
 298 11 E 19 
 299 11 E 18 
 300 11 E 23 
 301 11 F 23 
 302 11 F 23 
 303 11 F 23 
 304 11 F 24 
 305 11 G 32 
 306 11 G 37 
 307 11 G 37.0 
 308 11 G 39.0 
 309 12 A 72.5 
 310 12 A 76.0 
 311 12 A 69.5 
 312 12 A 70.5 
 313 12 B 133.0 
 314 12 B 133.0 
 315 12 B 128.5 
 316 12 B 128.5 
 317 12 C 32.5 
 318 12 C 32.8 
 319 12 C 32.9 
 320 12 C 34.6 
 321 12 D 63.0 
 322 12 D 64.5 
 323 12 D 61.5 
 324 12 D 62.7 
 325 12 E 31.2 
 326 12 E 30.2 
 327 12 E 29.0 
 328 12 E 29.7 
 329 12 F 30.7 
 330 12 F 30.8 
 331 12 F 30.0 
 332 12 F 29.5 
 333 12 G 45.8 
 334 12 G 45.2 
 335 12 G 43.5 
 336 12 G 46.5 
 337 13 A 51.0 
 338 13 A 50.0 
 339 13 A 49.0 
 340 13 A 49.0 
 341 13 B 86.0 
 342 13 B 84.0 
 343 13 B 96.0 
 344 13 B 81.0 
 345 13 C 24.0 
 346 13 C 24.0 
 347 13 C 24.0 
 348 13 C 26.0 
 349 13 D 45.0 
 350 13 D 43.0 
 351 13 D 42.0 
 352 13 D 45.0 
 353 13 E 21.8 
 354 13 E 21.8 
 355 13 E 24.0 
 356 13 E 22.0 
 357 13 F 24.0 
 358 13 F 24 
 359 13 F 22 
 360 13 F 24 
 361 13 G 33 
 362 13 G 33 
 363 13 G 31 
 364 13 G 31
", sep=" ")

options(scipen=999) # suppressing scientific notation

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

attach(rubber)

# The data frame called rubber is now created, 
# with four variables, obs, lab, material, and stretching.
##
#########
####################################################################

############################################################################### 
# 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 lab and material are factors with the factor() function:

lab <- factor(lab)
material <- factor(material) 

rubber.RBDsamp.fit <- lm(stretching ~ lab + material + lab:material);
anova(rubber.RBDsamp.fit)
   
# stretching is the response variable, material is the treatment factor,  
# and lab represents the blocks. 

# We MUST specify the lab:material interaction in the lm statement!!!    
# This will serve as our "experimental error" term.                         
# The regular "Residuals" (i.e., "Error") line in R will serve as 
# our "sampling error" term.  


# This code performs the F-test for a MATERIAL (treatment) effect, where
# the correct denominator involves the LAB*MATERIAL (experimental error) term:

MS.Trt <- anova(rubber.RBDsamp.fit)["material","Mean Sq"]
MS.Exper.Error <- anova(rubber.RBDsamp.fit)["lab:material","Mean Sq"]
df.Trt <- anova(rubber.RBDsamp.fit)["material","Df"]
df.Exper.Error <- anova(rubber.RBDsamp.fit)["lab:material","Df"]

F.statistic <- MS.Trt/MS.Exper.Error; print(paste("F = ", round(F.statistic,4)))
P.Value <- pf(F.statistic, df1=df.Trt, df2=df.Exper.Error, lower.tail=F); print(paste("P-value of F-test = ", round(P.Value,4)))

# Clearly there is a significant MATERIAL effect (F = 135.48, P-value < .0001). 

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

# The sample mean stretching amounts for each material, listed from smallest to largest:

sort( tapply(stretching, material, mean) )

# Now, which of these means are significantly different?

# Some code for the
# Tukey procedure:

# Specifying the experimentwise significance level:
alpha <- 0.05

df.Blk <- anova(rubber.RBDsamp.fit)["lab","Df"]

n <- 4 # (Since there are 4 observations per treatment-block combination here!)


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

my.std.error <- sqrt(2*MS.Exper.Error/((df.Blk+1)*n))

my.T <- qtukey(1-alpha, df.Trt+1, df.Trt*df.Blk) / 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) 

# The CIs that DO NOT contain zero indicate that those pairs of means are significantly different.