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
>
> # a simple function that works like "predict"
>
> CI.lm <- function(cur.lm, a, level=0.95, extra=0) {
+ sigma.sq <- sum(resid(cur.lm)^2) / cur.lm$df.resid
+ se <- sqrt(extra * sigma.sq + sum((a %*% vcov(cur.lm)) * a))
+ center <- sum(a*cur.lm$coef)
+ df <- cur.lm$df
+ U <- center + se * qt((1 - level)/2, df, lower.tail=F)
+ L <- center - se * qt((1 - level)/2, df, lower.tail=F)
+ return(data.frame(center, L, U))
+ }
>
> print(CI.lm(supervisor.lm, c(1,65,50,55,64,75,40)))
center L U
1 64.03723 61.14249 66.93197
> print(predict(supervisor.lm,
+ list(X1=65,X2=50,X3=55,X4=64,X5=75,X6=40),
+ interval='confidence'))
fit lwr upr
1 64.03723 61.14249 66.93197
> # these two are identical, but you can also look at intervals
> # without intercept, or where the intercept coef != 1
>
> print(CI.lm(supervisor.lm, c(1,65,50,55,64,75,40), extra=1))
center L U
1 64.03723 49.13217 78.94229
> print(predict(supervisor.lm,
+ list(X1=65,X2=50,X3=55,X4=64,X5=75,X6=40),
+ interval='prediction'))
fit lwr upr
1 64.03723 49.13217 78.94229
>
> print(CI.lm(supervisor.lm, c(0,5,5,0,0,0,0))) # how much supervisor
center L U
1 2.700687 0.9703106 4.431064
> # rating changes on
> # average if you
> # increase X1 and X2 by 5
>
> print(CI.lm(supervisor.lm, c(2,5,5,0,0,0,0))) # prediction interval
center L U
1 24.27484 -23.42877 71.97845
> # for difference in rating
> # for two supervisors who
> # differ on X1 and X2 by 5
>
>
>
> proc.time()
user system elapsed
0.600 0.024 0.628
R script