data <- read.table("data.dat",header=F)
attach(data)
Y <- as.matrix(data[,1])
X <- as.matrix(cbind(1, data[,2:ncol(data)]))
XTX <- t(X) %*% X
XTXINV <- solve(XTX)
XTY <- t(X) %*% Y
betahat <- XTXINV %*% XTY
Xbetahat <- X %*% betahat
sigsq <- as.vector(t(Y - Xbetahat) %*% (Y - Xbetahat) / (nrow(X) - ncol(X)))
VCOV <- sigsq * XTXINV
J <- 2
R <- matrix(0:0, J, ncol(X))
r <- matrix(0:0,J,1)
R[1,2] <- 1
R[2,3] <- 1
Rbetahat <- R %*% betahat
F <- t(Rbetahat - r) %*% solve(R %*% XTXINV %*% t(R)) %*% (Rbetahat - r) / (J * sigsq)
F
# Compare with the lm() command (linear regression model)
fitlm <- lm(Y ~ X[,2:ncol(X)])
fitreslm <- lm(Y ~ X[,4:ncol(X)])
anova(fitlm, fitreslm)