[Previously saved workspace restored] > dwaine <- read.table("F:\\sta4210\\dwaine_studios.txt",header=F,col.names=c("Y","X1","X2")) > > attach(dwaine) > > dwaine Y X1 X2 1 174.4 68.5 16.7 2 164.4 45.2 16.8 3 244.2 91.3 18.2 4 154.6 47.8 16.3 5 181.6 46.9 17.3 6 207.5 66.1 18.2 7 152.8 49.5 15.9 8 163.2 52.0 17.2 9 145.4 48.9 16.6 10 137.2 38.4 16.0 11 241.9 87.9 18.3 12 191.1 72.8 17.1 13 232.0 88.4 17.4 14 145.3 42.9 15.8 15 161.1 52.5 17.8 16 209.7 85.7 18.4 17 146.4 41.3 16.5 18 144.0 51.7 16.3 19 232.6 89.6 18.1 20 224.1 82.7 19.1 21 166.5 52.3 16.0 > > plot(dwaine) > > (n <- length(Y)) [1] 21 > > X0 <- rep(1,n) > > X <- as.matrix(cbind(X0,X1,X2)) > Y <- as.matrix(Y) > > (XX <- t(X) %*% X) # Compute X'X (t(X) is transpose of X, %*% performs matrix multiplication) X0 X1 X2 X0 21.0 1302.40 360.00 X1 1302.4 87707.94 22609.19 X2 360.0 22609.19 6190.26 > (XY <- t(X) %*% Y) # Compute X'Y [,1] X0 3820.00 X1 249643.35 X2 66072.75 > (XXI <- solve(XX)) # Compute (X'X)^(-1) (solve obtains inverse of square, full rank matrix) X0 X1 X2 X0 29.72892348 0.0721834719 -1.992553186 X1 0.07218347 0.0003701761 -0.005549917 X2 -1.99255319 -0.0055499173 0.136310637 > (b <- XXI %*% XY) # Compute b = (X'X)^(-1)X'Y [,1] X0 -68.85707 X1 1.45456 X2 9.36550 > > J_n <- matrix(rep(1/n,n^2),byrow=T,ncol=n) # Obtain The (1/n)J matrix > I_n <- diag(n) # Obtain the identity matrix > H <- X %*% XXI %*% t(X) # Compute the hat matrix > > (SSTO <- t(Y) %*% (I_n - J_n) %*% Y) # Compute the Total (Corrected) Sum of Squares [,1] [1,] 26196.21 > (SSE <- t(Y) %*% (I_n - H) %*% Y) # Compute the Error Sum of Squares [,1] [1,] 2180.927 > (SSR <- t(Y) %*% (H - J_n) %*% Y) # Compute the Regression Sum of Squares [,1] [1,] 24015.28 > > (df_TO <- n-1) # Obtain the Total Degrees of Freedom [1] 20 > (df_E <- n-ncol(X)) # Obtain the Error Degrees of Freedom [1] 18 > (df_R <- ncol(X)-1) # Obtain the Regression Degrees of Freedom [1] 2 > > (MSE <- SSE/df_E) # Compute Mean square Error [,1] [1,] 121.1626 > (MSR <- SSR/df_R) # Compute Mean Square Regression [,1] [1,] 12007.64 > > (F_obs <- MSR/MSE) # Compute F* for testing H0: beta1=beta2=0 [,1] [1,] 99.1035 > (F_05 <- qf(.95,df_E,df_R)) # Obtain critical F-value (alpha=0.05) [1] 19.44022 > (P_F_obs <- 1-pf(F_obs,df_E,df_R)) # Obtain the P-value for F-test [,1] [1,] 0.01003413 > > (R2 <- SSR/SSTO) # Compute the Coefficient of Multiple Determination [,1] [1,] 0.9167465 > > (s2_b <- MSE[1,1]*XXI) # Compute the Variance-Covariance Matrix of Regression Coefficients X0 X1 X2 X0 3602.03467 8.74593958 -241.4229923 X1 8.74594 0.04485151 -0.6724426 X2 -241.42299 -0.67244260 16.5157558 > > b1 <- b[2] # Retrieve b1 from Vector b > b2 <- b[3] # Retrieve b2 from Vector b > > (b1_95CI <- c(b1-qt(.975,df_E)*sqrt(s2_b[2,2]),b1+qt(.975,df_E)*sqrt(s2_b[2,2]))) # 95% CI for beta1 [1] 1.009623 1.899497 > (b2_95CI <- c(b2-qt(.975,df_E)*sqrt(s2_b[3,3]),b2+qt(.975,df_E)*sqrt(s2_b[3,3]))) # 95% CI for beta2 [1] 0.8274411 17.9035596 > > xh1 <- as.matrix(c(1,65.4,17.6),ncol=1) # Set up X_h vector when X1=65.4, X2=17.6 > (yh1 <- t(xh1) %*% b) # Compute Y-hat when X1=65.4, X2=17.6 [,1] [1,] 191.1039 > (s2_yh1 <- t(xh1) %*% s2_b %*% xh1) # Compute s2{Y-hat} when X1=65.4, X2=17.6 [,1] [1,] 7.65517 > (s_yh1 <- sqrt(s2_yh1)) # Compute s{Y-hat} when X1=65.4, X2=17.6 [,1] [1,] 2.766798 > > (yh1_95CI <- c(yh1-qt(.975,df_E)*s_yh1,yh1+qt(.975,df_E)*s_yh1)) # 95% CI for mean when X1=65.4, X2=16 [1] 185.2911 196.9168 > > xhpredA <- as.matrix(c(1,65.4,17.6),ncol=1) # Set up X_h vector when X1=65.4, X2=17.6 for City A > xhpredB <- as.matrix(c(1,53.1,17.7),ncol=1) # Set up X_h vector when X1=53.1, X2=17.7 for City B > (yhpredA <- t(xhpredA) %*% b) # Obtain predicted value for City A [,1] [1,] 191.1039 > (yhpredB <- t(xhpredB) %*% b) # Obtain predicted value for City B [,1] [1,] 174.1494 > (s2_yhpredA <- MSE + (t(xhpredA) %*% s2_b %*% xhpredA)) # Compute s2{yhpredA} [,1] [1,] 128.8178 > (s2_yhpredB <- MSE + (t(xhpredB) %*% s2_b %*% xhpredB)) # Compute s2{yhpredB} [,1] [1,] 142.3098 > (s_yhpredA <- sqrt(s2_yhpredA)) # Compute s{yhpredA} [,1] [1,] 11.34979 > (s_yhpredB <- sqrt(s2_yhpredB)) # Compute s{yhpredB} [,1] [1,] 11.92937 > > (Bon_90PI_A <- c(yhpredA-qt(1-.10/4,df_E)*s_yhpredA,yhpredA+qt(1-.10/4,df_E)*s_yhpredA)) # Bonferroni PI for City A (Joint 90% PI for both cities) [1] 167.2589 214.9490 > (Bon_90PI_B <- c(yhpredB-qt(1-.10/4,df_E)*s_yhpredB,yhpredB+qt(1-.10/4,df_E)*s_yhpredB)) # Bonferroni PI for City B (Joint 90% PI for both cities) [1] 149.0867 199.2121 > > (Scheffe_90_PI_A <- c(yhpredA-sqrt(2*qf(.90,2,df_E))*s_yhpredA,yhpredA+sqrt(2*qf(.90,2,df_E))*s_yhpredA)) # Scheffe PI for City A (Joint 90% PI for both cities) [1] 165.1035 217.1044 > (Scheffe_90_PI_B <- c(yhpredB-sqrt(2*qf(.90,2,df_E))*s_yhpredB,yhpredB+sqrt(2*qf(.90,2,df_E))*s_yhpredB)) # Scheffe PI for City B (Joint 90% PI for both cities) [1] 146.8213 201.4775 > > ############### Regression with lm Function > > dwaine.reg1 <- lm(Y~X1+X2) > summary(dwaine.reg1) Call: lm(formula = Y ~ X1 + X2) Residuals: Min 1Q Median 3Q Max -18.4239 -6.2161 0.7449 9.4356 20.2151 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -68.8571 60.0170 -1.147 0.2663 X1 1.4546 0.2118 6.868 2e-06 *** X2 9.3655 4.0640 2.305 0.0333 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 11.01 on 18 degrees of freedom Multiple R-squared: 0.9167, Adjusted R-squared: 0.9075 F-statistic: 99.1 on 2 and 18 DF, p-value: 1.921e-10 > anova(dwaine.reg1) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 23371.8 23371.8 192.8962 4.64e-11 *** X2 1 643.5 643.5 5.3108 0.03332 * Residuals 18 2180.9 121.2 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > drop1(dwaine.reg1) Single term deletions Model: Y ~ X1 + X2 Df Sum of Sq RSS AIC 2180.9 103.50 X1 1 5715.5 7896.4 128.52 X2 1 643.5 2824.4 106.93 > >