# R-code for E-learning Week notes:

# Prepare the data
x11 = c(30.9,31.9,31.3,32.1,30.9,31.3,31.3,32.1,30.3,32.2)
x12 = c(30.7,31.6,31.1,31.0,31.2,31.7,31.8,33.0,30.9,32.1)
x13 = c(30.9,31.6,31.0,31.7,30.5,31.4,31.8,31.7,30.8,32.2)
x14 = c(30.9,31.7,31.3,31.3,30.8,31.2,31.7,31.5,30.6,32.4)

x21 = c(31.5,31.2,31.3,30.4,30.7,29.8,31.4,30.9,31.1,31.3)
x22 = c(30.6,31.2,31.3,30.8,30.9,30.8,32.0,32.4,31.3,31.5)
x23 = c(30.8,31.1,31.5,30.4,30.9,30.9,31.7,31.8,31.2,31.6)
x24 = c(31.0,31.3,31.4,30.2,30.9,30.8,31.6,31.9,31.2,31.7)

x1 = c(x11,x12,x13,x14)
x2 = c(x21,x22,x23,x24)
x = c(x1, x2)

# Compute the pooled variance matrix
X1=cbind(x11,x12,x13,x14)
X2=cbind(x21,x22,x23,x24)
S1 = var(X1)
S2= var(X2)
S = (S1+S2)/2

# Compute coorelation matrix. 
sd.inv = diag( sqrt(diag(V))^{-1} )
R = sd.inv%*%V%*%sd.inv

# Using lm to compute NNOVA table

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

group = factor(c(rep(1,40),rep(2,40))) 
time = factor(rep(c(rep(1,10), rep(2,10),rep(3,10),rep(4,10)),2))
lm.fit = lm(x~group+time+group*time)
anova(lm.fit)

time = factor(c(rep(1,10), rep(2,10),rep(3,10),rep(4,10)))
subject = factor(rep(c(1:10),4))
lm.fit1 = lm(x1~time+subject)
anova(lm.fit1)
lm.fit2 = lm(x2~time+subject)
anova(lm.fit2)

gtss = 1.2751+0.4104+0.3364 
ttss = gts+22.4070

sss = 10.0860+7.0902
rss = ttss - gtss-sss

# Compute test statistic for group by time interaction

X1.bar = apply(X1, 2, mean)
X2.bar = apply(X2, 2, mean)

c1 = c(-3,-1,1,3)
c2 = c(1,-1,-1,1)
c3 = c(-1,3,-3,1)

C = rbind(c1,c2,c3)

X.d = C%*%(X1.bar-X2.bar)

S.c = C%*%S%*%t(C)

T2.i = 10*10/(10+10) * t(X.d) %*% solve(S.c) %*% X.d 

F.i = (10+10-4)/(10+10-2)/(4-1)*T2.i

# Compute test statistic for overall time trends

Y = C%*%(X1.bar+X2.bar)/2 

T2.t = 20 * t(Y) %*% solve(S.c) %*% Y

F.t = (20-4)/(20-2)/(4-1)*T2.i

#  the three contrasts are the components of Y above.