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) X012 <- cbind(X0,XY1[,1],XY1[,2]) # Bind the columns for Intercept (X0), X1 (XY[,1]), X2 (XY[,2]) into X012 Y <- as.matrix(XY1[,3]) # Create Y from the third column of XY1 X01 <- X012[,1:2] X02 <- X012[,-2] P0 <- X0 %*% solve(t(X0) %*% X0) %*% t(X0) P01 <- X01 %*% solve(t(X01) %*% X01) %*% t(X01) P02 <- X02 %*% solve(t(X02) %*% X02) %*% t(X02) P012 <- X012 %*% solve(t(X012) %*% X012) %*% t(X012) betahat <- solve(t(X012) %*% X012) %*% t(X012) %*% Y s2 <- (t(Y) %*% (diag(nrow(Y))-P012) %*% Y)/(nrow(Y)-nrow(betahat)) V_betahat <- s2[1,1] * solve(t(X012) %*% X012) s_betahat <- sqrt(diag(V_betahat)) t_betahat <- betahat/s_betahat p_betahat <- 2*(1-pt(abs(t_betahat),nrow(Y)-nrow(betahat))) print(round(cbind(betahat,s_betahat,t_betahat,p_betahat),4)) R0 <- t(Y) %*% P0 %*% Y R01 <- t(Y) %*% P01 %*% Y R02 <- t(Y) %*% P02 %*% Y R012 <- t(Y) %*% P012 %*% Y R_1_0 <- R01-R0 R_2_01 <- R012-R01 R_2_0 <- R02-R0 R_1_02 <- R012-R02 print(round(c(R_1_0,R_2_01),3)) print(round(c(R_2_0,R_1_02),3)) dev.off()