tests/testthat/data/rapp.test.1/packrat/lib-R/MASS/scripts/ch08.R
In rappster/rapp: Runtime environment for package development and application deployment in the context of the rapp framework
#-*- R -*-
## Script from Fourth Edition of `Modern Applied Statistics with S'
# Chapter 8 Non-linear and Smooth Regression
library(MASS)
library(lattice)
options(width=65, digits=5, height=9999)
pdf(file="ch08.pdf", width=8, height=6, pointsize=9)
# From Chapter 6, for comparisons
set.seed(123)
cpus.samp <-
c(3, 5, 6, 7, 8, 10, 11, 16, 20, 21, 22, 23, 24, 25, 29, 33, 39, 41, 44, 45,
46, 49, 57, 58, 62, 63, 65, 66, 68, 69, 73, 74, 75, 76, 78, 83, 86,
88, 98, 99, 100, 103, 107, 110, 112, 113, 115, 118, 119, 120, 122,
124, 125, 126, 127, 132, 136, 141, 144, 146, 147, 148, 149, 150, 151,
152, 154, 156, 157, 158, 159, 160, 161, 163, 166, 167, 169, 170, 173,
174, 175, 176, 177, 183, 184, 187, 188, 189, 194, 195, 196, 197, 198,
199, 202, 204, 205, 206, 208, 209)
cpus1 <- cpus
attach(cpus)
for(v in names(cpus)[2:7])
cpus1[[v]] <- cut(cpus[[v]], unique(quantile(cpus[[v]])),
include.lowest = TRUE)
detach()
cpus.lm <- lm(log10(perf) ~ ., data=cpus1[cpus.samp, 2:8])
cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
res2 <- log10(cpus1[-cpus.samp, "perf"]) -
predict(cpus.lm2, cpus1[-cpus.samp,])
cpus2 <- cpus[, 2:8] # excludes names, authors' predictions
cpus2[, 1:3] <- log10(cpus2[, 1:3])
test.cpus <- function(fit)
sqrt(sum((log10(cpus2[-cpus.samp, "perf"]) -
predict(fit, cpus2[-cpus.samp,]))^2)/109)
# 8.1 An introductory example
attach(wtloss)
# alter margin 4; others are default
oldpar <- par(mar = c(5.1, 4.1, 4.1, 4.1))
plot(Days, Weight, type = "p", ylab = "Weight (kg)")
Wt.lbs <- pretty(range(Weight*2.205))
axis(side = 4, at = Wt.lbs/2.205, lab = Wt.lbs, srt = 90)
mtext("Weight (lb)", side = 4, line = 3)
par(oldpar) # restore settings
detach()
# 8.2 Fitting non-linear regression models
wtloss.st <- c(b0 = 90, b1 = 95, th = 120)
wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th),
data = wtloss, start = wtloss.st, trace = TRUE)
wtloss.fm
expn <- function(b0, b1, th, x) {
temp <- 2^(-x/th)
model.func <- b0 + b1 * temp
Z <- cbind(1, temp, (b1 * x * temp * log(2))/th^2)
dimnames(Z) <- list(NULL, c("b0", "b1", "th"))
attr(model.func, "gradient") <- Z
model.func
}
wtloss.gr <- nls(Weight ~ expn(b0, b1, th, Days),
data = wtloss, start = wtloss.st, trace = TRUE)
expn1 <- deriv(y ~ b0 + b1 * 2^(-x/th), c("b0", "b1", "th"),
function(b0, b1, th, x) {})
negexp <- selfStart(model = ~ b0 + b1*exp(-x/th),
initial = negexp.SSival, parameters = c("b0", "b1", "th"),
template = function(x, b0, b1, th) {})
wtloss.ss <- nls(Weight ~ negexp(Days, B0, B1, theta),
data = wtloss, trace = TRUE)
# 8.3 Non-linear fitted model objects and method functions
summary(wtloss.gr)
deviance(wtloss.gr)
vcov(wtloss.gr)
A <- model.matrix(~ Strip - 1, data = muscle)
rats.nls1 <- nls(log(Length) ~ cbind(A, rho^Conc),
data = muscle, start = c(rho = 0.1), algorithm = "plinear")
(B <- coef(rats.nls1))
st <- list(alpha = B[2:22], beta = B[23], rho = B[1])
rats.nls2 <- nls(log(Length) ~ alpha[Strip] + beta*rho^Conc,
data = muscle, start = st)
attach(muscle)
Muscle <- expand.grid(Conc = sort(unique(Conc)),
Strip = levels(Strip))
Muscle$Yhat <- predict(rats.nls2, Muscle)
Muscle$logLength <- rep(NA, nrow(Muscle))
ind <- match(paste(Strip, Conc),
paste(Muscle$Strip, Muscle$Conc))
Muscle$logLength[ind] <- log(Length)
detach()
xyplot(Yhat ~ Conc | Strip, Muscle, as.table = TRUE,
ylim = range(c(Muscle$Yhat, Muscle$logLength), na.rm = TRUE),
subscripts = TRUE, xlab = "Calcium Chloride concentration (mM)",
ylab = "log(Length in mm)", panel =
function(x, y, subscripts, ...) {
panel.xyplot(x, Muscle$logLength[subscripts], ...)
llines(spline(x, y))
})
# 8.5 Confidence intervals for parameters
expn2 <- deriv(~b0 + b1*((w0 - b0)/b1)^(x/d0),
c("b0","b1","d0"), function(b0, b1, d0, x, w0) {})
wtloss.init <- function(obj, w0) {
p <- coef(obj)
d0 <- - log((w0 - p["b0"])/p["b1"])/log(2) * p["th"]
c(p[c("b0", "b1")], d0 = as.vector(d0))
}
out <- NULL
w0s <- c(110, 100, 90)
for(w0 in w0s) {
fm <- nls(Weight ~ expn2(b0, b1, d0, Days, w0),
wtloss, start = wtloss.init(wtloss.gr, w0))
out <- rbind(out, c(coef(fm)["d0"], confint(fm, "d0")))
}
dimnames(out)[[1]] <- paste(w0s,"kg:")
out
fm0 <- lm(Wt*Time ~ Viscosity + Time - 1, data = stormer)
b0 <- coef(fm0)
names(b0) <- c("b1", "b2")
b0
storm.fm <- nls(Time ~ b1*Viscosity/(Wt-b2), data = stormer,
start = b0, trace = TRUE)
bc <- coef(storm.fm)
se <- sqrt(diag(vcov(storm.fm)))
dv <- deviance(storm.fm)
par(pty = "s")
b1 <- bc[1] + seq(-3*se[1], 3*se[1], length = 51)
b2 <- bc[2] + seq(-3*se[2], 3*se[2], length = 51)
bv <- expand.grid(b1, b2)
attach(stormer)
ssq <- function(b)
sum((Time - b[1] * Viscosity/(Wt-b[2]))^2)
dbetas <- apply(bv, 1, ssq)
cc <- matrix(Time - rep(bv[,1],rep(23, 2601)) *
Viscosity/(Wt - rep(bv[,2], rep(23, 2601))), 23)
dbetas <- matrix(drop(rep(1, 23) %*% cc^2), 51)
fstat <- matrix( ((dbetas - dv)/2) / (dv/21), 51, 51)
qf(0.95, 2, 21)
plot(b1, b2, type = "n")
lev <- c(1, 2, 5, 7, 10, 15, 20)
contour(b1, b2, fstat, levels = lev, labex = 0.75, lty = 2, add = TRUE)
contour(b1, b2, fstat, levels = qf(0.95,2,21), add = TRUE, labex = 0)
text(31.6, 0.3, labels = "95% CR", adj = 0, cex = 0.75)
points(bc[1], bc[2], pch = 3, mkh = 0.1)
detach()
par(pty = "m")
library(boot)
storm.fm <- nls(Time ~ b*Viscosity/(Wt - c), stormer,
start = c(b=29.401, c=2.2183))
summary(storm.fm)$parameters
st <- cbind(stormer, fit=fitted(storm.fm))
storm.bf <- function(rs, i) {
# st <- st # for S-PLUS
st$Time <- st$fit + rs[i]
coef(nls(Time ~ b * Viscosity/(Wt - c), st,
start = coef(storm.fm)))
}
rs <- scale(resid(storm.fm), scale = FALSE) # remove the mean
(storm.boot <- boot(rs, storm.bf, R = 9999)) ## slow
boot.ci(storm.boot, index = 1,
type = c("norm", "basic", "perc", "bca"))
boot.ci(storm.boot, index = 2,
type = c("norm", "basic", "perc", "bca"))
# 8.5 Assessing the linear approximation
opar <- par(pty = "m", mfrow = c(1, 3))
plot(profile(update(wtloss.gr, trace = FALSE)))
par(opar)
# 8.7 One-dimensional curve fitting
attach(GAGurine)
par(mfrow = c(3, 2))
plot(Age, GAG, main = "Degree 6 polynomial")
GAG.lm <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) +
I(Age^5) + I(Age^6) + I(Age^7) + I(Age^8))
anova(GAG.lm)
GAG.lm2 <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) +
I(Age^5) + I(Age^6))
xx <- seq(0, 17, len = 200)
lines(xx, predict(GAG.lm2, data.frame(Age = xx)))
library(splines)
plot(Age, GAG, type = "n", main = "Splines")
lines(Age, fitted(lm(GAG ~ ns(Age, df = 5))))
lines(Age, fitted(lm(GAG ~ ns(Age, df = 10))), lty = 3)
lines(Age, fitted(lm(GAG ~ ns(Age, df = 20))), lty = 4)
lines(smooth.spline(Age, GAG), lwd = 3)
legend(12, 50, c("df=5", "df=10", "df=20", "Smoothing"),
lty = c(1, 3, 4, 1), lwd = c(1,1,1,3), bty = "n")
plot(Age, GAG, type = "n", main = "loess")
lines(loess.smooth(Age, GAG))
plot(Age, GAG, type = "n", main = "supsmu")
lines(supsmu(Age, GAG))
lines(supsmu(Age, GAG, bass = 3), lty = 3)
lines(supsmu(Age, GAG, bass = 10), lty = 4)
legend(12, 50, c("default", "base = 3", "base = 10"),
lty = c(1, 3, 4), bty = "n")
plot(Age, GAG, type = "n", main = "ksmooth")
lines(ksmooth(Age, GAG, "normal", bandwidth = 1))
lines(ksmooth(Age, GAG, "normal", bandwidth = 5), lty = 3)
legend(12, 50, c("width = 1", "width = 5"), lty = c(1, 3), bty = "n")
library(KernSmooth)
plot(Age, GAG, type = "n", main = "locpoly")
(h <- dpill(Age, GAG))
lines(locpoly(Age, GAG, degree = 0, bandwidth = h))
lines(locpoly(Age, GAG, degree = 1, bandwidth = h), lty = 3)
lines(locpoly(Age, GAG, degree = 2, bandwidth = h), lty = 4)
legend(12, 50, c("const", "linear", "quadratic"),
lty = c(1, 3, 4), bty = "n")
detach()
# 8.8 Additive models
## R has a different gam() in package mgcv
library(mgcv)
rock.lm <- lm(log(perm) ~ area + peri + shape, data = rock)
summary(rock.lm)
(rock.gam <- gam(log(perm) ~ s(area) + s(peri) + s(shape), data=rock))
#summary(rock.gam)
#anova(rock.lm, rock.gam)
par(mfrow = c(2, 3), pty = "s")
plot(rock.gam, se = TRUE, pages = 0)
rock.gam1 <- gam(log(perm) ~ area + peri + s(shape), data = rock)
plot(rock.gam1, se = TRUE)
par(pty="m")
#anova(rock.lm, rock.gam1, rock.gam)
library(mda)
rock.bruto <- bruto(rock[, -4], rock[, 4])
rock.bruto$type
rock.bruto$df
Xin <- as.matrix(cpus2[cpus.samp, 1:6])
test2 <- function(fit) {
Xp <- as.matrix(cpus2[-cpus.samp, 1:6])
sqrt(sum((log10(cpus2[-cpus.samp, "perf"]) -
predict(fit, Xp))^2)/109)
}
cpus.bruto <- bruto(Xin, log10(cpus2[cpus.samp, 7]))
test2(cpus.bruto)
cpus.bruto$type
cpus.bruto$df
# examine the fitted functions
par(mfrow = c(3, 2))
Xp <- matrix(sapply(cpus2[cpus.samp, 1:6], mean), 100, 6, byrow = TRUE)
for(i in 1:6) {
xr <- sapply(cpus2, range)
Xp1 <- Xp; Xp1[, i] <- seq(xr[1, i], xr[2, i], len = 100)
Xf <- predict(cpus.bruto, Xp1)
plot(Xp1[ ,i], Xf, xlab=names(cpus2)[i], ylab= "", type = "l")
}
cpus.mars <- mars(Xin, log10(cpus2[cpus.samp,7]))
showcuts <- function(obj)
{
tmp <- obj$cuts[obj$sel, ]
dimnames(tmp) <- list(NULL, dimnames(Xin)[[2]])
tmp
}
showcuts(cpus.mars)
test2(cpus.mars)
# examine the fitted functions
Xp <- matrix(sapply(cpus2[cpus.samp, 1:6], mean), 100, 6, byrow = TRUE)
for(i in 1:6) {
xr <- sapply(cpus2, range)
Xp1 <- Xp; Xp1[, i] <- seq(xr[1, i], xr[2, i], len = 100)
Xf <- predict(cpus.mars, Xp1)
plot(Xp1[ ,i], Xf, xlab = names(cpus2)[i], ylab = "", type = "l")
}
cpus.mars2 <- mars(Xin, log10(cpus2[cpus.samp,7]), degree = 2)
showcuts(cpus.mars2)
test2(cpus.mars2)
cpus.mars6 <- mars(Xin, log10(cpus2[cpus.samp,7]), degree = 6)
showcuts(cpus.mars6)
test2(cpus.mars6)
if(require(acepack)) {
attach(cpus2)
cpus.avas <- avas(cpus2[, 1:6], perf)
plot(log10(perf), cpus.avas$ty)
par(mfrow = c(2, 3))
for(i in 1:6) {
o <- order(cpus2[, i])
plot(cpus2[o, i], cpus.avas$tx[o, i], type = "l",
xlab = names(cpus2[i]), ylab = "")
}
detach()
}
# 8.9 Projection-pursuit regression
attach(rock)
rock1 <- data.frame(area = area/10000, peri = peri/10000,
shape = shape, perm = perm)
detach()
(rock.ppr <- ppr(log(perm) ~ area + peri + shape, data = rock1,
nterms = 2, max.terms = 5))
rock.ppr
summary(rock.ppr)
par(mfrow = c(3, 2))
plot(rock.ppr)
plot(update(rock.ppr, bass = 5))
plot(rock.ppr2 <- update(rock.ppr, sm.method = "gcv", gcvpen = 2))
par(mfrow = c(1, 1))
summary(rock.ppr2)
summary(rock1) # to find the ranges of the variables
Xp <- expand.grid(area = seq(0.1, 1.2, 0.05),
peri = seq(0, 0.5, 0.02), shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.ppr2, Xp))
wireframe(fit ~ area + peri, rock.grid, screen = list(z=160,x=-60),
aspect = c(1, 0.5), drape = TRUE)
# or
persp(seq(0.1, 1.2, 0.05), seq(0, 0.5, 0.02), matrix(rock.grid$fit, 23),
d = 5, theta = -160, phi = 30, zlim = c(-1, 15))
(cpus.ppr <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 2, max.terms = 10, bass = 5))
cpus.ppr <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 8, max.terms = 10, bass = 5)
test.cpus(cpus.ppr)
ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 2, max.terms = 10, sm.method = "spline")
cpus.ppr2 <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 7, max.terms = 10, sm.method = "spline")
test.cpus(cpus.ppr2)
res3 <- log10(cpus2[-cpus.samp, "perf"]) -
predict(cpus.ppr, cpus2[-cpus.samp,])
wilcox.test(res2^2, res3^2, paired = TRUE, alternative = "greater")
# 8.10 Neural networks
library(nnet)
attach(rock)
area1 <- area/10000; peri1 <- peri/10000
rock1 <- data.frame(perm, area = area1, peri = peri1, shape)
rock.nn <- nnet(log(perm) ~ area + peri + shape, rock1,
size = 3, decay = 1e-3, linout = TRUE, skip = TRUE,
maxit = 1000, Hess = TRUE)
sum((log(perm) - predict(rock.nn))^2)
detach()
eigen(rock.nn$Hessian, TRUE)$values # rock.nn$Hessian in R
Xp <- expand.grid(area = seq(0.1, 1.2, 0.05),
peri = seq(0, 0.5, 0.02), shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.nn, Xp))
wireframe(fit ~ area + peri, rock.grid, screen = list(z=160, x=-60),
aspect = c(1, 0.5), drape = TRUE)
# or
persp(seq(0.1, 1.2, 0.05), seq(0, 0.5, 0.02), matrix(rock.grid$fit, 23),
d = 5, theta = -160, phi = 30, zlim = c(-1, 15))
attach(cpus2)
cpus3 <-
data.frame(syct= syct-2, mmin=mmin-3, mmax=mmax-4, cach=cach/256,
chmin=chmin/100, chmax=chmax/100, perf=perf)
detach()
test.cpus <- function(fit)
sqrt(sum((log10(cpus3[-cpus.samp, "perf"]) -
predict(fit, cpus3[-cpus.samp,]))^2)/109)
cpus.nn1 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,],
linout = TRUE, skip = TRUE, size = 0)
test.cpus(cpus.nn1)
cpus.nn2 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 4, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn2)
cpus.nn3 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 10, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn3)
cpus.nn4 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 25, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn4)
CVnn.cpus <- function(formula, data = cpus3[cpus.samp, ],
size = c(0, 4, 4, 10, 10),
lambda = c(0, rep(c(0.003, 0.01), 2)),
nreps = 5, nifold = 10, ...)
{
CVnn1 <- function(formula, data, nreps=1, ri, ...)
{
truth <- log10(data$perf)
res <- numeric(length(truth))
cat(" fold")
for (i in sort(unique(ri))) {
cat(" ", i, sep="")
for(rep in 1:nreps) {
learn <- nnet(formula, data[ri !=i,], trace=FALSE, ...)
res[ri == i] <- res[ri == i] +
predict(learn, data[ri == i,])
}
}
cat("\n")
sum((truth - res/nreps)^2)
}
choice <- numeric(length(lambda))
ri <- sample(nifold, nrow(data), replace = TRUE)
for(j in seq(along=lambda)) {
cat(" size =", size[j], "decay =", lambda[j], "\n")
choice[j] <- CVnn1(formula, data, nreps=nreps, ri=ri,
size=size[j], decay=lambda[j], ...)
}
cbind(size=size, decay=lambda, fit=sqrt(choice/100))
}
CVnn.cpus(log10(perf) ~ ., data = cpus3[cpus.samp,],
linout = TRUE, skip = TRUE, maxit = 1000)
# End of ch08
rappster/rapp documentation built on May 26, 2019, 11:56 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#-*- R -*-
## Script from Fourth Edition of `Modern Applied Statistics with S'
# Chapter 8 Non-linear and Smooth Regression
library(MASS)
library(lattice)
options(width=65, digits=5, height=9999)
pdf(file="ch08.pdf", width=8, height=6, pointsize=9)
# From Chapter 6, for comparisons
set.seed(123)
cpus.samp <-
c(3, 5, 6, 7, 8, 10, 11, 16, 20, 21, 22, 23, 24, 25, 29, 33, 39, 41, 44, 45,
46, 49, 57, 58, 62, 63, 65, 66, 68, 69, 73, 74, 75, 76, 78, 83, 86,
88, 98, 99, 100, 103, 107, 110, 112, 113, 115, 118, 119, 120, 122,
124, 125, 126, 127, 132, 136, 141, 144, 146, 147, 148, 149, 150, 151,
152, 154, 156, 157, 158, 159, 160, 161, 163, 166, 167, 169, 170, 173,
174, 175, 176, 177, 183, 184, 187, 188, 189, 194, 195, 196, 197, 198,
199, 202, 204, 205, 206, 208, 209)
cpus1 <- cpus
attach(cpus)
for(v in names(cpus)[2:7])
cpus1[[v]] <- cut(cpus[[v]], unique(quantile(cpus[[v]])),
include.lowest = TRUE)
detach()
cpus.lm <- lm(log10(perf) ~ ., data=cpus1[cpus.samp, 2:8])
cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
res2 <- log10(cpus1[-cpus.samp, "perf"]) -
predict(cpus.lm2, cpus1[-cpus.samp,])
cpus2 <- cpus[, 2:8] # excludes names, authors' predictions
cpus2[, 1:3] <- log10(cpus2[, 1:3])
test.cpus <- function(fit)
sqrt(sum((log10(cpus2[-cpus.samp, "perf"]) -
predict(fit, cpus2[-cpus.samp,]))^2)/109)
# 8.1 An introductory example
attach(wtloss)
# alter margin 4; others are default
oldpar <- par(mar = c(5.1, 4.1, 4.1, 4.1))
plot(Days, Weight, type = "p", ylab = "Weight (kg)")
Wt.lbs <- pretty(range(Weight*2.205))
axis(side = 4, at = Wt.lbs/2.205, lab = Wt.lbs, srt = 90)
mtext("Weight (lb)", side = 4, line = 3)
par(oldpar) # restore settings
detach()
# 8.2 Fitting non-linear regression models
wtloss.st <- c(b0 = 90, b1 = 95, th = 120)
wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th),
data = wtloss, start = wtloss.st, trace = TRUE)
wtloss.fm
expn <- function(b0, b1, th, x) {
temp <- 2^(-x/th)
model.func <- b0 + b1 * temp
Z <- cbind(1, temp, (b1 * x * temp * log(2))/th^2)
dimnames(Z) <- list(NULL, c("b0", "b1", "th"))
attr(model.func, "gradient") <- Z
model.func
}
wtloss.gr <- nls(Weight ~ expn(b0, b1, th, Days),
data = wtloss, start = wtloss.st, trace = TRUE)
expn1 <- deriv(y ~ b0 + b1 * 2^(-x/th), c("b0", "b1", "th"),
function(b0, b1, th, x) {})
negexp <- selfStart(model = ~ b0 + b1*exp(-x/th),
initial = negexp.SSival, parameters = c("b0", "b1", "th"),
template = function(x, b0, b1, th) {})
wtloss.ss <- nls(Weight ~ negexp(Days, B0, B1, theta),
data = wtloss, trace = TRUE)
# 8.3 Non-linear fitted model objects and method functions
summary(wtloss.gr)
deviance(wtloss.gr)
vcov(wtloss.gr)
A <- model.matrix(~ Strip - 1, data = muscle)
rats.nls1 <- nls(log(Length) ~ cbind(A, rho^Conc),
data = muscle, start = c(rho = 0.1), algorithm = "plinear")
(B <- coef(rats.nls1))
st <- list(alpha = B[2:22], beta = B[23], rho = B[1])
rats.nls2 <- nls(log(Length) ~ alpha[Strip] + beta*rho^Conc,
data = muscle, start = st)
attach(muscle)
Muscle <- expand.grid(Conc = sort(unique(Conc)),
Strip = levels(Strip))
Muscle$Yhat <- predict(rats.nls2, Muscle)
Muscle$logLength <- rep(NA, nrow(Muscle))
ind <- match(paste(Strip, Conc),
paste(Muscle$Strip, Muscle$Conc))
Muscle$logLength[ind] <- log(Length)
detach()
xyplot(Yhat ~ Conc | Strip, Muscle, as.table = TRUE,
ylim = range(c(Muscle$Yhat, Muscle$logLength), na.rm = TRUE),
subscripts = TRUE, xlab = "Calcium Chloride concentration (mM)",
ylab = "log(Length in mm)", panel =
function(x, y, subscripts, ...) {
panel.xyplot(x, Muscle$logLength[subscripts], ...)
llines(spline(x, y))
})
# 8.5 Confidence intervals for parameters
expn2 <- deriv(~b0 + b1*((w0 - b0)/b1)^(x/d0),
c("b0","b1","d0"), function(b0, b1, d0, x, w0) {})
wtloss.init <- function(obj, w0) {
p <- coef(obj)
d0 <- - log((w0 - p["b0"])/p["b1"])/log(2) * p["th"]
c(p[c("b0", "b1")], d0 = as.vector(d0))
}
out <- NULL
w0s <- c(110, 100, 90)
for(w0 in w0s) {
fm <- nls(Weight ~ expn2(b0, b1, d0, Days, w0),
wtloss, start = wtloss.init(wtloss.gr, w0))
out <- rbind(out, c(coef(fm)["d0"], confint(fm, "d0")))
}
dimnames(out)[[1]] <- paste(w0s,"kg:")
out
fm0 <- lm(Wt*Time ~ Viscosity + Time - 1, data = stormer)
b0 <- coef(fm0)
names(b0) <- c("b1", "b2")
b0
storm.fm <- nls(Time ~ b1*Viscosity/(Wt-b2), data = stormer,
start = b0, trace = TRUE)
bc <- coef(storm.fm)
se <- sqrt(diag(vcov(storm.fm)))
dv <- deviance(storm.fm)
par(pty = "s")
b1 <- bc[1] + seq(-3*se[1], 3*se[1], length = 51)
b2 <- bc[2] + seq(-3*se[2], 3*se[2], length = 51)
bv <- expand.grid(b1, b2)
attach(stormer)
ssq <- function(b)
sum((Time - b[1] * Viscosity/(Wt-b[2]))^2)
dbetas <- apply(bv, 1, ssq)
cc <- matrix(Time - rep(bv[,1],rep(23, 2601)) *
Viscosity/(Wt - rep(bv[,2], rep(23, 2601))), 23)
dbetas <- matrix(drop(rep(1, 23) %*% cc^2), 51)
fstat <- matrix( ((dbetas - dv)/2) / (dv/21), 51, 51)
qf(0.95, 2, 21)
plot(b1, b2, type = "n")
lev <- c(1, 2, 5, 7, 10, 15, 20)
contour(b1, b2, fstat, levels = lev, labex = 0.75, lty = 2, add = TRUE)
contour(b1, b2, fstat, levels = qf(0.95,2,21), add = TRUE, labex = 0)
text(31.6, 0.3, labels = "95% CR", adj = 0, cex = 0.75)
points(bc[1], bc[2], pch = 3, mkh = 0.1)
detach()
par(pty = "m")
library(boot)
storm.fm <- nls(Time ~ b*Viscosity/(Wt - c), stormer,
start = c(b=29.401, c=2.2183))
summary(storm.fm)$parameters
st <- cbind(stormer, fit=fitted(storm.fm))
storm.bf <- function(rs, i) {
# st <- st # for S-PLUS
st$Time <- st$fit + rs[i]
coef(nls(Time ~ b * Viscosity/(Wt - c), st,
start = coef(storm.fm)))
}
rs <- scale(resid(storm.fm), scale = FALSE) # remove the mean
(storm.boot <- boot(rs, storm.bf, R = 9999)) ## slow
boot.ci(storm.boot, index = 1,
type = c("norm", "basic", "perc", "bca"))
boot.ci(storm.boot, index = 2,
type = c("norm", "basic", "perc", "bca"))
# 8.5 Assessing the linear approximation
opar <- par(pty = "m", mfrow = c(1, 3))
plot(profile(update(wtloss.gr, trace = FALSE)))
par(opar)
# 8.7 One-dimensional curve fitting
attach(GAGurine)
par(mfrow = c(3, 2))
plot(Age, GAG, main = "Degree 6 polynomial")
GAG.lm <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) +
I(Age^5) + I(Age^6) + I(Age^7) + I(Age^8))
anova(GAG.lm)
GAG.lm2 <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) +
I(Age^5) + I(Age^6))
xx <- seq(0, 17, len = 200)
lines(xx, predict(GAG.lm2, data.frame(Age = xx)))
library(splines)
plot(Age, GAG, type = "n", main = "Splines")
lines(Age, fitted(lm(GAG ~ ns(Age, df = 5))))
lines(Age, fitted(lm(GAG ~ ns(Age, df = 10))), lty = 3)
lines(Age, fitted(lm(GAG ~ ns(Age, df = 20))), lty = 4)
lines(smooth.spline(Age, GAG), lwd = 3)
legend(12, 50, c("df=5", "df=10", "df=20", "Smoothing"),
lty = c(1, 3, 4, 1), lwd = c(1,1,1,3), bty = "n")
plot(Age, GAG, type = "n", main = "loess")
lines(loess.smooth(Age, GAG))
plot(Age, GAG, type = "n", main = "supsmu")
lines(supsmu(Age, GAG))
lines(supsmu(Age, GAG, bass = 3), lty = 3)
lines(supsmu(Age, GAG, bass = 10), lty = 4)
legend(12, 50, c("default", "base = 3", "base = 10"),
lty = c(1, 3, 4), bty = "n")
plot(Age, GAG, type = "n", main = "ksmooth")
lines(ksmooth(Age, GAG, "normal", bandwidth = 1))
lines(ksmooth(Age, GAG, "normal", bandwidth = 5), lty = 3)
legend(12, 50, c("width = 1", "width = 5"), lty = c(1, 3), bty = "n")
library(KernSmooth)
plot(Age, GAG, type = "n", main = "locpoly")
(h <- dpill(Age, GAG))
lines(locpoly(Age, GAG, degree = 0, bandwidth = h))
lines(locpoly(Age, GAG, degree = 1, bandwidth = h), lty = 3)
lines(locpoly(Age, GAG, degree = 2, bandwidth = h), lty = 4)
legend(12, 50, c("const", "linear", "quadratic"),
lty = c(1, 3, 4), bty = "n")
detach()
# 8.8 Additive models
## R has a different gam() in package mgcv
library(mgcv)
rock.lm <- lm(log(perm) ~ area + peri + shape, data = rock)
summary(rock.lm)
(rock.gam <- gam(log(perm) ~ s(area) + s(peri) + s(shape), data=rock))
#summary(rock.gam)
#anova(rock.lm, rock.gam)
par(mfrow = c(2, 3), pty = "s")
plot(rock.gam, se = TRUE, pages = 0)
rock.gam1 <- gam(log(perm) ~ area + peri + s(shape), data = rock)
plot(rock.gam1, se = TRUE)
par(pty="m")
#anova(rock.lm, rock.gam1, rock.gam)
library(mda)
rock.bruto <- bruto(rock[, -4], rock[, 4])
rock.bruto$type
rock.bruto$df
Xin <- as.matrix(cpus2[cpus.samp, 1:6])
test2 <- function(fit) {
Xp <- as.matrix(cpus2[-cpus.samp, 1:6])
sqrt(sum((log10(cpus2[-cpus.samp, "perf"]) -
predict(fit, Xp))^2)/109)
}
cpus.bruto <- bruto(Xin, log10(cpus2[cpus.samp, 7]))
test2(cpus.bruto)
cpus.bruto$type
cpus.bruto$df
# examine the fitted functions
par(mfrow = c(3, 2))
Xp <- matrix(sapply(cpus2[cpus.samp, 1:6], mean), 100, 6, byrow = TRUE)
for(i in 1:6) {
xr <- sapply(cpus2, range)
Xp1 <- Xp; Xp1[, i] <- seq(xr[1, i], xr[2, i], len = 100)
Xf <- predict(cpus.bruto, Xp1)
plot(Xp1[ ,i], Xf, xlab=names(cpus2)[i], ylab= "", type = "l")
}
cpus.mars <- mars(Xin, log10(cpus2[cpus.samp,7]))
showcuts <- function(obj)
{
tmp <- obj$cuts[obj$sel, ]
dimnames(tmp) <- list(NULL, dimnames(Xin)[[2]])
tmp
}
showcuts(cpus.mars)
test2(cpus.mars)
# examine the fitted functions
Xp <- matrix(sapply(cpus2[cpus.samp, 1:6], mean), 100, 6, byrow = TRUE)
for(i in 1:6) {
xr <- sapply(cpus2, range)
Xp1 <- Xp; Xp1[, i] <- seq(xr[1, i], xr[2, i], len = 100)
Xf <- predict(cpus.mars, Xp1)
plot(Xp1[ ,i], Xf, xlab = names(cpus2)[i], ylab = "", type = "l")
}
cpus.mars2 <- mars(Xin, log10(cpus2[cpus.samp,7]), degree = 2)
showcuts(cpus.mars2)
test2(cpus.mars2)
cpus.mars6 <- mars(Xin, log10(cpus2[cpus.samp,7]), degree = 6)
showcuts(cpus.mars6)
test2(cpus.mars6)
if(require(acepack)) {
attach(cpus2)
cpus.avas <- avas(cpus2[, 1:6], perf)
plot(log10(perf), cpus.avas$ty)
par(mfrow = c(2, 3))
for(i in 1:6) {
o <- order(cpus2[, i])
plot(cpus2[o, i], cpus.avas$tx[o, i], type = "l",
xlab = names(cpus2[i]), ylab = "")
}
detach()
}
# 8.9 Projection-pursuit regression
attach(rock)
rock1 <- data.frame(area = area/10000, peri = peri/10000,
shape = shape, perm = perm)
detach()
(rock.ppr <- ppr(log(perm) ~ area + peri + shape, data = rock1,
nterms = 2, max.terms = 5))
rock.ppr
summary(rock.ppr)
par(mfrow = c(3, 2))
plot(rock.ppr)
plot(update(rock.ppr, bass = 5))
plot(rock.ppr2 <- update(rock.ppr, sm.method = "gcv", gcvpen = 2))
par(mfrow = c(1, 1))
summary(rock.ppr2)
summary(rock1) # to find the ranges of the variables
Xp <- expand.grid(area = seq(0.1, 1.2, 0.05),
peri = seq(0, 0.5, 0.02), shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.ppr2, Xp))
wireframe(fit ~ area + peri, rock.grid, screen = list(z=160,x=-60),
aspect = c(1, 0.5), drape = TRUE)
# or
persp(seq(0.1, 1.2, 0.05), seq(0, 0.5, 0.02), matrix(rock.grid$fit, 23),
d = 5, theta = -160, phi = 30, zlim = c(-1, 15))
(cpus.ppr <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 2, max.terms = 10, bass = 5))
cpus.ppr <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 8, max.terms = 10, bass = 5)
test.cpus(cpus.ppr)
ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 2, max.terms = 10, sm.method = "spline")
cpus.ppr2 <- ppr(log10(perf) ~ ., data = cpus2[cpus.samp,],
nterms = 7, max.terms = 10, sm.method = "spline")
test.cpus(cpus.ppr2)
res3 <- log10(cpus2[-cpus.samp, "perf"]) -
predict(cpus.ppr, cpus2[-cpus.samp,])
wilcox.test(res2^2, res3^2, paired = TRUE, alternative = "greater")
# 8.10 Neural networks
library(nnet)
attach(rock)
area1 <- area/10000; peri1 <- peri/10000
rock1 <- data.frame(perm, area = area1, peri = peri1, shape)
rock.nn <- nnet(log(perm) ~ area + peri + shape, rock1,
size = 3, decay = 1e-3, linout = TRUE, skip = TRUE,
maxit = 1000, Hess = TRUE)
sum((log(perm) - predict(rock.nn))^2)
detach()
eigen(rock.nn$Hessian, TRUE)$values # rock.nn$Hessian in R
Xp <- expand.grid(area = seq(0.1, 1.2, 0.05),
peri = seq(0, 0.5, 0.02), shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.nn, Xp))
wireframe(fit ~ area + peri, rock.grid, screen = list(z=160, x=-60),
aspect = c(1, 0.5), drape = TRUE)
# or
persp(seq(0.1, 1.2, 0.05), seq(0, 0.5, 0.02), matrix(rock.grid$fit, 23),
d = 5, theta = -160, phi = 30, zlim = c(-1, 15))
attach(cpus2)
cpus3 <-
data.frame(syct= syct-2, mmin=mmin-3, mmax=mmax-4, cach=cach/256,
chmin=chmin/100, chmax=chmax/100, perf=perf)
detach()
test.cpus <- function(fit)
sqrt(sum((log10(cpus3[-cpus.samp, "perf"]) -
predict(fit, cpus3[-cpus.samp,]))^2)/109)
cpus.nn1 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,],
linout = TRUE, skip = TRUE, size = 0)
test.cpus(cpus.nn1)
cpus.nn2 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 4, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn2)
cpus.nn3 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 10, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn3)
cpus.nn4 <- nnet(log10(perf) ~ ., cpus3[cpus.samp,], linout = TRUE,
skip = TRUE, size = 25, decay = 0.01, maxit = 1000)
test.cpus(cpus.nn4)
CVnn.cpus <- function(formula, data = cpus3[cpus.samp, ],
size = c(0, 4, 4, 10, 10),
lambda = c(0, rep(c(0.003, 0.01), 2)),
nreps = 5, nifold = 10, ...)
{
CVnn1 <- function(formula, data, nreps=1, ri, ...)
{
truth <- log10(data$perf)
res <- numeric(length(truth))
cat(" fold")
for (i in sort(unique(ri))) {
cat(" ", i, sep="")
for(rep in 1:nreps) {
learn <- nnet(formula, data[ri !=i,], trace=FALSE, ...)
res[ri == i] <- res[ri == i] +
predict(learn, data[ri == i,])
}
}
cat("\n")
sum((truth - res/nreps)^2)
}
choice <- numeric(length(lambda))
ri <- sample(nifold, nrow(data), replace = TRUE)
for(j in seq(along=lambda)) {
cat(" size =", size[j], "decay =", lambda[j], "\n")
choice[j] <- CVnn1(formula, data, nreps=nreps, ri=ri,
size=size[j], decay=lambda[j], ...)
}
cbind(size=size, decay=lambda, fit=sqrt(choice/100))
}
CVnn.cpus(log10(perf) ~ ., data = cpus3[cpus.samp,],
linout = TRUE, skip = TRUE, maxit = 1000)
# End of ch08
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.