Tutorial 7

Q1

x1 = c(130,155,74,180,34,40,80,75,20,70,82,58)
x2 = c(150,188,159,126,126,122,106,115,25,70,58,45)
x3 = c(138,110,168,160,174,120,150,139,96,104,82,60)

x=c(x1,x2,x3)
S = factor(c(rep(1,12),rep(2,12),rep(3,12)))
T = factor(rep(c(rep(1,4),rep(2,4),rep(3,4)),3))

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

lm.fit=lm(x~T*S)
anova(lm.fit)

V=vcov(lm.fit)[-1,-1]
b = lm.fit$coef[-1]
c11 = c(1,0,0,0,0,0,0,0)
c12 = c(0,1,0,0,0,0,0,0)
c13 = c11-c12
c21 = c(1,0,0,0,1,0,0,0)
c22 = c(0,1,0,0,0,1,0,0)
c23 = c21-c22
c31 = c(1,0,0,0,0,0,1,0)
c32 = c(0,1,0,0,0,0,0,1)
c33=c31-c32
C = rbind(c11,c12,c13,c21,c22,c23,c31,c32,c33)
l = C%*%b
s = sqrt(diag(C%*%V%*%t(C)))
L = l/s

round(t(L),3)

c.value = sqrt(2*qf(0.95,2,31))


lm.fit2=lm(x~T+S)
b=lm.fit2$coef[-c(1,4,5)]
V = vcov(lm.fit2)[-c(1,4,5),-c(1,4,5)]
c1 = c(1,0)
c2 = c(0,1)
c3 = c(1/2,1/2)
C=rbind(c1,c2,c3)
l = C%*%b
s = sqrt(diag(C%*%V%*%t(C)))
L = l/s

anova(lm.fit2)
round(t(L),3)

Q3 

n=c(19, 20, 20)
x.m=c(5.3158, 8.3000, 8.5500)
s.x=c(67.1053, 76.2000, 139.9500)
z.m=c(2.3158, 2.4500, 3.1500)
s.z =c(23.1050, 15.9500, 22.5500)
s.xz=c(27.1053, 18.3000, 34.3500)
b = s.xz/s.z
mu = x.m-b*z.m

s2 = (s.x - s.xz^2/s.z)/(n-2) 
ss2 = sum(s2*(n-2))/sum(n-2)
V.b = ss2*diag(1/s.z)
c1 = c(1,-1,0)
c2 = c(0,1,-1)
C = rbind(c1,c2)
F = t(C%*%b)%*%solve(C%*%V.b%*%t(C))%*% C%*%b

bb = sum(b*s.z)/sum(s.z)
mmu = x.m-bb*z.m

nn=sum(n)
ss.x=sum(s.x)
ss.z=sum(s.z)
ss.xz =sum(s.xz)
ss2=(ss.x-ss.xz^2/ss.z)/(nn-3-1)
S = ss2*(z.m%*%t(z.m)/ss.z + diag(1/n) )
c1 = c(1,-1,0)
c2=c(1,0,-1)
c3 =c(1,-1/2,-1/2)
C=rbind(c1,c2,c3)
l=C%*%mmu 
s = sqrt(diag(C%*%S%*%t(C)))
L=l/s
round(t(L),3)

sqrt(2*qf(0.95,2,55))