This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
data(globwarm, package = "faraway")
lmod <- lm(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals +
mongolia + tasman, globwarm)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.24256 0.02701 -8.98 2e-15
wusa 0.07738 0.04293 1.80 0.07365
jasper -0.22879 0.07811 -2.93 0.00399
westgreen 0.00958 0.04184 0.23 0.81917
chesapeake -0.03211 0.03405 -0.94 0.34735
tornetrask 0.09267 0.04505 2.06 0.04161
urals 0.18537 0.09143 2.03 0.04457
mongolia 0.04197 0.04579 0.92 0.36100
tasman 0.11545 0.03011 3.83 0.00019
n = 145, p = 9, Residual SE = 0.18, R-Squared = 0.48
cor(residuals(lmod)[-1], residuals(lmod)[-length(residuals(lmod))])
[1] 0.5833
require(nlme)
Loading required package: nlme
glmod <- gls(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask +
urals + mongolia + tasman, correlation = corAR1(form = ~year), data = na.omit(globwarm))
summary(glmod)
Generalized least squares fit by REML
Model: nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals + mongolia + tasman
Data: na.omit(globwarm)
AIC BIC logLik
-108.2 -76.17 65.1
Correlation Structure: AR(1)
Formula: ~year
Parameter estimate(s):
Phi
0.711
Coefficients:
Value Std.Error t-value p-value
(Intercept) -0.23011 0.06702 -3.433 0.0008
wusa 0.06674 0.09877 0.676 0.5004
jasper -0.20244 0.18803 -1.077 0.2835
westgreen -0.00440 0.08985 -0.049 0.9610
chesapeake -0.00735 0.07350 -0.100 0.9205
tornetrask 0.03835 0.09483 0.404 0.6865
urals 0.24142 0.22871 1.056 0.2930
mongolia 0.05695 0.10490 0.543 0.5881
tasman 0.12035 0.07457 1.614 0.1089
Correlation:
(Intr) wusa jasper wstgrn chespk trntrs urals mongol
wusa -0.517
jasper -0.058 -0.299
westgreen 0.330 -0.533 0.121
chesapeake 0.090 -0.314 0.230 0.147
tornetrask -0.430 0.499 -0.197 -0.328 -0.441
urals -0.110 -0.142 -0.265 0.075 -0.064 -0.346
mongolia 0.459 -0.437 -0.205 0.217 0.449 -0.343 -0.371
tasman 0.037 -0.322 0.065 0.134 0.116 -0.434 0.416 -0.017
Standardized residuals:
Min Q1 Med Q3 Max
-2.31123 -0.53484 0.02343 0.50016 2.97225
Residual standard error: 0.2046
Degrees of freedom: 145 total; 136 residual
intervals(glmod, which = "var-cov")
Approximate 95% confidence intervals
Correlation structure:
lower est. upper
Phi 0.51 0.711 0.8384
attr(,"label")
[1] "Correlation structure:"
Residual standard error:
lower est. upper
0.1541 0.2046 0.2716
glmod <- gls(yield ~ variety, oatvar, correlation = corCompSymm(form = ~1 |
block))
intervals(glmod)
Approximate 95% confidence intervals
Coefficients:
lower est. upper
(Intercept) 291.543 334.4 377.2570
variety2 -4.904 42.2 89.3039
variety3 -18.904 28.2 75.3039
variety4 -94.704 -47.6 -0.4961
variety5 57.896 105.0 152.1039
variety6 -50.904 -3.8 43.3039
variety7 -63.104 -16.0 31.1039
variety8 2.696 49.8 96.9039
attr(,"label")
[1] "Coefficients:"
Correlation structure:
lower est. upper
Rho 0.06596 0.396 0.7494
attr(,"label")
[1] "Correlation structure:"
Residual standard error:
lower est. upper
33.39 47.05 66.28
data(fpe, package = "faraway")
fpe
EI A B C D E F G H J K A2 B2 N
Ain 260 51 64 36 23 9 5 4 4 3 3 105 114 17
Alpes 75 14 17 9 9 3 1 2 1 1 1 32 31 5
Ariege 107 27 18 13 17 2 2 2 1 1 1 57 33 6
Bouches.du.Rhone 1036 191 204 119 205 29 13 13 10 10 6 466 364 30
Charente.Maritime 367 71 76 47 37 8 34 5 4 4 2 163 142 17
Cotes.du.Nord 396 93 90 57 54 13 5 9 4 3 5 193 155 15
Drome 257 57 55 31 30 10 4 5 4 3 3 116 99 13
Finistere 595 132 149 95 49 21 9 11 6 5 10 249 259 21
Gironde 735 195 137 98 83 20 16 13 13 8 5 356 261 29
Indre 181 34 39 28 28 4 3 4 3 2 1 82 72 8
Landes 219 62 47 31 26 5 3 3 3 2 1 107 84 8
Loire.Atlantique 653 149 156 94 49 23 15 13 10 7 8 274 275 25
Lozere 58 10 18 9 4 2 0 1 1 1 1 20 29 2
Marne 145 32 33 20 15 4 2 3 2 2 1 63 59 8
Morbihan 414 86 117 65 33 14 6 8 5 4 4 162 190 10
Oise 416 87 88 59 62 13 7 10 6 5 3 192 160 12
Pyrenees.Atlantique 391 90 91 66 33 12 6 6 5 4 3 165 168 17
Rhin 413 75 125 58 19 17 6 8 6 6 4 138 204 18
Sarthe 346 72 87 49 40 10 6 8 4 3 3 149 145 12
Seine.Maritime 783 171 181 91 123 24 13 18 10 7 6 370 297 23
Sevres 240 54 66 34 16 8 7 5 3 4 2 98 108 7
Val.D.Oise 533 111 100 74 81 22 12 10 7 7 6 252 192 14
Vendee 336 61 105 59 19 10 11 6 5 4 3 115 176 8
Yonne 216 44 52 31 24 7 4 4 3 3 2 91 91 8
lmod <- lm(A2 ~ A + B + C + D + E + F + G + H + J + K + N - 1, fpe, weights = 1/EI)
coef(lmod)
A B C D E F G H J
1.0671 -0.1051 0.2460 0.9262 0.2494 0.7551 1.9722 -0.5662 0.6116
K N
1.2107 0.5294
lm(A2 ~ A + B + C + D + E + F + G + H + J + K + N - 1, fpe)$coef
A B C D E F G H J
1.0751 -0.1246 0.2574 0.9045 0.6707 0.7825 2.1657 -0.8543 0.1444
K N
0.5181 0.5583
lm(A2 ~ A + B + C + D + E + F + G + H + J + K + N - 1, fpe, weights = 53/EI)$coef
A B C D E F G H J
1.0671 -0.1051 0.2460 0.9262 0.2494 0.7551 1.9722 -0.5662 0.6116
K N
1.2107 0.5294
lm(A2 ~ offset(A + G + K) + C + D + E + F + N - 1, fpe, weights = 1/EI)$coef
C D E F N
0.2258 0.9700 0.3902 0.7442 0.6085
require(mgcv)
Loading required package: mgcv
This is mgcv 1.7-29. For overview type 'help("mgcv-package")'.
M <- list(w = 1/fpe$EI, X = model.matrix(lmod), y = fpe$A2, Ain = rbind(diag(11),
-diag(11)), C = matrix(0, 0, 0), array(0, 0), S = list(), off = NULL, p = rep(0.5,
11), bin = c(rep(0, 11), rep(-1, 11)))
a <- pcls(M)
names(a) <- colnames(model.matrix(lmod))
round(a, 3)
A B C D E F G H J K N
1.000 0.000 0.208 0.969 0.359 0.743 1.000 0.367 0.000 1.000 0.575
lmod <- lm(dist ~ speed, cars)
plot(residuals(lmod) ~ speed, cars)
wlmod <- gls(dist ~ speed, data = cars, weight = varConstPower(1, form = ~speed))
summary(wlmod)
Generalized least squares fit by REML
Model: dist ~ speed
Data: cars
AIC BIC logLik
412.8 422.2 -201.4
Variance function:
Structure: Constant plus power of variance covariate
Formula: ~speed
Parameter estimates:
const power
3.160 1.022
Coefficients:
Value Std.Error t-value p-value
(Intercept) -11.085 4.052 -2.736 0.0087
speed 3.484 0.320 10.880 0.0000
Correlation:
(Intr)
speed -0.9
Standardized residuals:
Min Q1 Med Q3 Max
-1.4521 -0.6898 -0.1308 0.6375 3.0757
Residual standard error: 0.7637
Degrees of freedom: 50 total; 48 residual
data(corrosion, package = "faraway")
plot(loss ~ Fe, corrosion, xlab = "Iron content", ylab = "Weight loss")
lmod <- lm(loss ~ Fe, corrosion)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 129.79 1.40 92.5 < 2e-16
Fe -24.02 1.28 -18.8 1.1e-09
n = 13, p = 2, Residual SE = 3.06, R-Squared = 0.97
abline(coef(lmod))
lmoda <- lm(loss ~ factor(Fe), corrosion)
points(corrosion$Fe, fitted(lmoda), pch = 3)
anova(lmod, lmoda)
Analysis of Variance Table
Model 1: loss ~ Fe
Model 2: loss ~ factor(Fe)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11 102.9
2 6 11.8 5 91.1 9.28 0.0086
lmodp <- lm(loss ~ Fe + I(Fe^2) + I(Fe^3) + I(Fe^4) + I(Fe^5) + I(Fe^6), corrosion)
plot(loss ~ Fe, data = corrosion, ylim = c(60, 130))
points(corrosion$Fe, fitted(lmoda), pch = 3)
grid <- seq(0, 2, len = 50)
lines(grid, predict(lmodp, data.frame(Fe = grid)))
summary(lmodp)$r.squared
[1] 0.9965
data(gala, package = "faraway")
lsmod <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala)
sumary(lsmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.06822 19.15420 0.37 0.7154
Area -0.02394 0.02242 -1.07 0.2963
Elevation 0.31946 0.05366 5.95 3.8e-06
Nearest 0.00914 1.05414 0.01 0.9932
Scruz -0.24052 0.21540 -1.12 0.2752
Adjacent -0.07480 0.01770 -4.23 0.0003
n = 30, p = 6, Residual SE = 60.98, R-Squared = 0.77
require(MASS)
Loading required package: MASS
rlmod <- rlm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala)
summary(rlmod)
Call: rlm(formula = Species ~ Area + Elevation + Nearest + Scruz +
Adjacent, data = gala)
Residuals:
Min 1Q Median 3Q Max
-74.39 -18.35 -6.36 21.19 229.08
Coefficients:
Value Std. Error t value
(Intercept) 6.361 12.390 0.513
Area -0.006 0.015 -0.421
Elevation 0.248 0.035 7.132
Nearest 0.359 0.682 0.527
Scruz -0.195 0.139 -1.401
Adjacent -0.055 0.011 -4.765
Residual standard error: 29.7 on 24 degrees of freedom
wts <- rlmod$w
names(wts) <- row.names(gala)
head(sort(wts), 10)
SantaCruz SantaMaria SanCristobal Pinta Gardner1
0.1746 0.3078 0.4142 0.5376 0.6613
Espanola Gardner2 Baltra Bartolome Caldwell
0.6795 0.8500 1.0000 1.0000 1.0000
require(quantreg)
Loading required package: quantreg
Loading required package: SparseM
Attaching package: 'SparseM'
The following object is masked from 'package:base':
backsolve
l1mod <- rq(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data = gala)
summary(l1mod)
Call: rq(formula = Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
data = gala)
tau: [1] 0.5
Coefficients:
coefficients lower bd upper bd
(Intercept) 1.31445 -19.87777 24.37411
Area -0.00306 -0.03185 0.52800
Elevation 0.23211 0.12453 0.50196
Nearest 0.16366 -3.16339 2.98896
Scruz -0.12314 -0.47987 0.13476
Adjacent -0.05185 -0.10458 0.01739
set.seed(123)
ltsmod <- ltsreg(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala)
coef(ltsmod)
(Intercept) Area Elevation Nearest Scruz Adjacent
12.50668 1.54536 0.01673 0.52349 -0.09407 -0.14259
ltsmod <- ltsreg(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala,
nsamp = "exact")
coef(ltsmod)
(Intercept) Area Elevation Nearest Scruz Adjacent
9.38115 1.54366 0.02412 0.81111 -0.11773 -0.19792
bcoef <- matrix(0, 1000, 6)
for (i in 1:1000) {
newy <- predict(ltsmod) + residuals(ltsmod)[sample(30, rep = T)]
brg <- ltsreg(newy ~ Area + Elevation + Nearest + Scruz + Adjacent, gala,
nsamp = "best")
bcoef[i, ] <- brg$coef
}
colnames(bcoef) <- names(coef(ltsmod))
apply(bcoef, 2, function(x) quantile(x, c(0.025, 0.975)))
(Intercept) Area Elevation Nearest Scruz Adjacent
2.5% 1.918 1.494 -0.01462 0.1588 -0.26238 -0.2203
97.5% 21.467 1.607 0.07333 1.9147 0.09822 -0.1270
require(ggplot2)
bcoef <- data.frame(bcoef)
p1 <- ggplot(bcoef, aes(x = Area)) + geom_density() + xlim(1.45, 1.65)
p1 + geom_vline(xintercept = c(1.4976, 1.623), linetype = "dashed")
Warning: Removed 25 rows containing non-finite values (stat_density).
p2 <- ggplot(bcoef, aes(x = Adjacent)) + geom_density() + xlim(-0.25, -0.13)
p2 + geom_vline(xintercept = c(-0.23375, -0.15138), linetype = "dashed")
Warning: Removed 39 rows containing non-finite values (stat_density).
limod <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala, subset = (row.names(gala) !=
"Isabela"))
sumary(limod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.5861 13.4019 1.69 0.1055
Area 0.2957 0.0619 4.78 8.0e-05
Elevation 0.1404 0.0497 2.82 0.0096
Nearest -0.2552 0.7217 -0.35 0.7269
Scruz -0.0901 0.1498 -0.60 0.5534
Adjacent -0.0650 0.0122 -5.32 2.1e-05
n = 29, p = 6, Residual SE = 41.65, R-Squared = 0.87
data(star, package = "faraway")
plot(light ~ temp, star)
gs1 <- lm(light ~ temp, star)
abline(coef(gs1))
gs2 <- rlm(light ~ temp, star)
abline(coef(gs2), lty = 2)
gs3 <- ltsreg(light ~ temp, star, nsamp = "exact")
abline(coef(gs3), lty = 5)
i <- order(pipeline$Field)
npipe <- pipeline[i, ]
ff <- gl(12, 9)[-108]
meanfield <- unlist(lapply(split(npipe$Field, ff), mean))
varlab <- unlist(lapply(split(npipe$Lab, ff), var))
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] quantreg_5.05 SparseM_1.03 MASS_7.3-31 mgcv_1.7-29
[5] nlme_3.1-117 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 lattice_0.20-29
[10] Matrix_1.1-3 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:02:09 BST"