questions from Wichern - Multivariate Stat analysis

wichern6.34.R

data1 <-read.table("C:/Users/DASA0/Desktop/Stat 524/wichern data/T6-5.dat", sep="")
data2 <-read.table("C:/Users/DASA0/Desktop/Stat 524/wichern data/T6-6.dat", sep="")
data1 <- as.data.frame(data1) 
data2 <- as.data.frame(data2)

xbar1 <- apply(data1, 2, mean)
xbar2 <- apply(data2, 2, mean)
t <- c(0, 1, 2, 3)
df <- data.frame(x=t, y=xbar1, type='Profile 1') 
df <- rbind(df, data.frame(x=t, y=xbar2, type='Profile 2'))
plot(y~x, pch=16, xlim=c(0,3.5), ylim=c(63, 74), xlab = 'Time, in years', ylab = 'Profiles: xbar1, xbar2', main = 'Profiles for xbar1, xbar2', data = df)
with(df, text(y~x, labels = row.names(df), pos = 4))

S1 <- cov(data1)
S2 <- cov(data2)
n1 <- dim(data1)[1]
n2 <- dim(data2)[1]
Spooled <- ((n1-1)*S1 + (n2-1)*S2 )/(n1+n2-2)
Spoolinv <- solve(Spooled)
B <- matrix(c(1, 1, 1, 1, 0, 1, 2, 3), nc=2)
temp <- t(B)%*%Spoolinv%*%B
Xlbar <- matrix(cbind(xbar1, xbar2), nc=2)
MLE <- solve(temp)%*%t(B)%*%Spoolinv%*%Xlbar
N <- n1 +n2
g <- dim(B)[2]
W <- (N-g)*Spooled
p <- dim(data1)[2]
q <- 1 # linear growth model
k <- ((N-g)*(N-g-1))/((N-g-p+q)*(N-g-p+q+1))
CovBetaL1 <- (k/n1)*solve(temp)
CovBetaL2 <- (k/n2)*solve(temp) 

beta1<-  solve(t(B)%*%Spoolinv%*%B)%*%t(B)%*%Spoolinv%*%xbar1
beta2<-  solve(t(B)%*%Spoolinv%*%B)%*%t(B)%*%Spoolinv%*%xbar2
pred1 <- B%*%beta1
onevec1 <- rep(1,n1)
pred2 <- B%*%beta2
onevec2 <- rep(1,n2)
diff1 <- as.matrix(data1 - outer(onevec1, pred1))
diff2 <- as.matrix( data2 - outer(onevec2, pred2))
Wq <-  t(diff1)%*%diff1 + t(diff2)%*%diff2

WilkSLambda <- det(W)/det(Wq)

teststat <- -(N -0.5*(p-q+g))*log(WilkSLambda)
alpha <- 0.01
critical <- qchisq(1-alpha, (p-q-1)*g)