This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
data(savings, package = "faraway")
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
plot(fitted(lmod), residuals(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
plot(fitted(lmod), sqrt(abs(residuals(lmod))), xlab = "Fitted", ylab = expression(sqrt(hat(epsilon))))
sumary(lm(sqrt(abs(residuals(lmod))) ~ fitted(lmod)))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.1622 0.3479 6.22 1.2e-07
fitted(lmod) -0.0614 0.0348 -1.77 0.084
n = 50, p = 2, Residual SE = 0.63, R-Squared = 0.06
par(mfrow = c(3, 3))
n <- 50
for (i in 1:9) {
x <- runif(n)
plot(x, rnorm(n))
}
for (i in 1:9) {
x <- runif(n)
plot(x, x * rnorm(n))
}
for (i in 1:9) {
x <- runif(n)
plot(x, sqrt((x)) * rnorm(n))
}
for (i in 1:9) {
x <- runif(n)
plot(x, cos(x * pi/25) + rnorm(n, sd = 1))
}
par(mfrow = c(1, 1))
plot(savings$pop15, residuals(lmod), xlab = "Population under 15", ylab = "Residuals")
abline(h = 0)
plot(savings$pop75, residuals(lmod), xlab = "Population over 75", ylab = "Residuals")
abline(h = 0)
var.test(residuals(lmod)[savings$pop15 > 35], residuals(lmod)[savings$pop15 <
35])
F test to compare two variances
data: residuals(lmod)[savings$pop15 > 35] and residuals(lmod)[savings$pop15 < 35]
F = 2.785, num df = 22, denom df = 26, p-value = 0.01358
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.241 6.430
sample estimates:
ratio of variances
2.785
data(gala, package = "faraway")
lmod <- lm(Species ~ Area + Elevation + Scruz + Nearest + Adjacent, gala)
plot(fitted(lmod), residuals(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
lmod <- lm(sqrt(Species) ~ Area + Elevation + Scruz + Nearest + Adjacent, gala)
plot(fitted(lmod), residuals(lmod), xlab = "Fitted", ylab = "Residuals")
abline(h = 0)
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
qqnorm(residuals(lmod), ylab = "Residuals", main = "")
qqline(residuals(lmod))
hist(residuals(lmod), xlab = "Residuals", main = "")
par(mfrow = c(3, 3))
n <- 50
for (i in 1:9) {
x <- rnorm(n)
qqnorm(x)
qqline(x)
}
for (i in 1:9) {
x <- exp(rnorm(n))
qqnorm(x)
qqline(x)
}
for (i in 1:9) {
x <- rcauchy(n)
qqnorm(x)
qqline(x)
}
for (i in 1:9) {
x <- runif(n)
qqnorm(x)
qqline(x)
}
par(mfrow = c(1, 1))
shapiro.test(residuals(lmod))
Shapiro-Wilk normality test
data: residuals(lmod)
W = 0.987, p-value = 0.8524
data(globwarm, package = "faraway")
lmod <- lm(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals +
mongolia + tasman, globwarm)
plot(residuals(lmod) ~ year, na.omit(globwarm), ylab = "Residuals")
abline(h = 0)
n <- length(residuals(lmod))
plot(tail(residuals(lmod), n - 1) ~ head(residuals(lmod), n - 1), xlab = expression(hat(epsilon)[i]),
ylab = expression(hat(epsilon)[i + 1]))
abline(h = 0, v = 0, col = grey(0.75))
sumary(lm(tail(residuals(lmod), n - 1) ~ head(residuals(lmod), n - 1) - 1))
Estimate Std. Error t value Pr(>|t|)
head(residuals(lmod), n - 1) 0.5951 0.0693 8.59 1.4e-14
n = 144, p = 1, Residual SE = 0.14, R-Squared = 0.34
require(lmtest)
Loading required package: lmtest
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
dwtest(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals +
mongolia + tasman, data = globwarm)
Durbin-Watson test
data: nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals + mongolia + tasman
DW = 0.8166, p-value = 1.402e-15
alternative hypothesis: true autocorrelation is greater than 0
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
hatv <- hatvalues(lmod)
head(hatv)
Australia Austria Belgium Bolivia Brazil Canada
0.06771 0.12038 0.08748 0.08947 0.06956 0.15840
sum(hatv)
[1] 5
countries <- row.names(savings)
halfnorm(hatv, labs = countries, ylab = "Leverages")
qqnorm(rstandard(lmod))
abline(0, 1)
set.seed(123)
testdata <- data.frame(x = 1:10, y = 1:10 + rnorm(10))
lmod <- lm(y ~ x, testdata)
p1 <- c(5.5, 12)
lmod1 <- lm(y ~ x, rbind(testdata, p1))
plot(y ~ x, rbind(testdata, p1))
points(5.5, 12, pch = 4, cex = 2)
abline(lmod)
abline(lmod1, lty = 2)
p2 <- c(15, 15.1)
lmod2 <- lm(y ~ x, rbind(testdata, p2))
plot(y ~ x, rbind(testdata, p2))
points(15, 15.1, pch = 4, cex = 2)
abline(lmod)
abline(lmod2, lty = 2)
p3 <- c(15, 5.1)
lmod3 <- lm(y ~ x, rbind(testdata, p3))
plot(y ~ x, rbind(testdata, p3))
points(15, 5.1, pch = 4, cex = 2)
abline(lmod)
abline(lmod3, lty = 2)
stud <- rstudent(lmod)
stud[which.max(abs(stud))]
6
2.219
qt(0.05/(50 * 2), 44)
[1] -3.526
data(star, package = "faraway")
plot(star$temp, star$light, xlab = "log(Temperature)", ylab = "log(Light Intensity)")
lmod <- lm(light ~ temp, star)
abline(lmod)
range(rstudent(lmod))
[1] -2.049 1.906
lmodi <- lm(light ~ temp, data = star, subset = (temp > 3.6))
abline(lmodi, lty = 2)
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
cook <- cooks.distance(lmod)
halfnorm(cook, 3, labs = countries, ylab = "Cook's distances")
lmodi <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings, subset = (cook < max(cook)))
sumary(lmodi)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.524046 8.224026 2.98 0.0047
pop15 -0.391440 0.157909 -2.48 0.0171
pop75 -1.280867 1.145182 -1.12 0.2694
dpi -0.000319 0.000929 -0.34 0.7331
ddpi 0.610279 0.268778 2.27 0.0281
n = 49, p = 5, Residual SE = 3.79, R-Squared = 0.36
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
plot(dfbeta(lmod)[, 2], ylab = "Change in pop15 coef")
abline(h = 0)
lmodj <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings, subset = (countries !=
"Japan"))
sumary(lmodj)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.940171 7.783997 3.08 0.0036
pop15 -0.367901 0.153630 -2.39 0.0210
pop75 -0.973674 1.155450 -0.84 0.4040
dpi -0.000471 0.000919 -0.51 0.6112
ddpi 0.334749 0.198446 1.69 0.0987
n = 49, p = 5, Residual SE = 3.74, R-Squared = 0.28
plot(lmod)
d <- residuals(lm(sr ~ pop75 + dpi + ddpi, savings))
m <- residuals(lm(pop15 ~ pop75 + dpi + ddpi, savings))
plot(m, d, xlab = "pop15 residuals", ylab = "Savings residuals")
coef(lm(d ~ m))
(Intercept) m
9.908e-17 -4.612e-01
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
coef(lmod)
(Intercept) pop15 pop75 dpi ddpi
28.5660865 -0.4611931 -1.6914977 -0.0003369 0.4096949
abline(0, coef(lmod)["pop15"])
termplot(lmod, partial.resid = TRUE, terms = 1)
mod1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings, subset = (pop15 > 35))
mod2 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings, subset = (pop15 < 35))
sumary(mod1)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.433969 21.155028 -0.12 0.91
pop15 0.273854 0.439191 0.62 0.54
pop75 -3.548477 3.033281 -1.17 0.26
dpi 0.000421 0.005000 0.08 0.93
ddpi 0.395474 0.290101 1.36 0.19
n = 23, p = 5, Residual SE = 4.45, R-Squared = 0.16
sumary(mod2)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.961795 8.083750 2.96 0.0072
pop15 -0.385898 0.195369 -1.98 0.0609
pop75 -1.327742 0.926063 -1.43 0.1657
dpi -0.000459 0.000724 -0.63 0.5326
ddpi 0.884394 0.295341 2.99 0.0067
n = 27, p = 5, Residual SE = 2.77, R-Squared = 0.51
savings$status <- ifelse(savings$pop15 > 35, "young", "old")
require(ggplot2)
ggplot(savings, aes(x = ddpi, y = sr, shape = status)) + geom_point()
ggplot(savings, aes(x = ddpi, y = sr)) + geom_point() + facet_grid(~status) +
stat_smooth(method = "lm")
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] lmtest_0.9-33 zoo_1.7-11 faraway_1.0.6 knitr_1.5
[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 lattice_0.20-29
[10] MASS_7.3-31 munsell_0.4.2 plyr_1.8.1
[13] proto_0.3-10 RColorBrewer_1.0-5 Rcpp_0.11.1
[16] reshape2_1.2.2 scales_0.2.3 stringr_0.6.2
[19] tools_3.1.0
Sys.time()
[1] "2014-06-16 14:01:38 BST"