sleepdep <- read.table("http://www.stat.ufl.edu/~winner/data/caffeine_mixed.dat", header=F,col.names=c("subj_id","caff.trt","day0.x","day1.x","day2.x","day3.x")) attach(sleepdep) X.c <- cbind(day0.x[caff.trt==1],day1.x[caff.trt==1],day2.x[caff.trt==1],day3.x[caff.trt==1]) X.p <- cbind(day0.x[caff.trt==2],day1.x[caff.trt==2],day2.x[caff.trt==2],day3.x[caff.trt==2]) n.c <- nrow(X.c); n.p <- nrow(X.p); n <- n.c+n.p; g <- 2; p <- ncol(X.c) I.c <- diag(n.c); J.c <- matrix(rep(1,n.c^2),n.c,n.c) I.p <- diag(n.p); J.p <- matrix(rep(1,n.p^2),n.p,n.p) (xbar.c <- (1/n.c) * (t(X.c) %*% rep(1,n.c))) (xbar.p <- (1/n.p) * (t(X.p) %*% rep(1,n.p))) (S.c <- (1/(n.c-1)) * (t(X.c) %*% (I.c - (1/n.c)*J.c) %*% X.c)) (S.p <- (1/(n.p-1)) * (t(X.p) %*% (I.p - (1/n.p)*J.p) %*% X.p)) (Spooled <- ((n.c-1)*S.c + (n.p-1)*S.p) / (n.c+n.p-2)) B.l <- matrix(c(1,0,1,1,1,2,1,3),byrow=T,ncol=2) B.q <- matrix(c(1,0,0,1,1,1,1,2,4,1,3,9),byrow=T,ncol=3) B.c <- matrix(c(1,0,0,0,1,1,1,1,1,2,4,8,1,3,9,27),byrow=T,ncol=4) (beta.l.c <- solve(t(B.l) %*% solve(Spooled) %*% B.l) %*% t(B.l) %*% solve(Spooled) %*% xbar.c) (beta.l.p <- solve(t(B.l) %*% solve(Spooled) %*% B.l) %*% t(B.l) %*% solve(Spooled) %*% xbar.p) (beta.q.c <- solve(t(B.q) %*% solve(Spooled) %*% B.q) %*% t(B.q) %*% solve(Spooled) %*% xbar.c) (beta.q.p <- solve(t(B.q) %*% solve(Spooled) %*% B.q) %*% t(B.q) %*% solve(Spooled) %*% xbar.p) k.l <- ((n-g)*(n-g-1)) / ((n-g-p+1)*(n-g-p+1+1)) (cov.beta.l.c <- (k.l/n.c) * solve(t(B.l) %*% solve(Spooled) %*% B.l)) (cov.beta.l.p <- (k.l/n.p) * solve(t(B.l) %*% solve(Spooled) %*% B.l)) se.beta.l.c <- sqrt(diag(cov.beta.l.c)); t.beta.l.c <- beta.l.c/se.beta.l.c se.beta.l.p <- sqrt(diag(cov.beta.l.p)); t.beta.l.p <- beta.l.p/se.beta.l.p beta.l.out <- cbind(beta.l.c,se.beta.l.c,t.beta.l.c,beta.l.p,se.beta.l.p,t.beta.l.p) colnames(beta.l.out) <- c("beta.l.c","se.l.c","t.l.c", "beta.l.p","se.l.p","t.l.p") round(beta.l.out,3) k.q <- ((n-g)*(n-g-1)) / ((n-g-p+2)*(n-g-p+2+1)) (cov.beta.q.c <- (k.q/n.c) * solve(t(B.q) %*% solve(Spooled) %*% B.q)) (cov.beta.q.p <- (k.q/n.p) * solve(t(B.q) %*% solve(Spooled) %*% B.q)) se.beta.q.c <- sqrt(diag(cov.beta.q.c)); t.beta.q.c <- beta.q.c/se.beta.q.c se.beta.q.p <- sqrt(diag(cov.beta.q.p)); t.beta.q.p <- beta.q.p/se.beta.q.p beta.q.out <- cbind(beta.q.c,se.beta.q.c,t.beta.q.c,beta.q.p,se.beta.q.p,t.beta.q.p) colnames(beta.q.out) <- c("beta.q.c","se.q.c","t.q.c", "beta.q.p","se.q.p","t.q.p") round(beta.q.out,3) (Xhat.l.c <- t(B.l %*% beta.l.c)) (Xhat.l.p <- t(B.l %*% beta.l.p)) Xhatm.l.c <- matrix(rep(Xhat.l.c,n.c),byrow=T,ncol=p) Xhatm.l.p <- matrix(rep(Xhat.l.p,n.p),byrow=T,ncol=p) Res.l.c <- X.c - Xhatm.l.c Res.l.p <- X.p - Xhatm.l.p W.l <- (t(Res.l.c) %*% Res.l.c) + (t(Res.l.p) %*% Res.l.p) W <- (n.c+n.p-2) * Spooled Lambda.star.l <- det(W) / det(W.l) GF.TS.l <- -(n-(p-1+g)/2) * log(Lambda.star.l) GF.df.l <- (p-1-1)*g GF.CV.l <- qchisq(.95,GF.df.l) GF.PV.l <- 1-pchisq(GF.TS.l,GF.df.l) L.l.out <- cbind(Lambda.star.l, GF.TS.l, GF.df.l, GF.CV.l, GF.PV.l) colnames(L.l.out) <- c("Lambda", "X2 stat", "df", "X2(.05)", "P-value") round(L.l.out,3) (Xhat.q.c <- t(B.q %*% beta.q.c)) (Xhat.q.p <- t(B.q %*% beta.q.p)) Xhatm.q.c <- matrix(rep(Xhat.q.c,n.c),byrow=T,ncol=p) Xhatm.q.p <- matrix(rep(Xhat.q.p,n.p),byrow=T,ncol=p) Res.q.c <- X.c - Xhatm.q.c Res.q.p <- X.p - Xhatm.q.p W.q <- (t(Res.q.c) %*% Res.q.c) + (t(Res.q.p) %*% Res.q.p) W <- (n.c+n.p-2) * Spooled Lambda.star.q <- det(W) / det(W.q) GF.TS.q <- -(n-(p-2+g)/2) * log(Lambda.star.q) GF.df.q <- (p-2-1)*g GF.CV.q <- qchisq(.95,GF.df.q) GF.PV.q <- 1-pchisq(GF.TS.q,GF.df.q) L.q.out <- cbind(Lambda.star.q, GF.TS.q, GF.df.q, GF.CV.q, GF.PV.q) colnames(L.q.out) <- c("Lambda", "X2 stat", "df", "X2(.05)", "P-value") round(L.q.out,3)