Nothing
#-*- R -*-
## Script from Fourth Edition of `Modern Applied Statistics with S'
# Chapter 4 Graphical Output
library(MASS)
library(lattice)
pdf(file="ch04.pdf", width=8, height=6, pointsize=9)
options(width=65, digits=5)
# 4.2 Basic plotting functions
lung.deaths <- aggregate(ts.union(mdeaths, fdeaths), 1)
barplot(t(lung.deaths), names = dimnames(lung.deaths)[[1]],
main = "UK deaths from lung disease")
if(interactive())
legend(locator(1), c("Males", "Females"), fill = c(2, 3))
loc <- barplot(t(lung.deaths), names = dimnames(lung.deaths)[[1]],
angle = c(45, 135), density = 10, col = 1)
total <- rowSums(lung.deaths)
text(loc, total + par("cxy")[2], total, cex = 0.7, xpd = TRUE)
# S: if(interactive()) brush(hills)
topo.loess <- loess(z ~ x * y, topo, degree = 2, span = 0.25)
topo.mar <- list(x = seq(0, 6.5, 0.2), y=seq(0, 6.5, 0.2))
topo.lo <- predict(topo.loess, expand.grid(topo.mar))
par(pty = "s") # square plot
contour(topo.mar$x, topo.mar$y, topo.lo, xlab = "", ylab = "",
levels = seq(700,1000,25), cex = 0.7)
points(topo$x, topo$y)
par(pty = "m")
topo.lo1 <- cbind(expand.grid(x=topo.mar$x, y=topo.mar$y),
z=as.vector(topo.lo))
contourplot(z ~ x * y, topo.lo1, aspect = 1,
at = seq(700, 1000, 25), xlab = "", ylab = "",
panel = function(x, y, subscripts, ...) {
panel.levelplot(x, y, subscripts, ...)
panel.xyplot(topo$x,topo$y, cex = 0.5)
}
)
# see help(Skye)
# ternary(Skye/100, ord = c(1, 3, 2))
# 4.3 Enhancing plots
attach(wtloss)
oldpar <- par(no.readonly = TRUE)
# alter margin 4; others are default
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, las = 0)
mtext("Weight (lb)", side = 4, line = 3)
detach()
par(oldpar)
x <- 0:100
plik <- function(lambda)
sum(dpois(x, lambda) * 2 * ( (lambda - x) +
x * log(pmax(1, x)/lambda)))
lambda <- c(1e-8, 0.05, seq(0.1, 5, 0.1))
plot(lambda, sapply(lambda, plik), type = "l", ylim = c(0, 1.4),
xlab = expression(lambda),
ylab = expression(paste(E[lambda], "(deviance)")))
abline(h = 1, lty = 3)
# 4.4 Fine control of graphics
## in R just use swiss
# swiss <- data.frame(Fertility = swiss.fertility, swiss.x)
attach(swiss)
qqnorm(Infant.Mortality)
qqline(Infant.Mortality)
samp <- cbind(Infant.Mortality, matrix(rnorm(47*19), 47, 19))
samp <- apply(scale(samp), 2, sort)
rs <- samp[, 1]
xs <- qqnorm(rs, plot = FALSE)$x
env <- t(apply(samp[, -1], 1, range))
matplot(xs, cbind(rs, env), type = "pnn",
pch = 4, mkh = 0.06, axes = FALSE, xlab = "", ylab = "")
xyul <- par("usr")
smidge <- min(diff(c(xyul[1], xs, xyul[2])))/2
segments(xs - smidge, env[, 1], xs + smidge, env[, 1])
segments(xs - smidge, env[, 2], xs + smidge, env[, 2])
xul <- trunc(10*xyul[1:2])/10
axis(1, at=seq(xul[1], xul[2], by=0.1), labels = FALSE, tck=0.01)
xi <- trunc(xyul[1:2])
axis(1, at = seq(xi[1], xi[2], by = 0.5), tck = 0.02)
yul <- trunc(5*xyul[3:4])/5
axis(2, at=seq(yul[1], yul[2], by=0.2), labels = FALSE, tck=0.01)
yi <- trunc(xyul[3:4])
axis(2, at = yi[1]:yi[2], tck = 0.02)
box(bty = "l") # lower case "L"
# ps.options()$fonts
# R cannot change font family in a plot.
mtext("Quantiles of Standard Normal", side=1, line=2.5, font=3)
mtext(expression(R[i]), side = 2, line = 2, at = yul[2])
detach()
# 4.5 Trellis graphics
xyplot(time ~ dist, data = hills,
panel = function(x, y, ...) {
panel.xyplot(x, y, ...)
panel.lmline(x, y, type = "l")
panel.abline(lqs(y ~ x), lty = 3)
# identify(x, y, row.names(hills))
}
)
## note: don't use separate title() call
bwplot(Expt ~ Speed, data = michelson, ylab = "Experiment No.",
main = "Speed of Light Data")
splom(~ swiss, aspect = "fill",
panel = function(x, y, ...) {
panel.xyplot(x, y, ...); panel.loess(x, y, ...)
}
)
sps <- trellis.par.get("superpose.symbol")
sps$pch <- 1:7
trellis.par.set("superpose.symbol", sps)
xyplot(Time ~ Viscosity, data = stormer, groups = Wt,
panel = panel.superpose, type = "b",
key = list(columns = 3,
text = list(paste(c("Weight: ", "", ""),
unique(stormer$Wt), "gms")),
points = Rows(sps, 1:3)
)
)
rm(sps)
topo.plt <- expand.grid(topo.mar)
topo.plt$pred <- as.vector(predict(topo.loess, topo.plt))
levelplot(pred ~ x * y, topo.plt, aspect = 1,
at = seq(690, 960, 10), xlab = "", ylab = "",
panel = function(x, y, subscripts, ...) {
panel.levelplot(x, y, subscripts, ...)
panel.xyplot(topo$x,topo$y, cex = 0.5, col = 1)
}
)
wireframe(pred ~ x * y, topo.plt, aspect = c(1, 0.5),
drape = TRUE, screen = list(z = -150, x = -60),
colorkey = list(space="right", height=0.6))
lcrabs.pc <- predict(princomp(log(crabs[,4:8])))
crabs.grp <- c("B", "b", "O", "o")[rep(1:4, each = 50)]
splom(~ lcrabs.pc[, 1:3], groups = crabs.grp,
panel = panel.superpose,
key = list(text = list(c("Blue male", "Blue female",
"Orange Male", "Orange female")),
points = Rows(trellis.par.get("superpose.symbol"), 1:4),
columns = 4)
)
sex <- crabs$sex
levels(sex) <- c("Female", "Male")
sp <- crabs$sp
levels(sp) <- c("Blue", "Orange")
splom(~ lcrabs.pc[, 1:3] | sp*sex, cex = 0.5, pscales = 0)
Quine <- quine
levels(Quine$Eth) <- c("Aboriginal", "Non-aboriginal")
levels(Quine$Sex) <- c("Female", "Male")
levels(Quine$Age) <- c("primary", "first form",
"second form", "third form")
levels(Quine$Lrn) <- c("Average learner", "Slow learner")
bwplot(Age ~ Days | Sex*Lrn*Eth, data = Quine)
bwplot(Age ~ Days | Sex*Lrn*Eth, data = Quine, layout = c(4, 2),
strip = function(...) strip.default(..., style = 1))
stripplot(Age ~ Days | Sex*Lrn*Eth, data = Quine,
jitter = TRUE, layout = c(4, 2))
stripplot(Age ~ Days | Eth*Sex, data = Quine,
groups = Lrn, jitter = TRUE,
panel = function(x, y, subscripts, jitter.data = FALSE, ...) {
if(jitter.data) y <- jitter(as.numeric(y))
panel.superpose(x, y, subscripts, ...)
},
xlab = "Days of absence",
between = list(y = 1), par.strip.text = list(cex = 0.7),
key = list(columns = 2, text = list(levels(Quine$Lrn)),
points = Rows(trellis.par.get("superpose.symbol"), 1:2)
),
strip = function(...)
strip.default(..., strip.names = c(TRUE, TRUE), style = 1)
)
fgl0 <- fgl[ ,-10] # omit type.
fgl.df <- data.frame(type = rep(fgl$type, 9),
y = as.vector(as.matrix(fgl0)),
meas = factor(rep(1:9, each = 214), labels = names(fgl0)))
stripplot(type ~ y | meas, data = fgl.df,
scales = list(x = "free"), xlab = "", cex = 0.5,
strip = function(...) strip.default(style = 1, ...))
if(FALSE) { # no data supplied
xyplot(ratio ~ scant | subject, data = A5,
xlab = "scan interval (years)",
ylab = "ventricle/brain volume normalized to 1 at start",
subscripts = TRUE, ID = A5$ID,
strip = function(factor, ...)
strip.default(..., factor.levels = labs, style = 1),
layout = c(8, 5, 1),
skip = c(rep(FALSE, 37), rep(TRUE, 1), rep(FALSE, 1)),
panel = function(x, y, subscripts, ID) {
panel.xyplot(x, y, type = "b", cex = 0.5)
which <- unique(ID[subscripts])
panel.xyplot(c(0, 1.5), pr3[names(pr3) == which],
type = "l", lty = 3)
if(which == 303 || which == 341) points(1.4, 1.3)
})
}
Cath <- equal.count(swiss$Catholic, number = 6, overlap = 0.25)
xyplot(Fertility ~ Education | Cath, data = swiss,
span = 1, layout = c(6, 1), aspect = 1,
panel = function(x, y, span) {
panel.xyplot(x, y); panel.loess(x, y, span)
}
)
Cath2 <- equal.count(swiss$Catholic, number = 2, overlap = 0)
Agr <- equal.count(swiss$Agric, number = 3, overlap = 0.25)
xyplot(Fertility ~ Education | Agr * Cath2, data = swiss,
span = 1, aspect = "xy",
panel = function(x, y, span) {
panel.xyplot(x, y); panel.loess(x, y, span)
}
)
Cath
levels(Cath)
plot(Cath, aspect = 0.3)
# End of ch04
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.