R version 2.9.2 (2009-08-24)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>
> # Read in new data
>
> supervisor.table <- read.table('http://www-stat.stanford.edu/~jtaylo/courses/stats191/data/supervisor.table', header=T)
> attach(supervisor.table)
>
> # original regression model
>
> supervisor.lm <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6)
> print(summary(supervisor.lm))
Call:
lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6)
Residuals:
Min 1Q Median 3Q Max
-10.9418 -4.3555 0.3158 5.5425 11.5990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.78708 11.58926 0.931 0.361634
X1 0.61319 0.16098 3.809 0.000903 ***
X2 -0.07305 0.13572 -0.538 0.595594
X3 0.32033 0.16852 1.901 0.069925 .
X4 0.08173 0.22148 0.369 0.715480
X5 0.03838 0.14700 0.261 0.796334
X6 -0.21706 0.17821 -1.218 0.235577
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.240e-05
>
> n <- length(Y)
> X <- cbind(rep(1,n), X1, X2, X3, X4, X5, X6)
> colnames(X)[1] <- '(Intercept)'
>
> # solve for beta.hat
>
> beta <- solve(t(X) %*% X) %*% t(X) %*% Y
>
> # covariance matrix
>
> Y.hat <- X %*% beta
>
> sigma.hat <- sqrt(sum((Y - Y.hat)^2) / (n - ncol(X)))
> cov.beta <- sigma.hat^2 * solve(t(X) %*% X)
> print(cov.beta)
(Intercept) X1 X2 X3 X4
(Intercept) 134.31088348 -1.270544e-01 -0.181523035 -0.307411947 -0.216292873
X1 -0.12705442 2.591556e-02 -0.008174628 -0.008184570 -0.018572460
X2 -0.18152304 -8.174628e-03 0.018421192 -0.003122857 0.002306358
X3 -0.30741195 -8.184570e-03 -0.003122857 0.028399098 -0.007689886
X4 -0.21629287 -1.857246e-02 0.002306358 -0.007689886 0.049052362
X5 -1.13204004 2.336757e-05 -0.000457418 0.004847088 -0.009130728
X6 0.03278243 1.153912e-02 -0.004464340 -0.010427982 -0.016854411
X5 X6
(Intercept) -1.132040e+00 0.032782433
X1 2.336757e-05 0.011539122
X2 -4.574180e-04 -0.004464340
X3 4.847088e-03 -0.010427982
X4 -9.130728e-03 -0.016854411
X5 2.160766e-02 -0.003349602
X6 -3.349602e-03 0.031758616
> print(sqrt(diag(cov.beta)))
(Intercept) X1 X2 X3 X4 X5
11.5892572 0.1609831 0.1357247 0.1685203 0.2214777 0.1469954
X6
0.1782095
>
> # an easier way to get the covariance matrix
>
> print(vcov(supervisor.lm))
(Intercept) X1 X2 X3 X4
(Intercept) 134.31088348 -1.270544e-01 -0.181523035 -0.307411947 -0.216292873
X1 -0.12705442 2.591556e-02 -0.008174628 -0.008184570 -0.018572460
X2 -0.18152304 -8.174628e-03 0.018421192 -0.003122857 0.002306358
X3 -0.30741195 -8.184570e-03 -0.003122857 0.028399098 -0.007689886
X4 -0.21629287 -1.857246e-02 0.002306358 -0.007689886 0.049052362
X5 -1.13204004 2.336757e-05 -0.000457418 0.004847088 -0.009130728
X6 0.03278243 1.153912e-02 -0.004464340 -0.010427982 -0.016854411
X5 X6
(Intercept) -1.132040e+00 0.032782433
X1 2.336757e-05 0.011539122
X2 -4.574180e-04 -0.004464340
X3 4.847088e-03 -0.010427982
X4 -9.130728e-03 -0.016854411
X5 2.160766e-02 -0.003349602
X6 -3.349602e-03 0.031758616
>
>
>
> proc.time()
user system elapsed
0.572 0.016 0.586
R script