# Tutorial 5


# Q1
x1 = c(3,2,2.5,3,3,3.5)
x2=c(5.5,4,5,4.5,4,4.5)
x3=c(5,4.5,4,4.5,2.5,4.5)
x4=c(6,5.5,5,6,5.5,5.5)

r1 =c(1,1,1,1,2,1)
r2= c(3,2,3.5,2.5,3,2.5)
r3=c(2,3,2,2.5,1,2.5)
r4=c(4,4,3.5,4,4,4)

Rk = cbind(r1,r2,r3,r4)

Rmean=apply(Rk, 2, mean)

rbind(Rk, Rmean)

g=4
n=6 

Xsq = 12*n/g/(g+1)*sum((Rmean-(g+1)/2)^2)

1-pchisq(Xsq,3)

#Q2 

Imi = c(6,4,6,7,5,6,8,7,8,3,9,4,8,11,12,6,10,3,5,4,6,7,5,6,3,10,5,4,4,7)
Pla = c(4,7,12,10,2,11,9,5,11,8,7,6,8,9,9,8,10,9,8,5,8,7,6,9,3,5,11,7,3,10)

d=Imi-Pla

y=c(Imi, Pla)

y.m=mean(y)
Imi.m=mean(Imi)
Pla.m=mean(Pla)

X = cbind(Imi,Pla)
b.mean = apply(X,1,mean)

TSS = 30*((Imi.m-y.m)^2 + (Pla.m-y.m)^2)
BSS = 2*sum((b.mean-y.m)^2)

X1 = b.mean%*%t(c(1,1))
X2 = rep(1,30)%*%t(c(Imi.m,Pla.m))
X3 = rep(1,30)%*%t(c(1,1))*y.m

RSS=sum(apply((X-X1-X2+X3)^2,2,sum))

# another way to compute RSS
TTSS = sum((y-y.m)^2)
RSS = TTSS-TSS-BSS

# 
c(TSS,BSS,RSS)/29

F= TSS/RSS*29 

d.m = mean(d)

s.d = sqrt(var(d))

(d.m/s.d * sqrt(30) )^2


#Q3

y1 = c(8.4,12.8,9.6,9.8,8.4,8.6,8.9,7.9)
y2 = c(9.4,15.2,9.1,8.8,8.2,9.9,9.0,8.1)
y3 = c(9.8,12.9,11.2,9.9,8.5,9.8,9.2,8.2)
y4 = c(12.2,14.4,9.8,12.0,8.5,10.9,10.4,10)
y=c(y1,y2,y3,y4)

tmt=factor(c(rep(1,8),rep(2,8),rep(3,8),rep(4,8)))
blk=factor(rep(1:8,4))
X=data.frame(blk,tmt)

options(contrasts=c("contr.treatment","contr.poly"))
attach(X)


lm.fit= lm(y~tmt+blk,data=X) 

anova(lm.fit)

vv=vcov(lm.fit)
V=vv[2:4,2:4]
bb=lm.fit$coef[2:4]

c1=c(0,1,0)
c2=c(1,-1/2,0)
c3=c(1/3,1/3,-1)
C1 =  t(c1)%*%bb                
C2 =  t(c2)%*%bb                
C3 =  t(c3)%*%bb                

v1 = V[2,2]
v2 = t(c2)%*%V%*%c2
v3 = t(c3)%*%V%*%c3

L1 = C1/sqrt(v1)
L2 = C2/sqrt(v2)
L3 = C3/sqrt(v3)

c(L1,L2,L3)

scheffe = sqrt(3*qf(0.95,3,21))
bonferroni = qt(1-0.05/6, 21)