> surgunit1 <- read.table("F:\\data\\CH09TA01.txt",header=F, + col.names=c("X1","X2","X3","X4","X5","X6","X7","X8","Y","lnY")) > > attach(surgunit1) > > > > su.mod0t<- lm(lnY ~ X1+X2+X3+X4+X5+X6+X7+X8) > (MSE.0t <- deviance(su.mod0t)/df.residual(su.mod0t)) [1] 0.04379427 > > > su.mod1t<- lm(lnY ~ X1+X2+X3+X8) > summary(su.mod1t) Call: lm(formula = lnY ~ X1 + X2 + X3 + X8) Residuals: Min 1Q Median 3Q Max -0.45307 -0.16149 -0.02779 0.12073 0.59524 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.852419 0.192695 19.992 < 2e-16 *** X1 0.073323 0.018973 3.865 0.000327 *** X2 0.014185 0.001731 8.196 9.58e-11 *** X3 0.015453 0.001396 11.072 6.15e-15 *** X8 0.352968 0.077191 4.573 3.29e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2109 on 49 degrees of freedom Multiple R-squared: 0.8299, Adjusted R-squared: 0.816 F-statistic: 59.76 on 4 and 49 DF, p-value: < 2.2e-16 > b.1t <- coef(su.mod1t) > X.1t <- model.matrix(su.mod1t) > (SSE.1t <- deviance(su.mod1t)) [1] 2.178799 > (MSE.1t <- SSE.1t/df.residual(su.mod1t)) [1] 0.04446529 > (Cp.1t <- (SSE.1t/MSE.0t) - (length(lnY)-2*ncol(X.1t))) [1] 5.750774 > H.1t <- diag(X.1t %*% solve(t(X.1t) %*% X.1t) %*% t(X.1t)) > d_res.1t <- residuals(su.mod1t)/(1-H.1t) > (PRESS.1t <- sum(d_res.1t^2)) [1] 2.737771 > (R2adj.1t <- 1-(MSE.1t/var(lnY))) [1] 0.815997 > > su.mod2t<- lm(lnY ~ X1+X2+X3+X6+X8) > summary(su.mod2t) Call: lm(formula = lnY ~ X1 + X2 + X3 + X6 + X8) Residuals: Min 1Q Median 3Q Max -0.40657 -0.13030 -0.02359 0.14138 0.54064 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.867095 0.190572 20.292 < 2e-16 *** X1 0.071241 0.018791 3.791 0.000419 *** X2 0.013890 0.001721 8.073 1.71e-10 *** X3 0.015115 0.001397 10.821 1.80e-14 *** X6 0.086910 0.058180 1.494 0.141768 X8 0.362677 0.076515 4.740 1.94e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2083 on 48 degrees of freedom Multiple R-squared: 0.8374, Adjusted R-squared: 0.8205 F-statistic: 49.46 on 5 and 48 DF, p-value: < 2.2e-16 > b.2t <- coef(su.mod2t) > X.2t <- model.matrix(su.mod2t) > (SSE.2t <- deviance(su.mod2t)) [1] 2.082008 > (MSE.2t <- SSE.2t/df.residual(su.mod2t)) [1] 0.04337516 > (Cp.2t <- (SSE.2t/MSE.0t) - (length(lnY)-2*ncol(X.2t))) [1] 5.540639 > H.2t <- diag(X.2t %*% solve(t(X.2t) %*% X.2t) %*% t(X.2t)) > d_res.1t <- residuals(su.mod2t)/(1-H.2t) > (PRESS.2t <- sum(d_res.2t^2)) Error: object 'd_res.2t' not found > (R2adj.2t <- 1-(MSE.2t/var(lnY))) [1] 0.8205081 > > su.mod3t<- lm(lnY ~ X1+X2+X3+X5+X6+X8) > summary(su.mod3t) Call: lm(formula = lnY ~ X1 + X2 + X3 + X5 + X6 + X8) Residuals: Min 1Q Median 3Q Max -0.34608 -0.13506 -0.02851 0.13239 0.48733 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.053974 0.234794 17.266 < 2e-16 *** X1 0.071517 0.018637 3.837 0.00037 *** X2 0.013755 0.001709 8.047 2.17e-10 *** X3 0.015116 0.001385 10.912 1.78e-14 *** X5 -0.003450 0.002572 -1.342 0.18620 X6 0.087317 0.057702 1.513 0.13691 X8 0.350904 0.076391 4.594 3.28e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2066 on 47 degrees of freedom Multiple R-squared: 0.8434, Adjusted R-squared: 0.8234 F-statistic: 42.2 on 6 and 47 DF, p-value: < 2.2e-16 > b.3t <- coef(su.mod3t) > X.3t <- model.matrix(su.mod3t) > (SSE.3t <- deviance(su.mod3t)) [1] 2.005225 > (MSE.3t <- SSE.3t/df.residual(su.mod3t)) [1] 0.04266437 > (Cp.3t <- (SSE.3t/MSE.0t) - (length(lnY)-2*ncol(X.3t))) [1] 5.787389 > H.3t <- diag(X.3t %*% solve(t(X.3t) %*% X.3t) %*% t(X.3t)) > d_res.3t <- residuals(su.mod3t)/(1-H.3t) > (PRESS.3t <- sum(d_res.3t^2)) [1] 2.772325 > (R2adj.3t <- 1-(MSE.3t/var(lnY))) [1] 0.8234494 > > detach(surgunit1) > > ######################### Read in Validation Sample > > surgunit2 <- read.table("F:\\data\\CH09TA05.txt",header=F, + col.names=c("X1","X2","X3","X4","X5","X6","X7","X8","Y","lnY")) > > attach(surgunit2) > > ################## Compute MSPR = mean squared prediction error for each model > > X.1v <- cbind(rep(1,length(lnY)),X1,X2,X3,X8) > yhat.1v <- X.1v %*% b.1t > (MSPR.1v <- sum((lnY-yhat.1v)^2)/length(lnY)) [1] 0.07731876 > > X.2v <- cbind(rep(1,length(lnY)),X1,X2,X3,X6,X8) > yhat.2v <- X.2v %*% b.2t > (MSPR.2v <- sum((lnY-yhat.2v)^2)/length(lnY)) [1] 0.0764202 > > X.3v <- cbind(rep(1,length(lnY)),X1,X2,X3,X5,X6,X8) > yhat.3v <- X.3v %*% b.3t > (MSPR.3v <- sum((lnY-yhat.3v)^2)/length(lnY)) [1] 0.07941142 > > #################### Re-Fit models 1-3 for Validation Sample > > > su.mod0v<- lm(lnY ~ X1+X2+X3+X4+X5+X6+X7+X8) > (MSE.0v <- deviance(su.mod0v)/df.residual(su.mod0v)) [1] 0.07558445 > > > su.mod1v<- lm(lnY ~ X1+X2+X3+X8) > summary(su.mod1v) Call: lm(formula = lnY ~ X1 + X2 + X3 + X8) Residuals: Min 1Q Median 3Q Max -0.61149 -0.13057 0.01547 0.17307 0.53690 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.635003 0.289364 12.562 < 2e-16 *** X1 0.095780 0.031939 2.999 0.00425 ** X2 0.016352 0.002322 7.044 5.68e-09 *** X3 0.015647 0.002027 7.720 5.14e-10 *** X8 0.185962 0.096435 1.928 0.05961 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2783 on 49 degrees of freedom Multiple R-squared: 0.7064, Adjusted R-squared: 0.6824 F-statistic: 29.47 on 4 and 49 DF, p-value: 1.672e-12 > b.1v <- coef(su.mod1v) > X.1v <- model.matrix(su.mod1v) > (SSE.1v <- deviance(su.mod1v)) [1] 3.795053 > (MSE.1v <- SSE.1v/df.residual(su.mod1v)) [1] 0.07745006 > (Cp.1v <- (SSE.1v/MSE.0v) - (length(lnY)-2*ncol(X.1v))) [1] 6.209438 > H.1v <- diag(X.1v %*% solve(t(X.1v) %*% X.1v) %*% t(X.1v)) > d_res.1v <- residuals(su.mod1v)/(1-H.1v) > (PRESS.1v <- sum(d_res.1v^2)) [1] 4.521911 > (R2adj.1v <- 1-(MSE.1v/var(lnY))) [1] 0.6824148 > > su.mod2v<- lm(lnY ~ X1+X2+X3+X6+X8) > summary(su.mod2v) Call: lm(formula = lnY ~ X1 + X2 + X3 + X6 + X8) Residuals: Min 1Q Median 3Q Max -0.57757 -0.13516 0.01021 0.13439 0.49539 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.614319 0.290665 12.435 < 2e-16 *** X1 0.099885 0.032295 3.093 0.0033 ** X2 0.015910 0.002374 6.702 2.1e-08 *** X3 0.015447 0.002042 7.566 1.0e-09 *** X6 0.073091 0.079154 0.923 0.3604 X8 0.188575 0.096622 1.952 0.0568 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2787 on 48 degrees of freedom Multiple R-squared: 0.7115, Adjusted R-squared: 0.6815 F-statistic: 23.68 on 5 and 48 DF, p-value: 6.482e-12 > b.2v <- coef(su.mod2v) > X.2v <- model.matrix(su.mod2v) > (SSE.2v <- deviance(su.mod2v)) [1] 3.728815 > (MSE.2v <- SSE.2v/df.residual(su.mod2v)) [1] 0.07768364 > (Cp.2v <- (SSE.2v/MSE.0v) - (length(lnY)-2*ncol(X.2v))) [1] 7.333092 > H.2v <- diag(X.2v %*% solve(t(X.2v) %*% X.2v) %*% t(X.2v)) > d_res.1v <- residuals(su.mod2v)/(1-H.2v) > (PRESS.2v <- sum(d_res.2v^2)) Error: object 'd_res.2v' not found > (R2adj.2v <- 1-(MSE.2v/var(lnY))) [1] 0.681457 > > su.mod3v<- lm(lnY ~ X1+X2+X3+X5+X6+X8) > summary(su.mod3v) Call: lm(formula = lnY ~ X1 + X2 + X3 + X5 + X6 + X8) Residuals: Min 1Q Median 3Q Max -0.60092 -0.12580 0.01123 0.13493 0.49742 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.469928 0.346781 10.006 3.14e-13 *** X1 0.098737 0.032466 3.041 0.00385 ** X2 0.016203 0.002414 6.712 2.23e-08 *** X3 0.015582 0.002058 7.572 1.12e-09 *** X5 0.002540 0.003293 0.771 0.44443 X6 0.072651 0.079493 0.914 0.36542 X8 0.193054 0.097206 1.986 0.05288 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2799 on 47 degrees of freedom Multiple R-squared: 0.7151, Adjusted R-squared: 0.6787 F-statistic: 19.66 on 6 and 47 DF, p-value: 2.525e-11 > b.3v <- coef(su.mod3v) > X.3v <- model.matrix(su.mod3v) > (SSE.3v <- deviance(su.mod3v)) [1] 3.682215 > (MSE.3v <- SSE.3v/df.residual(su.mod3v)) [1] 0.07834499 > (Cp.3v <- (SSE.3v/MSE.0v) - (length(lnY)-2*ncol(X.3v))) [1] 8.716564 > H.3v <- diag(X.3v %*% solve(t(X.3v) %*% X.3v) %*% t(X.3v)) > d_res.3v <- residuals(su.mod3v)/(1-H.3v) > (PRESS.3v <- sum(d_res.3v^2)) [1] 4.898106 > (R2adj.3v <- 1-(MSE.3v/var(lnY))) [1] 0.6787451 >