############Randomized complete block design############
# # page 895, in an experiment on decision making, executives
# # were exposed to one of three methods of quantifying the
# # maximum risk premium they would be willing to pay to
# avoid uncertaint in a business decision. The three methods
# are the utility method, the worry method, and the comparison method.
# After using the assigned methods, the subjects were asked
# to state their degree of confidence in the method of
# quantifying the risk premium on a scale from 0 (no confidence)
# to 20 (highest confidence).
# Fifteen subjects are used in the study. They were grouped
# into five blocks of three excutives, according to age.
# Block 1 contained the three oldest excutives, and so on.
# Five independent random permutations of three were
# employed within each block.
riskp<-read.table("~/Desktop/jenn/teaching/stat445545/data/CH21TA01.txt",col.names=c("y","block","method"))
nt<-nrow(riskp)
nt
## [1] 15
attach(riskp)
riskp
## y block method
## 1 1 1 1
## 2 5 1 2
## 3 8 1 3
## 4 2 2 1
## 5 8 2 2
## 6 14 2 3
## 7 7 3 1
## 8 9 3 2
## 9 16 3 3
## 10 6 4 1
## 11 13 4 2
## 12 18 4 3
## 13 12 5 1
## 14 14 5 2
## 15 17 5 3
#summary statistics
riskp.xtabs <-xtabs(y~method+block)
riskp.xtabs
## block
## method 1 2 3 4 5
## 1 1 2 7 6 12
## 2 5 8 9 13 14
## 3 8 14 16 18 17
tapply(y,method, data=riskp, mean) #row means
## 1 2 3
## 5.6 9.8 14.6
tapply(y,block, data=riskp, mean) #column means
## 1 2 3 4 5
## 4.666667 8.000000 10.666667 12.333333 14.333333
##View interaction plots, also called profile plots
interaction.plot(method,block,y,type='b',
col=1:5, pch=1:5)
##fit anova model (treat block as fixed)
myfit <- aov(y~factor(block) +factor(method),data=riskp)
summary(myfit)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(block) 4 171.33 42.83 14.36 0.001008 **
## factor(method) 2 202.80 101.40 33.99 0.000123 ***
## Residuals 8 23.87 2.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tukey test for additivity
#install.packages("asbio")
library(asbio)
## Loading required package: tcltk
tukey.add.test(y, block, method)
##
## Tukey's one df test for additivity
## F = 0.0778959 Denom df = 7 p-value = 0.7882351
##residual analysis
myfit$resid
## 1 2 3 4 5 6
## 0.7333333 0.5333333 -1.2666667 -1.6000000 0.2000000 1.4000000
## 7 8 9 10 11 12
## 0.7333333 -1.4666667 0.7333333 -1.9333333 0.8666667 1.0666667
## 13 14 15
## 2.0666667 -0.1333333 -1.9333333
qqnorm(myfit$resid)
qqline(myfit$resid)
plot(myfit$fitted,myfit$res,xlab="Fitted",ylab="Residuals")
