pdf("muscle1.pdf") XY1 <- matrix(c(43.7,19,177, 43.7,43,279, 43.7,56,346, 54.6,13,160, 54.6,19,193, 54.6,43,280, 54.6,56,335, 55.7,13,169, 55.7,26,212, 55.7,34.5,244, 55.7,43,285, 58.8,13,181, 58.8,43,298, 60.5,19,212, 60.5,43,317, 60.5,56,347, 61.9,13,186, 61.9,19,216, 61.9,34.5,265, 61.9,43,306, 61.9,56,348, 66.7,13,209, 66.7,43,324, 66.7,56,352), byrow=T,ncol=3) # XY1 is a matrix containing n=24 rows and 3 columns (X1,X2,Y) we must form the X and Y matrices from it # This is similar to when we input a raw data frame (dataset) and create a matrix X0 <- matrix(rep(1,nrow(XY1)),ncol=1) # Create a Column of 1s for the intercept (same number of rows as XY1) X <- cbind(X0,XY1[,1],XY1[,2]) # Bind the columns for Intercept (X0), X1 (XY[,1]), X2 (XY[,2]) into X Y <- as.matrix(XY1[,3]) # Create Y from the third column of XY1 XXinv <- solve(t(X) %*% X) # Compute X'X-Inverse = ((X'X)^(-1)) betahat <- as.matrix(XXinv %*% t(X) %*% Y) # Compute betahat = ((X'X)^(-1))X'Y P <- X %*% XXinv %*% t(X) # Compute P = X((X'X)^(-1))X IMP <- diag(nrow(X))-P # Compute IMP = I-P Yhat <- as.matrix(P %*% Y) # Compute Yhat = PY e <- as.matrix(IMP %*% Y) # Compute e = (I-P)Y s2 <- (t(e) %*% e)/(nrow(Y)-nrow(betahat)) # s2 = (e'e)/(n-(p+1)) V_betahat <- s2[1,1]* XXinv # V(betahat) = s2((X'X)^(-1)) Note how s2 is treated as indexed scalar print(cbind(X,Y,Yhat,e)) XXinv betahat s2 V_betahat dev.off()