This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
require(MASS)
Loading required package: MASS
data(savings, package = "faraway")
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
boxcox(lmod, plotit = T)
boxcox(lmod, plotit = T, lambda = seq(0.5, 1.5, by = 0.1))
data(gala, package = "faraway")
lmod <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, gala)
boxcox(lmod, lambda = seq(-0.25, 0.75, by = 0.05), plotit = T)
lmod <- lm(burntime ~ nitrogen + chlorine + potassium, leafburn)
logtrans(lmod, plotit = TRUE, alpha = seq(-min(leafburn$burntime) + 0.001, 0,
by = 0.01))
lmod1 <- lm(sr ~ pop15, savings, subset = (pop15 < 35))
lmod2 <- lm(sr ~ pop15, savings, subset = (pop15 > 35))
plot(sr ~ pop15, savings, xlab = "Pop'n under 15", ylab = "Savings Rate")
abline(v = 35, lty = 5)
segments(20, lmod1$coef[1] + lmod1$coef[2] * 20, 35, lmod1$coef[1] + lmod1$coef[2] *
35)
segments(48, lmod2$coef[1] + lmod2$coef[2] * 48, 35, lmod2$coef[1] + lmod2$coef[2] *
35)
lhs <- function(x) ifelse(x < 35, 35 - x, 0)
rhs <- function(x) ifelse(x < 35, 0, x - 35)
lmod <- lm(sr ~ lhs(pop15) + rhs(pop15), savings)
x <- seq(20, 48, by = 1)
py <- lmod$coef[1] + lmod$coef[2] * lhs(x) + lmod$coef[3] * rhs(x)
lines(x, py, lty = 2)
summary(lm(sr ~ ddpi, savings))
Call:
lm(formula = sr ~ ddpi, data = savings)
Residuals:
Min 1Q Median 3Q Max
-8.553 -3.735 0.984 2.772 9.310
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.883 1.011 7.80 4.5e-10
ddpi 0.476 0.215 2.22 0.031
Residual standard error: 4.31 on 48 degrees of freedom
Multiple R-squared: 0.0929, Adjusted R-squared: 0.074
F-statistic: 4.92 on 1 and 48 DF, p-value: 0.0314
summary(lm(sr ~ ddpi + I(ddpi^2), savings))
Call:
lm(formula = sr ~ ddpi + I(ddpi^2), data = savings)
Residuals:
Min 1Q Median 3Q Max
-8.560 -2.561 0.555 2.573 7.808
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.1304 1.4347 3.58 0.00082
ddpi 1.7575 0.5377 3.27 0.00203
I(ddpi^2) -0.0930 0.0361 -2.57 0.01326
Residual standard error: 4.08 on 47 degrees of freedom
Multiple R-squared: 0.205, Adjusted R-squared: 0.171
F-statistic: 6.06 on 2 and 47 DF, p-value: 0.00456
summary(lm(sr ~ ddpi + I(ddpi^2) + I(ddpi^3), savings))
Call:
lm(formula = sr ~ ddpi + I(ddpi^2) + I(ddpi^3), data = savings)
Residuals:
Min 1Q Median 3Q Max
-8.557 -2.557 0.562 2.576 7.798
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.145360 2.198606 2.34 0.024
ddpi 1.746017 1.380455 1.26 0.212
I(ddpi^2) -0.090967 0.225598 -0.40 0.689
I(ddpi^3) -0.000085 0.009374 -0.01 0.993
Residual standard error: 4.12 on 46 degrees of freedom
Multiple R-squared: 0.205, Adjusted R-squared: 0.153
F-statistic: 3.95 on 3 and 46 DF, p-value: 0.0137
savings <- data.frame(savings, mddpi = savings$ddpi - 10)
summary(lm(sr ~ mddpi + I(mddpi^2), savings))
Call:
lm(formula = sr ~ mddpi + I(mddpi^2), data = savings)
Residuals:
Min 1Q Median 3Q Max
-8.560 -2.561 0.555 2.573 7.808
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.4070 1.4240 9.41 2.2e-12
mddpi -0.1022 0.3027 -0.34 0.737
I(mddpi^2) -0.0930 0.0361 -2.57 0.013
Residual standard error: 4.08 on 47 degrees of freedom
Multiple R-squared: 0.205, Adjusted R-squared: 0.171
F-statistic: 6.06 on 2 and 47 DF, p-value: 0.00456
lmod <- lm(sr ~ poly(ddpi, 4), savings)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.6710 0.5846 16.54 <2e-16
poly(ddpi, 4)1 9.5590 4.1338 2.31 0.025
poly(ddpi, 4)2 -10.4999 4.1338 -2.54 0.015
poly(ddpi, 4)3 -0.0374 4.1338 -0.01 0.993
poly(ddpi, 4)4 3.6120 4.1338 0.87 0.387
n = 50, p = 5, Residual SE = 4.13, R-Squared = 0.22
lmod <- lm(sr ~ polym(pop15, ddpi, degree = 2), savings)
pop15r <- seq(20, 50, len = 10)
ddpir <- seq(0, 20, len = 10)
pgrid <- expand.grid(pop15 = pop15r, ddpi = ddpir)
pv <- predict(lmod, pgrid)
persp(pop15r, ddpir, matrix(pv, 10, 10), theta = 45, xlab = "Pop under 15",
ylab = "Growth", zlab = "Savings rate", ticktype = "detailed", shade = 0.25)
funky <- function(x) sin(2 * pi * x^3)^3
x <- seq(0, 1, by = 0.01)
y <- funky(x) + 0.1 * rnorm(101)
matplot(x, cbind(y, funky(x)), type = "pl", ylab = "y", pch = 20, lty = 1, col = 1)
g4 <- lm(y ~ poly(x, 4))
g12 <- lm(y ~ poly(x, 12))
matplot(x, cbind(y, g4$fit, g12$fit), type = "pll", ylab = "y", lty = c(1, 2),
pch = 20, col = 1)
require(splines)
Loading required package: splines
knots <- c(0, 0, 0, 0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 1, 1, 1, 1)
bx <- splineDesign(knots, x)
lmodb <- lm(y ~ bx - 1)
matplot(x, bx, type = "l", col = 1)
matplot(x, cbind(y, lmodb$fit), type = "pl", ylab = "y", pch = 20, lty = 1,
col = 1)
ssf <- smooth.spline(x, y)
matplot(x, cbind(y, ssf$y), type = "pl", ylab = "y", lty = 1, pch = 20, col = 1)
require(mgcv)
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.7-29. For overview type 'help("mgcv-package")'.
gamod <- gam(sr ~ s(pop15) + s(pop75) + s(dpi) + s(ddpi), data = savings)
plot(gamod)
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] splines graphics grDevices utils datasets methods stats
[8] base
other attached packages:
[1] mgcv_1.7-29 nlme_3.1-117 MASS_7.3-31 faraway_1.0.6
[5] 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:12 BST"