This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
data(cars)
plot(dist ~ speed, cars, ylab = "distance")
lmod <- lm(dist ~ speed, cars)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.579 6.758 -2.60 0.012
speed 3.932 0.416 9.46 1.5e-12
n = 50, p = 2, Residual SE = 15.38, R-Squared = 0.65
abline(lmod)
lmod1 <- lm(dist ~ I(speed + rnorm(50)), cars)
coef(lmod1)
(Intercept) I(speed + rnorm(50))
-17.940 3.986
abline(lmod1, lty = 2)
lmod2 <- lm(dist ~ I(speed + 2 * rnorm(50)), cars)
coef(lmod2)
(Intercept) I(speed + 2 * rnorm(50))
-19.558 3.932
abline(lmod2, lty = 3)
lmod5 <- lm(dist ~ I(speed + 5 * rnorm(50)), cars)
coef(lmod5)
(Intercept) I(speed + 5 * rnorm(50))
13.588 1.933
abline(lmod5, lty = 4)
vv <- rep(1:5/10, each = 1000)
slopes <- numeric(5000)
for (i in 1:5000) slopes[i] <- lm(dist ~ I(speed + sqrt(vv[i]) * rnorm(50)),
cars)$coef[2]
betas <- c(coef(lmod)[2], colMeans(matrix(slopes, nrow = 1000)))
variances <- c(0, 1:5/10) + 0.5
plot(variances, betas, xlim = c(0, 1), ylim = c(3.86, 4))
gv <- lm(betas ~ variances)
coef(gv)
(Intercept) variances
3.9920 -0.1224
points(0, gv$coef[1], pch = 3)
require(simex)
Loading required package: simex
set.seed(123)
lmod <- lm(dist ~ speed, cars, x = TRUE)
simout <- simex(lmod, "speed", 0.5, B = 1000)
simout
Naive model:
lm(formula = dist ~ speed, data = cars, x = TRUE)
SIMEX-Variables: speed
Number of Simulations: 1000
Coefficients:
(Intercept) speed
-18.01 3.96
data(savings, package = "faraway")
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.566087 7.354516 3.88 0.00033
pop15 -0.461193 0.144642 -3.19 0.00260
pop75 -1.691498 1.083599 -1.56 0.12553
dpi -0.000337 0.000931 -0.36 0.71917
ddpi 0.409695 0.196197 2.09 0.04247
n = 50, p = 5, Residual SE = 3.80, R-Squared = 0.34
lmod <- lm(sr ~ pop15 + pop75 + I(dpi/1000) + ddpi, savings)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.566 7.355 3.88 0.00033
pop15 -0.461 0.145 -3.19 0.00260
pop75 -1.691 1.084 -1.56 0.12553
I(dpi/1000) -0.337 0.931 -0.36 0.71917
ddpi 0.410 0.196 2.09 0.04247
n = 50, p = 5, Residual SE = 3.80, R-Squared = 0.34
scsav <- data.frame(scale(savings))
lmod <- lm(sr ~ ., scsav)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.01e-16 1.20e-01 0.00 1.0000
pop15 -9.42e-01 2.95e-01 -3.19 0.0026
pop75 -4.87e-01 3.12e-01 -1.56 0.1255
dpi -7.45e-02 2.06e-01 -0.36 0.7192
ddpi 2.62e-01 1.26e-01 2.09 0.0425
n = 50, p = 5, Residual SE = 0.85, R-Squared = 0.34
edf <- data.frame(coef(lmod), confint(lmod))[-1, ]
names(edf) <- c("Estimate", "lb", "ub")
require(ggplot2)
p <- ggplot(aes(y = Estimate, ymin = lb, ymax = ub, x = row.names(edf)), data = edf) +
geom_pointrange()
p + coord_flip() + xlab("Predictor") + geom_hline(xint = 0, col = gray(0.75))
Error: object 'edf' not found
savings$age <- ifelse(savings$pop15 > 35, 0, 1)
savings$dpis <- (savings$dpi - mean(savings$dpi))/(2 * sd(savings$dpi))
savings$ddpis <- (savings$ddpi - mean(savings$ddpi))/(2 * sd(savings$ddpi))
sumary(lm(sr ~ age + dpis + ddpis, savings))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.82 1.01 6.75 2.2e-08
age 5.28 1.58 3.33 0.0017
dpis -1.56 1.61 -0.97 0.3361
ddpis 2.47 1.11 2.23 0.0309
n = 50, p = 4, Residual SE = 3.80, R-Squared = 0.32
data(seatpos, package = "faraway")
lmod <- lm(hipcenter ~ ., seatpos)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 436.4321 166.5716 2.62 0.014
Age 0.7757 0.5703 1.36 0.184
Weight 0.0263 0.3310 0.08 0.937
HtShoes -2.6924 9.7530 -0.28 0.784
Ht 0.6013 10.1299 0.06 0.953
Seated 0.5338 3.7619 0.14 0.888
Arm -1.3281 3.9002 -0.34 0.736
Thigh -1.1431 2.6600 -0.43 0.671
Leg -6.4390 4.7139 -1.37 0.182
n = 38, p = 9, Residual SE = 37.72, R-Squared = 0.69
round(cor(seatpos[, -9]), 2)
Age Weight HtShoes Ht Seated Arm Thigh Leg
Age 1.00 0.08 -0.08 -0.09 -0.17 0.36 0.09 -0.04
Weight 0.08 1.00 0.83 0.83 0.78 0.70 0.57 0.78
HtShoes -0.08 0.83 1.00 1.00 0.93 0.75 0.72 0.91
Ht -0.09 0.83 1.00 1.00 0.93 0.75 0.73 0.91
Seated -0.17 0.78 0.93 0.93 1.00 0.63 0.61 0.81
Arm 0.36 0.70 0.75 0.75 0.63 1.00 0.67 0.75
Thigh 0.09 0.57 0.72 0.73 0.61 0.67 1.00 0.65
Leg -0.04 0.78 0.91 0.91 0.81 0.75 0.65 1.00
x <- model.matrix(lmod)[, -1]
e <- eigen(t(x) %*% x)
e$val
[1] 3.654e+06 2.148e+04 9.043e+03 2.990e+02 1.484e+02 8.117e+01 5.336e+01
[8] 7.298e+00
sqrt(e$val[1]/e$val)
[1] 1.00 13.04 20.10 110.55 156.91 212.16 261.67 707.55
summary(lm(x[, 1] ~ x[, -1]))$r.squared
[1] 0.4995
1/(1 - 0.49948)
[1] 1.998
require(faraway)
vif(x)
Age Weight HtShoes Ht Seated Arm Thigh Leg
1.998 3.647 307.429 333.138 8.951 4.496 2.763 6.694
lmod <- lm(hipcenter + 10 * rnorm(38) ~ ., seatpos)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 456.6518 178.1780 2.56 0.016
Age 0.8215 0.6101 1.35 0.189
Weight 0.0195 0.3540 0.06 0.956
HtShoes -3.4522 10.4326 -0.33 0.743
Ht 1.0121 10.8357 0.09 0.926
Seated 0.5625 4.0240 0.14 0.890
Arm -1.0400 4.1720 -0.25 0.805
Thigh -1.8641 2.8454 -0.66 0.518
Leg -4.8880 5.0423 -0.97 0.340
n = 38, p = 9, Residual SE = 40.35, R-Squared = 0.66
round(cor(x[, 3:8]), 2)
HtShoes Ht Seated Arm Thigh Leg
HtShoes 1.00 1.00 0.93 0.75 0.72 0.91
Ht 1.00 1.00 0.93 0.75 0.73 0.91
Seated 0.93 0.93 1.00 0.63 0.61 0.81
Arm 0.75 0.75 0.63 1.00 0.67 0.75
Thigh 0.72 0.73 0.61 0.67 1.00 0.65
Leg 0.91 0.91 0.81 0.75 0.65 1.00
lmod2 <- lm(hipcenter ~ Age + Weight + Ht, seatpos)
sumary(lmod2)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 528.29773 135.31295 3.90 0.00043
Age 0.51950 0.40804 1.27 0.21159
Weight 0.00427 0.31172 0.01 0.98915
Ht -4.21190 0.99906 -4.22 0.00017
n = 38, p = 4, Residual SE = 36.49, R-Squared = 0.66
The R session information (including the OS info, R version and all packages used):
sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] graphics grDevices utils datasets methods stats base
other attached packages:
[1] simex_1.5 faraway_1.0.6 knitr_1.5 ggplot2_0.9.3.1
loaded via a namespace (and not attached):
[1] colorspace_1.2-4 dichromat_2.0-0 digest_0.6.4
[4] evaluate_0.5.3 formatR_0.10 grid_3.1.0
[7] gtable_0.1.2 labeling_0.2 MASS_7.3-31
[10] munsell_0.4.2 plyr_1.8.1 proto_0.3-10
[13] RColorBrewer_1.0-5 Rcpp_0.11.1 reshape2_1.2.2
[16] scales_0.2.3 stringr_0.6.2 tools_3.1.0
Sys.time()
[1] "2014-06-16 14:01:54 BST"