#TWO-WAY ANOVA LINEAR MODELS
require(car)

#CELL MEANS MODEL IN R
#READ STRUCTURED DATA TABLE WITH CELL MEANS CODING
K=read.table("c:/2008LinearModelsData/KNNL23.8R.txt")
K
attach(K)
Y=Y
X=cbind(K$X1,K$X2,K$X3,K$X4,K$X5,K$X6)
X

FMc=lm(Y~X-1)
FMC=lm(Y~X1+X2+X3+X4+X5+X6-1)
summary(FMc)  #CELL MEANS AS ESTIMATES
summary(FMC)  #CELL MEANS AS ESTIMATES

#CENTERING CELL MEANS MODEL
Y=scale(Y,scale=F)
FMcs=lm(Y~X-1)
FMCs=lm(Y~X1+X2+X3+X4+X5+X6-1)

anova(FMcs) #TYPE 1 ANOVA SS KNNL CELL MEANS MODEL FOR LUMPED X VARIABLE
anova(FMCs) #TYPE 1 ANOVA SS KNNL CELL MEANS MODEL ITEMIZED FOR EACH X VARIABLE
Anova(FMCs) #TYPE 2 ANOVA SS KNNL CELL MEANS MODEL ITEMIZED...
Anova(FMCs,type="3") #TYPE 3 ANOVA SS KNNL CELL MEANS MODEL ITEMIZED...
drop1(FMCs)
detach(K)

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

#TREATMENTS MODEL IN R
#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)

#CROSS TABULATIONS VERIFYING CELL MEANS
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

FM=lm(Y~A*B)
summary(FM)
MM=data.frame(model.matrix(FM))
MM

anova(FM)  #TYPE 1 ANOVA SS TREATMENTS MODEL IN R
Anova(FM)  #TYPE 2 ANOVA SS TREATMENTS MODEL IN R
Anova(FM,type="III") #TYPE 3 ANOVA SS TREATMENTS MODEL IN R

#GLM TESTS FOR TREATMENTS MODEL IN R
attach(MM)
RMi=lm(Y~A2+B2+B3)
RMa=lm(Y~B2+B3+A2.B2+A2.B3)
RMb=lm(Y~A2+A2.B2+A2.B3)
anova(RMi,FM) #FOR INTERACTIONS
anova(RMa,FM) #FOR FACTOR A
anova(RMb,FM) #FOR FACTOR B   #NOTE DIFFERENT RESULT HERE vs KNNL

#EXPLICIT STATEMENT OF FM TO ITEMIZE FACTORS
FM2=lm(Y~A2+B2+B3+A2.B2+A2.B3)
Anova(FM2,type="III")
drop1(FM2)

Anova(RMi,type="3")
anova(RMi,FM2)

Anova(RMa,type="3")
anova(RMa,FM2)

Anova(FM2,type="3")
anova(RMb,FM2)
detach(MM)

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

#CONVERTING TO KNNL TREATMENT EFFECTS MODEL = contr.sum MODEL IN R
#FULL MODEL FROM contr.sum
contrasts(A)=contr.sum
contrasts(B)=contr.sum
FM=lm(Y~A*B)
summary(FM)
MM=data.frame(model.matrix(FM))
MM

anova(FM)  #TYPE 1 ANOVA SS KNNL TREATMENT EFFECTS MODEL
Anova(FM)  #TYPE 2 ANOVA SS KNNL TREATMENT EFFECTS MODEL
Anova(FM,type="III") #TYPE 3 ANOVA SS KNNL TREATMENT EFFECTS MODEL

#GLM TESTS FOR KNNL TREATMENTS EFFECTS MODEL = contr.sum in R

attach(MM)
RMi=lm(Y~A1+B1+B2)
RMa=lm(Y~B1+B2+A1.B1+A1.B2)
RMb=lm(Y~A1+A1.B1+A1.B2)
anova(RMi,FM) #FOR INTERACTIONS KNNL TREATMENT EFFECTS MODEL 
anova(RMa,FM) #FOR FACTOR A KNNL TREATMENT EFFECTS MODEL
anova(RMb,FM) #FOR FACTOR B KNNL TREATMENT EFFECTS MODEL = KNNL RESULTS

#EXPLICIT STATEMENT OF FM TO ITEMIZE FACTORS
FM2=lm(Y~A1+B1+B2+A1.B1+A1.B2)
Anova(FM2,type="III")
drop1(FM2)

Anova(RMi,type="3")
anova(RMi,FM2)

Anova(RMa,type="3")
anova(RMa,FM2)

Anova(FM2,type="3")
anova(RMb,FM2)
detach(MM)

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

#PAIRWISE CONTRASTS OF FACTOR LEVEL MEANS

N=length(Y)
attach(K)

#EXTRACTING DATA BLOCKS
Y11=Rate[GenderA=="1"&BoneDevB=="1"]
Y12=Rate[GenderA=="1"&BoneDevB=="2"]
Y13=Rate[GenderA=="1"&BoneDevB=="3"]
Y21=Rate[GenderA=="2"&BoneDevB=="1"]
Y22=Rate[GenderA=="2"&BoneDevB=="2"]
Y23=Rate[GenderA=="2"&BoneDevB=="3"]

#BLOCK MEANS
mean(Y11)
mean(Y12)
mean(Y13)
mean(Y21)
mean(Y22)
mean(Y23)

#UNWEIGHTED TREATMENT MEANS FOR FACTOR B
mu1=mean(c(mean(Y11),mean(Y21)))
mu2=mean(c(mean(Y12),mean(Y22)))
mu3=mean(c(mean(Y13),mean(Y23)))
mu1
mu2
mu3

#PAIRWISE CONTRASTS FOR FACTOR B
pairwise.t.test(Y,B,"none")  #SINGLE P
pairwise.t.test(Y,B,"bonferroni") #BONFERRONI P
pairwise.t.test(Y,B,"holm")  # A USEFUL BONFERRONI VARIANT

#CONFIDENCE INTERVALS FOR 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

#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 TEST
a
b
#NUMBERS in CLASSES
n13=length(Y13)
n23=length(Y23)

alpha=0.05 #TEST LEVEL

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

V3=(MSE/(a^2))*sum((1/n13)+(1/n23))
V3  #VARIANCE OF mu3
sd3=sqrt(V3)
sd3 #STANDARD DEVIATION OF mu3
t=mu3/sd3
t  #T STATISTIC
qt(1-alpha,N-a*b) #CRITICAL VALUE
1-pt(t,N-a*b) #PROBABILITY