#INTERVAL ESTIMATION IN 2-WAY ANOVA LINEAR MODEL
require(car)      #MUST LOAD {gmodels} PACKAGE FROM CRAN
require(gmodels)  #MUST LOAD {gmodels} PACKAGE FROM CRAN

#READ STRUCTURED DATA TABLE WITH NUMERIC CODED FACTOR
K=read.table("c:/2008LinearModelsData/GrowthHormoneR.txt")
K
attach(K)
Y=Rate
A=factor(GenderA)
B=factor(BoneDevB)
detach(K)

#CONVERTING TO TREATMENT EFFECTS = contr.sum MODEL IN R
#FULL MODEL FROM contr.sum
contrasts(A)=contr.sum
contrasts(B)=contr.sum
FM=lm(Y~A*B)

MM=data.frame(model.matrix(FM))
MM

summary(FM)
Anova(FM,type="III")  #TYPE 3 ANOVA SS

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#CROSS TABULATIONS FOR CELL MEANS AND COUNTS IN EACH BLOCK
XX=xtabs(Y~A+B) 
XX    #SUM OF Y VALUES IN EACH BLOCK
n=table(A,B)  
n     #COUNTS FOR EACH BLOCK
N=length(Y)
N     #TOTAL NUMBER OF OBJECTS
mu=XX/n   
mu    #MEANS FOR EACH BLOCK
N=length(Y)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#PAIRWISE CONTRASTS OF FACTOR LEVEL MEANS

#CONFIDENCE INTERVALS FOR SINGLE CONTRASTS

CS12=fit.contrast(FM,B,c(1,-1,0),conf.int=0.90)
CS13=fit.contrast(FM,B,c(1,0,-1),conf.int=0.90)
CS23=fit.contrast(FM,B,c(0,1,-1),conf.int=0.90)

CS12
CS13
CS23

#PAIRWISE.T.TEST GIVES WRONG ANSWERS
pairwise.t.test(Y,B,"none")

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#TUKEY CI'S
alpha=0.10
a=nlevels(A)
a
b=nlevels(B)
b

Q=qtukey((1-alpha),(b),(N-a*b)) # 3 CONTRASTS KNNL P961
T=(1/sqrt(2))*Q

D=CS12[1]
lwr=CS12[1]-CS12[2]*T
upr=CS12[1]+CS12[2]*T
CT12=cbind(lwr,D,upr)

D=CS13[1]
lwr=CS13[1]-CS13[2]*T
upr=CS13[1]+CS13[2]*T
CT13=cbind(lwr,D,upr)

D=CS23[1]
lwr=CS23[1]-CS23[2]*T
upr=CS23[1]+CS23[2]*T
CT23=cbind(lwr,D,upr)

CT12
CT13
CT23

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

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

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

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

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

CB12
CB13
CB23

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

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

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

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

CF12
CF13
CF23

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#SINGLE DEGREE OF FREEDOM TESTS
a
b

alpha=0.05 #TEST LEVEL

MSE=anova(FM)[4,3]
MSE  #FROM anova()

#VARIANCES
V1=(MSE/(a^2))*sum((1/n[1,1])+(1/n[2,1]))
V2=(MSE/(a^2))*sum((1/n[1,2])+(1/n[2,2]))
V3=(MSE/(a^2))*sum((1/n[1,3])+(1/n[2,3]))
VAR=c(V1,V2,V3)
#STANDARD ERROR
sd1=sqrt(V1)
sd2=sqrt(V2)
sd3=sqrt(V3)
SE=c(sd1,sd2,sd3)
#MEANS
mu_dot1=mean(c(mu[1,1],mu[2,1]))
mu_dot2=mean(c(mu[1,2],mu[2,2]))
mu_dot3=mean(c(mu[1,3],mu[2,3]))
MEAN=c(mu_dot1,mu_dot2,mu_dot3)
# t STATISTICS
t1=mu_dot1/sd1
t2=mu_dot2/sd2
t3=mu_dot3/sd3
tSTAT=c(t1,t2,t3)
#CRITICAL VALUES
CV1=qt(1-alpha/2,N-a*b) #TWO-SIDED
CV2=qt(1-alpha/2,N-a*b) #TWO-SIDED
CV3=qt(1-alpha,N-a*b)   #ONE-SIDED
CV=c(CV1,CV2,CV3)
#PROBABILITIES
P1=1-pt(t1,N-a*b) #TWO-SIDED TEST
P2=1-pt(t2,N-a*b) #TWO-SIDED TEST
P3=1-pt(t3,N-a*b) #ONE-SIDED TEST
PROB=c(P1,P2,P3)
RESULT=cbind(VAR,SE,MEAN,tSTAT,CV,PROB)
RESULT