> ### ANOVA and F-test > (SSTO <- t(Y) %*% (I_n - J_n) %*% Y) ### Total Sum of Squares [,1] [1,] 13771.58 > (df.TO <- n-1) ### Total degrees of freedom [1] 368 > (SSE <- t(Y) %*% (I_n - H) %*% Y) ### Error Sum of Squares [,1] [1,] 614.4363 > (df.E <- n-p) ### Error degrees of freedom [1] 365 > (MSE <- SSE/df.E) ### s^2 = MSE [,1] [1,] 1.683387 > (SSR <- t(Y) %*% (H - J_n) %*% Y) ### Regression Sum of Squares [,1] [1,] 13157.14 > (df.R <- p-1) ### Regression degrees of freedom [1] 3 > (MSR <- SSR/df.R) ### Regression Mean Square [,1] [1,] 4385.714 > (F.star <- MSR/MSE) ### F* statistic [,1] [1,] 2605.292 > (F.95 <- qf(.95,df.R,df.E)) ### Critical F-value [1] 2.629365 > (F.pv <- 1-pf(F.star,df.R,df.E)) ### P-value [,1] [1,] 0 > > > ### Inference on Betas > (b <- solve(t(X) %*% X) %*% t(X) %*% Y) ### Reg Coeffs [,1] X0 124.92215741 ELEV.C -0.09216969 LAT -2.28717895 LONG -0.06222776 > (s2.b <- MSE[1,1] * solve(t(X) %*% X)) ### Var-Cov Matrix for b X0 ELEV.C LAT LONG X0 42.81601556 0.0872962403 -0.1696193390 -0.391521475 ELEV.C 0.08729624 0.0002058061 -0.0003714998 -0.000794204 LAT -0.16961934 -0.0003714998 0.0016047703 0.001259702 LONG -0.39152147 -0.0007942040 0.0012597017 0.003672220 > (s.b <- sqrt(diag(s2.b))) ### SE's of b X0 ELEV.C LAT LONG 6.54339480 0.01434595 0.04005958 0.06059885 > (t.b <- b/s.b) ### t-stats for Beta's [,1] X0 19.09134 ELEV.C -6.42479 LAT -57.09443 LONG -1.02688 > (pv.b <- 2*(1-pt(abs(t.b),n-p))) ### P-values for t-tests [,1] X0 0.000000e+00 ELEV.C 4.124359e-10 LAT 0.000000e+00 LONG 3.051569e-01 > print(round(cbind(b,s.b,t.b,pv.b),4)) s.b X0 124.9222 6.5434 19.0913 0.0000 ELEV.C -0.0922 0.0143 -6.4248 0.0000 LAT -2.2872 0.0401 -57.0944 0.0000 LONG -0.0622 0.0606 -1.0269 0.3052 > ci.b.lo <- b + qt(0.025,df.E)*s.b ### Lower Bound of CI for Beta > ci.b.hi <- b + qt(0.975,df.E)*s.b ### Upper Bound of CI for Beta > print(round(cbind(ci.b.lo,ci.b.hi),3)) [,1] [,2] X0 112.055 137.790 ELEV.C -0.120 -0.064 LAT -2.366 -2.208 LONG -0.181 0.057 > > ### CI for Mean and PI for Individual @ X-h' = [1,10,30,100] > x_h <- matrix(c(1,10,30,100),ncol=1) ### Generate x_h vector > (x_h.b <- t(x_h) %*% b) ### Predicted value [,1] [1,] 49.16232 > (s2.x_h.b <- t(x_h) %*% s2.b %*% x_h) ### s2{yhat_h) [,1] [1,] 0.0144642 > (s2.pred <- MSE + s2.x_h.b) ### s2{pred} [,1] [1,] 1.697851 > (s.x_h.b <- sqrt(s2.x_h.b)) ### s{yhat_h) [,1] [1,] 0.1202672 > (s.pred <- sqrt(s2.pred)) ### s{pred} [,1] [1,] 1.303016 > (x_h.beta.ci <- x_h.b + qt(c(.025,.975),df.E)*s.x_h.b) ### 95% CI For Mean [1] 48.92581 49.39882 > (x_h.beta.pi <- x_h.b + qt(c(.025,.975),df.E)*s.pred) ### 95% PI For Individual [1] 46.59996 51.72468 > > > #### Fit Regression Model using lm function > tx.mod1 <- lm(Mean.Jan ~ ELEV.C + LAT + LONG) > summary(tx.mod1) Call: lm(formula = Mean.Jan ~ ELEV.C + LAT + LONG) Residuals: Min 1Q Median 3Q Max -3.6464 -0.8679 -0.1357 0.7990 4.5759 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 124.92216 6.54339 19.091 < 2e-16 *** ELEV.C -0.09217 0.01435 -6.425 4.12e-10 *** LAT -2.28718 0.04006 -57.094 < 2e-16 *** LONG -0.06223 0.06060 -1.027 0.305 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.297 on 365 degrees of freedom Multiple R-squared: 0.9554, Adjusted R-squared: 0.955 F-statistic: 2605 on 3 and 365 DF, p-value: < 2.2e-16 > anova(tx.mod1) Analysis of Variance Table Response: Mean.Jan Df Sum Sq Mean Sq F value Pr(>F) ELEV.C 1 5785.3 5785.3 3436.6972 <2e-16 *** LAT 1 7370.1 7370.1 4378.1235 <2e-16 *** LONG 1 1.8 1.8 1.0545 0.3052 Residuals 365 614.4 1.7 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > drop1(tx.mod1) Single term deletions Model: Mean.Jan ~ ELEV.C + LAT + LONG Df Sum of Sq RSS AIC 614.4 196.16 ELEV.C 1 69.5 683.9 233.69 LAT 1 5487.5 6101.9 1041.25 LONG 1 1.8 616.2 195.22 > confint(tx.mod1) 2.5 % 97.5 % (Intercept) 112.0546723 137.78964250 ELEV.C -0.1203808 -0.06395861 LAT -2.3659555 -2.20840240 LONG -0.1813945 0.05693894 > >