rm(list=ls()) x1 <- round(runif(10),1) x2 <- round(runif(10),1) y <- 1 + x1 - x2 + rnorm(10,sd=.1) y <- round(y,1) summary(lm(y~x1+x2)) data <- data.frame(y=y,x1=x1,x2=x2) write.table(data,file="~/Oct11.csv",sep=",",row.names=F) # read in the data data <- read.csv("~/Oct11.csv") # fit a regression model # (x=T asks R to make the design matrix available) fit <- lm(y~x1+x2,data=data,x=T) # get coefficients and standard errors summary(fit) # get an anova table anova(fit) # extract y and X y <- data$y X <- fit$x Note (in R): XtX <- t(X)%*%(X) and solve(t(X)%*%(X)) computes an inverse eye <- diag(rep(1,n)) (gives an n by n identity matrix) J <- matrix(1,n,n) gives and n by n matrix of all 1s fracJ <- (1/n)*matrix(1,n,n) givens and n by n matrix of all 1/n s. Please also note that R uses more numerically stable linear algebra algorithms to compute. For our purposes though, it is fine to do lazy things like: betahat <- solve(t(X)%*%X)%*%t(X)%*%y To do (for each, verify that your computation agree with the ones done by R): 1. compute estimated betas 2. compute estimated std errors and t-statistics 3. compute the corresponding p-values (Hint: use pt() function with lower.tail=F. Remember to mult by 2 b/c of absolute value.) 4. Make P, the projection matrix, and compute SST, SSR, SSE 5. Use "sum(diag(P))" to find p and lenght(n) to find n. 5. Compute MSE and verify that sqrt(MSE) = "Residual standard error: 0.07093" 6. Compute R-squared = SSR/SST 7. Verify that adjusted R-squared is % reduction in unexplained variance: (variance of y - estimated variance of e)/(variance of y) 8. Verify that "F-statistic: 153.6" is (SSR/(p-1))/MSE. 9. What models are being compared by that F statistic? 10 Verify that pvalue is pf(F,p-1,n-p,lower.tail=F) 11. Use matrices to compute the "Sum Sq" column in anova(fit) 12. Verify that SST = SSR(x1,x2) + SSE = SSR(x1)+SSR(x2|x1)+SSE 13. Verify that the F value column is (SSR/rank of projection matrix)/MSE 14. Verify the p-values 15. What models are being compared by the F statistics in the avova table?