TREATMENTS MODEL IN R
#READ STRUCTURED DATA TABLE WITH NUMERIC CODED FACTOR
K=read.table("c:/2008LinearModelsData/KentonFoodR.txt")
K
attach(K)

Y=Sales
X=factor(Design) # factor() IN DEFAULT SETTING

FM=lm(Y~X)
model.matrix(FM)
summary(FM)
anova(FM)

MSE=summary(FM)$sigma^2
MSE

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
aov(FM)

#TUKEYHSD MULTIPLE COMPARISONS
TukeyHSD(aov(FM),conf.level=0.90)

#STUDENTIZED RANGE DISTRIBUTION FUNCTIONS
?ptukey
?qtukey
alpha=0.10
N=length(Y) #TOTAL NUMBER OF CASES
levels(X) 
r=4         #NUMBER OF LEVELS IN FACTOR X
Q=qtukey(1-alpha,r,N-r)
Q

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
?pairwise.t.test()

#PAIRWISE COMPARISONS USING BONFERRONI & HOLM
pairwise.t.test(Y,X,"none")  #SINGLE P
pairwise.t.test(Y,X,"bonferroni") #BONFERRONI P
pairwise.t.test(Y,X,"holm")  # A USEFUL BONFERRONI VARIANT

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#MAKING CONTRASTS IN R
require(gmodels)  #MUST LOAD {gmodels} PACKAGE FROM CRAN

#CONTRAST FOR KNNL P. 743T
fit.contrast(FM,X,c(0.5,0.5,-0.5,-0.5),conf.int=0.95)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#SINGLE CALCULATION OF SCHEFFE 

alpha=.10  #SET FAMILYWIDE CONFIDENCE ALPHA

N=length(Y) #TOTAL NUMBER OF CASES
levels(X) 
r=4         #NUMBER OF LEVELS IN FACTOR X
S=sqrt((r-1)*qf(1-alpha,r-1,N-r)) #SCHEFFE MULTIPLIER
S
#CALCULATING NECESSARY VALUES FROM CONTRAST
LFIT=fit.contrast(FM,X,c(0.5,0.5,-0.5,-0.5),conf.int=0.90)
LFIT

L=LFIT[1] # CALCULATION OF LINEAR CONTRAST L
L
sL=LFIT[2] # CALCULATION OF sL
sL
CI=LFIT[5:6]  # SINGLE CONFIDENCE INTERVAL
CI
CIScheffelower=L-S*sL
CIScheffeupper=L+S*sL
CIS=cbind(CIScheffelower,L,CIScheffeupper)
CIS  #SCHEFFE FAMILY CONFIDENCE INTERVAL

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#CONVERTING P VALUES IN pairwise.t.test() INTO CI'S
options(digits=10)
alpha=0.10
#SINGLE CI'S
C=abs(qt(alpha/2,N-r)) #SINGLE MULTIPLIER
C
summary(FM) #CALCULATE BY HAND: coef +/- C*Std.Error

CS12=fit.contrast(FM,X,c(1,-1,0,0),conf.int=0.90)
CS13=fit.contrast(FM,X,c(1,0,-1,0),conf.int=0.90)
CS14=fit.contrast(FM,X,c(1,0,0,-1),conf.int=0.90)
CS23=fit.contrast(FM,X,c(0,1,-1,0),conf.int=0.90)
CS24=fit.contrast(FM,X,c(0,1,0,-1),conf.int=0.90)
CS34=fit.contrast(FM,X,c(0,0,1,-1),conf.int=0.90)
CS12
CS13
CS14
CS23
CS24
CS34

#BONFERRONI CI'S
alpha=0.10
g=6  #NUMBER OF PAIRWISE CONTRASTS
B=qt(1-(alpha/(2*g)),N-r) #BONFERRONI MULTIPLIER
B

lwr=CS12[1]-CS12[2]*B
upr=CS12[1]+CS12[2]*B
CB12=cbind(lwr,upr)

lwr=CS13[1]-CS13[2]*B
upr=CS13[1]+CS13[2]*B
CB13=cbind(lwr,upr)

lwr=CS14[1]-CS14[2]*B
upr=CS14[1]+CS14[2]*B
CB14=cbind(lwr,upr)

lwr=CS23[1]-CS23[2]*B
upr=CS23[1]+CS23[2]*B
CB23=cbind(lwr,upr)

lwr=CS24[1]-CS24[2]*B
upr=CS24[1]+CS24[2]*B
CB24=cbind(lwr,upr)

lwr=CS34[1]-CS34[2]*B
upr=CS34[1]+CS34[2]*B
CB34=cbind(lwr,upr)

CB12
CB13
CB14
CB23
CB24
CB34

#SCHEFFE CI'S
alpha=0.10
S=sqrt((r-1)*qf(1-alpha,r-1,N-r)) #SCHEFFE MULTIPLIER
S

lwr=CS12[1]-CS12[2]*S
upr=CS12[1]+CS12[2]*S
CSF12=cbind(lwr,upr)

lwr=CS13[1]-CS13[2]*S
upr=CS13[1]+CS13[2]*S
CSF13=cbind(lwr,upr)

lwr=CS14[1]-CS14[2]*S
upr=CS14[1]+CS14[2]*S
CSF14=cbind(lwr,upr)

lwr=CS23[1]-CS23[2]*S
upr=CS23[1]+CS23[2]*S
CSF23=cbind(lwr,upr)

lwr=CS24[1]-CS24[2]*S
upr=CS24[1]+CS24[2]*S
CSF24=cbind(lwr,upr)

lwr=CS34[1]-CS34[2]*S
upr=CS34[1]+CS34[2]*S
CSF34=cbind(lwr,upr)

CSF12
CSF13
CSF14
CSF23
CSF24
CSF34