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?