TXnorm1a <- read.csv("http://www.stat.ufl.edu/~winner/data/TXnorm1a.csv", header=TRUE) attach(TXnorm1a); names(TXnorm1a) X <- as.matrix(cbind(X0,ELEV.C,LAT,LONG)) Y <- as.matrix(Mean.Jan) n <- length(Mean.Jan) p <- ncol(X) ### ncol(X) = # of columns for X J_n <- matrix(rep(1/n,n*n),ncol=n) ### (1/n)J H <- X %*% solve(t(X) %*% X) %*% t(X) ### H I_n <- diag(n) ### I ### ANOVA and F-test (SSTO <- t(Y) %*% (I_n - J_n) %*% Y) ### Total Sum of Squares (df.TO <- n-1) ### Total degrees of freedom (SSE <- t(Y) %*% (I_n - H) %*% Y) ### Error Sum of Squares (df.E <- n-p) ### Error degrees of freedom (MSE <- SSE/df.E) ### s^2 = MSE (SSR <- t(Y) %*% (H - J_n) %*% Y) ### Regression Sum of Squares (df.R <- p-1) ### Regression degrees of freedom (MSR <- SSR/df.R) ### Regression Mean Square (F.star <- MSR/MSE) ### F* statistic (F.95 <- qf(.95,df.R,df.E)) ### Critical F-value (F.pv <- 1-pf(F.star,df.R,df.E)) ### P-value ### Inference on Betas (b <- solve(t(X) %*% X) %*% t(X) %*% Y) ### Reg Coeffs (s2.b <- MSE[1,1] * solve(t(X) %*% X)) ### Var-Cov Matrix for b (s.b <- sqrt(diag(s2.b))) ### SE's of b (t.b <- b/s.b) ### t-stats for Beta's (pv.b <- 2*(1-pt(abs(t.b),n-p))) ### P-values for t-tests print(round(cbind(b,s.b,t.b,pv.b),4)) 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)) ### 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 (s2.x_h.b <- t(x_h) %*% s2.b %*% x_h) ### s2{yhat_h) (s2.pred <- MSE + s2.x_h.b) ### s2{pred} (s.x_h.b <- sqrt(s2.x_h.b)) ### s{yhat_h) (s.pred <- sqrt(s2.pred)) ### s{pred} (x_h.beta.ci <- x_h.b + qt(c(.025,.975),df.E)*s.x_h.b) ### 95% CI For Mean (x_h.beta.pi <- x_h.b + qt(c(.025,.975),df.E)*s.pred) ### 95% PI For Individual #### Fit Regression Model using lm function tx.mod1 <- lm(Mean.Jan ~ ELEV.C + LAT + LONG) summary(tx.mod1) anova(tx.mod1) drop1(tx.mod1) confint(tx.mod1)