> cracker.mod1 <- lm(y~xc) > summary(cracker.mod1) Call: lm(formula = y ~ xc) Residuals: Min 1Q Median 3Q Max -8.711 -5.481 1.289 3.975 9.017 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 33.8000 1.5287 22.110 1.07e-11 *** xc 0.7278 0.3121 2.332 0.0364 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 5.921 on 13 degrees of freedom Multiple R-squared: 0.295, Adjusted R-squared: 0.2408 F-statistic: 5.439 on 1 and 13 DF, p-value: 0.03641 > anova(cracker.mod1) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) xc 1 190.68 190.678 5.4393 0.03641 * Residuals 13 455.72 35.056 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > cracker.mod2 <- lm(y~xc + i1 + i2) > summary(cracker.mod2) Call: lm(formula = y ~ xc + i1 + i2) Residuals: Min 1Q Median 3Q Max -2.4348 -1.2739 -0.3362 1.6710 2.4869 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 33.8000 0.4835 69.908 6.37e-16 *** xc 0.8986 0.1026 8.759 2.73e-06 *** i1 6.0174 0.7083 8.496 3.67e-06 *** i2 0.9420 0.6987 1.348 0.205 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.873 on 11 degrees of freedom Multiple R-squared: 0.9403, Adjusted R-squared: 0.9241 F-statistic: 57.78 on 3 and 11 DF, p-value: 5.082e-07 > anova(cracker.mod2) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) xc 1 190.68 190.68 54.3786 1.405e-05 *** i1 1 410.78 410.78 117.1478 3.335e-07 *** i2 1 6.37 6.37 1.8178 0.2047 Residuals 11 38.57 3.51 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > cracker.mod3 <- lm(y~xc + i1 + i2 + i1xc + i2xc) > summary(cracker.mod3) Call: lm(formula = y ~ xc + i1 + i2 + i1xc + i2xc) Residuals: Min 1Q Median 3Q Max -2.2555 -0.7820 -0.5773 1.1295 2.3965 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 33.89433 0.51234 66.155 2.08e-13 *** xc 0.93874 0.11267 8.332 1.60e-05 *** i1 6.26990 0.75167 8.341 1.58e-05 *** i2 0.71791 0.71600 1.003 0.342 i1xc 0.15251 0.18438 0.827 0.430 i2xc 0.05252 0.14561 0.361 0.727 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.871 on 9 degrees of freedom Multiple R-squared: 0.9512, Adjusted R-squared: 0.9241 F-statistic: 35.11 on 5 and 9 DF, p-value: 1.216e-05 > anova(cracker.mod3) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) xc 1 190.68 190.68 54.4434 4.198e-05 *** i1 1 410.78 410.78 117.2872 1.836e-06 *** i2 1 6.37 6.37 1.8200 0.2103 i1xc 1 6.59 6.59 1.8830 0.2032 i2xc 1 0.46 0.46 0.1301 0.7267 Residuals 9 31.52 3.50 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > anova(cracker.mod1,cracker.mod2) Analysis of Variance Table Model 1: y ~ xc Model 2: y ~ xc + i1 + i2 Res.Df RSS Df Sum of Sq F Pr(>F) 1 13 455.72 2 11 38.57 2 417.15 59.483 1.264e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(cracker.mod2,cracker.mod3) Analysis of Variance Table Model 1: y ~ xc + i1 + i2 Model 2: y ~ xc + i1 + i2 + i1xc + i2xc Res.Df RSS Df Sum of Sq F Pr(>F) 1 11 38.571 2 9 31.521 2 7.0505 1.0065 0.4032 > > beta_hat <- as.vector(coef(cracker.mod2)) > v_beta_hat <- as.matrix(vcov(cracker.mod2)) > > vec12 <- as.vector(c(0, 0, 1, -1)) > vec13 <- as.vector(c(0, 0, 2, 1)) > vec23 <- as.vector(c(0, 0, 1, 2)) > > (tau1_2 <- t(vec12) %*% beta_hat) [,1] [1,] 5.07539 > (v_tau1_2 <- t(vec12) %*% v_beta_hat %*% vec12) [,1] [1,] 1.510355 > > (tau1_3 <- t(vec13) %*% beta_hat) [,1] [1,] 12.97683 > (v_tau1_3 <- t(vec13) %*% v_beta_hat %*% vec13) [,1] [1,] 1.453528 > > (tau2_3 <- t(vec23) %*% beta_hat) [,1] [1,] 7.90144 > (v_tau2_3 <- t(vec23) %*% v_beta_hat %*% vec23) [,1] [1,] 1.413117 > > > > ### Scheffe CI's > > scheffe_cv <- sqrt((3-1)*qf(.95,2,11)) > (tau1_2 + c(-scheffe_cv*sqrt(v_tau1_2),scheffe_cv*sqrt(v_tau1_2))) [1] 1.607052 8.543728 > (tau1_3 + c(-scheffe_cv*sqrt(v_tau1_3),scheffe_cv*sqrt(v_tau1_3))) [1] 9.574367 16.379294 > (tau2_3 + c(-scheffe_cv*sqrt(v_tau2_3),scheffe_cv*sqrt(v_tau2_3))) [1] 4.546608 11.256273 > > > ### Bonferroni CI's > > t_cv <- qt(1-.05/(2*3),11) > (tau1_2 + c(-t_cv*sqrt(v_tau1_2),t_cv*sqrt(v_tau1_2))) [1] 1.609667 8.541113 > (tau1_3 + c(-t_cv*sqrt(v_tau1_3),t_cv*sqrt(v_tau1_3))) [1] 9.576932 16.376729 > (tau2_3 + c(-t_cv*sqrt(v_tau2_3),t_cv*sqrt(v_tau2_3))) [1] 4.549137 11.253744 > > >