> > prob4_6 <- read.table("http://www.stat.ufl.edu/~winner/sta6208/rpd/PEAKFLOW_DTA.txt",header=F, + col.names=c("Qo","Qp","X1","X2","X3","X4")) > > attach(prob4_6) > > Y <- log(Qo/Qp) # Create Y > > (SSTOU <- sum(Y^2)) # Total (Uncorrected) SS [1] 0.7742164 > > mod1 <- lm(Y ~ X1 + X4 + X2 + X3) # Full model in order given in problem (includes intercept) > summary(mod1) # Regression coefficients and partial t-tests and Global F-test Call: lm(formula = Y ~ X1 + X4 + X2 + X3) Residuals: 1 2 3 4 5 6 7 8 -0.02800 0.03143 0.01892 0.01694 -0.01530 -0.02860 -0.02467 0.02501 9 10 -0.04129 0.04556 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.087321 0.311931 0.280 0.79074 X1 -0.035384 0.009336 -3.790 0.01276 * X4 -0.120080 0.023810 -5.043 0.00396 ** X2 0.004726 0.004765 0.992 0.36677 X3 -0.001913 0.004189 -0.457 0.66702 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04119 on 5 degrees of freedom Multiple R-squared: 0.9287, Adjusted R-squared: 0.8717 F-statistic: 16.29 on 4 and 5 DF, p-value: 0.004502 > anova(mod1) # Sequential SS: R(B1|B0), R(B4|B0,B1), R(B2|B0,B1,B4), R(B3|B0,B1,B4,B2) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 0.055536 0.055536 32.7313 0.002282 ** X4 1 0.048498 0.048498 28.5834 0.003073 ** X2 1 0.006189 0.006189 3.6477 0.114403 X3 1 0.000354 0.000354 0.2086 0.667019 Residuals 5 0.008484 0.001697 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > drop1(mod1) # Partial SS: R(B1|B0,B2,B3,B4), R(B4|B0,B1,B2,B3), R(B2|B0,B1,B3,B4), R(B3|B0,B1,B4,B2) Single term deletions Model: Y ~ X1 + X4 + X2 + X3 Df Sum of Sq RSS AIC 0.008484 -60.722 X1 1 0.024370 0.032854 -49.183 X4 1 0.043156 0.051640 -44.660 X2 1 0.001669 0.010153 -60.926 X3 1 0.000354 0.008838 -62.313 > > (SSRES1 <- deviance(mod1)) # SSRes(Mod1) where deviance(mod1)=SSRes(Mod1) [1] 0.008483568 > (dfRES1 <- df.residual(mod1)) # dfRes(Mod1) [1] 5 > > (SSMODEL <- SSTOU - SSRES1) # SSModel(Mod1) [1] 0.7657328 > > (SSREG <- SSMODEL - 10*(mean(Y)^2)) # SSReg(Mod1) [1] 0.1105765 > > > > mod2 <- lm(Y ~ X1 + X4 + X2 + X3 -1) # Model with all 4 predictors, no intercept > summary(mod2) Call: lm(formula = Y ~ X1 + X4 + X2 + X3 - 1) Residuals: Min 1Q Median 3Q Max -0.04001 -0.02546 -0.00017 0.02573 0.04582 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 -0.0330400 0.0037987 -8.698 0.000128 *** X4 -0.1215292 0.0213807 -5.684 0.001279 ** X2 0.0057518 0.0028025 2.052 0.085956 . X3 -0.0007525 0.0005484 -1.372 0.219124 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.0379 on 6 degrees of freedom Multiple R-squared: 0.9889, Adjusted R-squared: 0.9815 F-statistic: 133.3 on 4 and 6 DF, p-value: 5.468e-06 > anova(mod2) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 0.56463 0.56463 393.1687 1.067e-06 *** X4 1 0.19451 0.19451 135.4416 2.424e-05 *** X2 1 0.00377 0.00377 2.6217 0.1565 X3 1 0.00270 0.00270 1.8826 0.2191 Residuals 6 0.00862 0.00144 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > drop1(mod2) Single term deletions Model: Y ~ X1 + X4 + X2 + X3 - 1 Df Sum of Sq RSS AIC 0.008617 -62.567 X1 1 0.108640 0.117256 -38.460 X4 1 0.046398 0.055014 -46.027 X2 1 0.006049 0.014666 -59.248 X3 1 0.002704 0.011320 -61.838 > > (SSRES2 <- deviance(mod2)) [1] 0.00861653 > (dfRES2 <- df.residual(mod2)) [1] 6 > > (F_B0 <- ((SSRES2-SSRES1)/(dfRES2-dfRES1))/(SSRES1/dfRES1)) [1] 0.07836508 > > (F_05 <- qf(.95,dfRES2-dfRES1,dfRES1)) [1] 6.607891 > > anova(mod2,mod1) Analysis of Variance Table Model 1: Y ~ X1 + X4 + X2 + X3 - 1 Model 2: Y ~ X1 + X4 + X2 + X3 Res.Df RSS Df Sum of Sq F Pr(>F) 1 6 0.0086165 2 5 0.0084836 1 0.00013296 0.0784 0.7907 > > mod3 <- lm(Y ~ X1 + X4 -1) # Model with X1, X4, no intercept > summary(mod3) Call: lm(formula = Y ~ X1 + X4 - 1) Residuals: Min 1Q Median 3Q Max -0.057972 -0.037806 0.004037 0.025781 0.059843 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 -0.033182 0.003892 -8.525 2.76e-05 *** X4 -0.124274 0.012236 -10.156 7.56e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04342 on 8 degrees of freedom Multiple R-squared: 0.9805, Adjusted R-squared: 0.9756 F-statistic: 201.3 on 2 and 8 DF, p-value: 1.441e-07 > anova(mod3) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 0.56463 0.56463 299.43 1.267e-07 *** X4 1 0.19451 0.19451 103.15 7.559e-06 *** Residuals 8 0.01509 0.00189 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > drop1(mod3) Single term deletions Model: Y ~ X1 + X4 - 1 Df Sum of Sq RSS AIC 0.015085 -60.966 X1 1 0.13703 0.152119 -39.857 X4 1 0.19451 0.209591 -36.652 > > (SSRES3 <- deviance(mod3)) [1] 0.01508514 > (dfRES3 <- df.residual(mod3)) [1] 8 > > (F_B2B3 <- ((SSRES3-SSRES2)/(dfRES3-dfRES2))/(SSRES2/dfRES2)) [1] 2.252164 > > (F_05 <- qf(.95,dfRES3-dfRES2,dfRES2)) [1] 5.143253 > > anova(mod3,mod2) Analysis of Variance Table Model 1: Y ~ X1 + X4 - 1 Model 2: Y ~ X1 + X4 + X2 + X3 - 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 8 0.0150851 2 6 0.0086165 2 0.0064686 2.2522 0.1864 > > > ###### Matrix Approach for F-tests > n <- length(Y) > > (Y <- matrix(Y,ncol=1)) ## Y-vector [,1] [1,] -0.13353139 [2,] -0.23732819 [3,] -0.23214811 [4,] -0.02304840 [5,] -0.29300788 [6,] -0.26261197 [7,] -0.31260102 [8,] -0.34386172 [9,] -0.44618181 [10,] -0.27528156 > > X0 <- rep(1,length(Y)) ### Intercept Column > > (XM1 <- matrix(cbind(X0,X1,X2,X3,X4),ncol=5)) # X-matrix for Model 1 [,1] [,2] [,3] [,4] [,5] [1,] 1 0.03 3.0 70 0.6 [2,] 1 0.03 3.0 80 1.8 [3,] 1 0.13 6.5 65 2.0 [4,] 1 1.00 15.0 60 0.4 [5,] 1 1.00 15.0 65 2.3 [6,] 1 3.00 7.0 67 1.0 [7,] 1 5.00 6.0 62 0.9 [8,] 1 7.00 6.5 56 1.1 [9,] 1 7.00 6.5 56 1.4 [10,] 1 7.00 6.5 56 0.7 > (XMXMI1 <- solve(t(XM1) %*% XM1)) [,1] [,2] [,3] [,4] [,5] [1,] 57.3468515 -1.53947322 -0.67353914 -0.76224179 0.95173019 [2,] -1.5394732 0.05137534 0.01736300 0.01995079 -0.01791338 [3,] -0.6735391 0.01736300 0.01337964 0.00849132 -0.01791842 [4,] -0.7622418 0.01995079 0.00849132 0.01034099 -0.01835306 [5,] 0.9517302 -0.01791338 -0.01791842 -0.01835306 0.33411476 > > (XM2 <- matrix(cbind(X1,X2,X3,X4),ncol=4)) # X-matrix for Model 2 [,1] [,2] [,3] [,4] [1,] 0.03 3.0 70 0.6 [2,] 0.03 3.0 80 1.8 [3,] 0.13 6.5 65 2.0 [4,] 1.00 15.0 60 0.4 [5,] 1.00 15.0 65 2.3 [6,] 3.00 7.0 67 1.0 [7,] 5.00 6.0 62 0.9 [8,] 7.00 6.5 56 1.1 [9,] 7.00 6.5 56 1.4 [10,] 7.00 6.5 56 0.7 > (XMXMI2 <- solve(t(XM2) %*% XM2)) [,1] [,2] [,3] [,4] [1,] 0.0100482664 -0.0007181195 -0.0005115508 0.007635763 [2,] -0.0007181195 0.0054689183 -0.0004612143 -0.006740347 [3,] -0.0005115508 -0.0004612143 0.0002094399 -0.005702875 [4,] 0.0076357629 -0.0067403474 -0.0057028747 0.318319816 > > (betahat1 <- XMXMI1 %*% t(XM1) %*% Y) [,1] [1,] 0.087321331 [2,] -0.035384099 [3,] 0.004726196 [4,] -0.001913147 [5,] -0.120079975 > (betahat2 <- XMXMI2 %*% t(XM2) %*% Y) [,1] [1,] -0.0330399627 [2,] 0.0057517861 [3,] -0.0007524905 [4,] -0.1215291626 > > > (K1P <- matrix(c(1,0,0,0,0),ncol=5)) [,1] [,2] [,3] [,4] [,5] [1,] 1 0 0 0 0 > > (K2P <- matrix(c(0,1,0,0, 0,0,1,0),byrow=T,ncol=4)) [,1] [,2] [,3] [,4] [1,] 0 1 0 0 [2,] 0 0 1 0 > > k1 <- nrow(K1P) > k2 <- nrow(K2P) > > > ######### SSRes = Y'(I-P)Y s^2=MSRes > > (SSRes1 <- t(Y) %*% (diag(n) - (XM1 %*% XMXMI1 %*% t(XM1))) %*% Y) [,1] [1,] 0.008483568 > (s2_1 <- SSRes1/(n-ncol(XM1))) [,1] [1,] 0.001696714 > > (SSRes2 <- t(Y) %*% (diag(n) - (XM2 %*% XMXMI2 %*% t(XM2))) %*% Y) [,1] [1,] 0.00861653 > (s2_2 <- SSRes2/(n-ncol(XM2))) [,1] [1,] 0.001436088 > > > (F1 <- + (t(K1P %*% betahat1) %*% solve(K1P %*% XMXMI1 %*% t(K1P)) %*% (K1P %*% betahat1))/(k1*s2_1)) [,1] [1,] 0.07836508 > > (F2 <- + (t(K2P %*% betahat2) %*% solve(K2P %*% XMXMI2 %*% t(K2P)) %*% (K2P %*% betahat2))/(k2*s2_2)) [,1] [1,] 2.252164 > > >