Splus Textbook Examples: Visualizing Data, Chapter 4



Stat Computing > Splus > Textbook Examples > Visualizing Data

Splus Textbook Examples
Visualizing Data: Chapter 4, Trivariate Data

Figure 4.1

splom(~rubber[,1:3], 
  varnames = c("Hardness\n(degrees Shore)",
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"),
  sub = list("Figure 4.1",cex=.8))

Figure 4.2

splom(~rubber[, 1:3], 
  panel = function(x, y) {
       i <- rubber[,1] < 62
       panel.splom(x[!i], y[!i])
       panel.splom(x[i], y[i], pch = "+")
  },
  sub = list("Figure 4.2",cex=.8),
  varnames = c("Hardness\n(degrees Shore)", 
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"))

Figure 4.3

attach(rubber)
Hardness <- equal.count(hardness,6,3/4)
ans <- xyplot(abrasion.loss ~ tensile.strength | Hardness,
  prepanel=function(x,y) prepanel.loess(x, y, span = 3/4, degree = 1),
  panel = function(x, y) {
    panel.grid(h=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 3/4, degree = 1)
  },
  layout = c(3,2),
  xlab = "Tensile Strength (kg/square cm)",
  ylab = "Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber$hardness,6,3/4), 
  aspect=.2,
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  sub = list("Figure 4.3",cex=.8),
  xlab = "Hardness (degrees Shore)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.4

attach(rubber)
Tensile.strength <- equal.count(tensile.strength,6,3/4)
ans <- xyplot(abrasion.loss ~ hardness | Tensile.strength,
  prepanel=function(x,y) prepanel.loess(x, y, span = 3/4, degree = 1),
  panel = function(x, y){
    panel.grid()
      panel.xyplot(x,y)
      panel.loess(x, y, span = 3/4, degree = 1)
  },
  layout = c(3, 2),
  xlab = "Hardness (degrees Shore)",
  ylab = "Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber$tensile.strength,6,3/4), 
  aspect=.25,
  sub = list("Figure 4.4",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Tensile Strength (kg/square cm)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.5

splom(~ethanol, 
  sub = list("Figure 4.5",cex=.8),
  aspect = 1,
  varnames=c("NOx\n(micrograms/J)",
    "Compression\nRatio",
    "Equivalence\nRatio"))

Figure 4.6

attach(ethanol)
Equivalence.Ratio <- equal.count(E, number = 9, overlap = 0.25)
ans <- xyplot(NOx ~ C | Equivalence.Ratio,
  prepanel= function(x,y) prepanel.loess(x, y, span = 1),
  panel = function(x, y) {
    panel.grid(v = 2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1)
  },
  aspect = 2.5, # banking to 45 degrees results in aspect that is too big 
  layout = c(5, 2),
  xlab = "Compression Ratio",
  ylab = "NOx (micrograms/J)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(ethanol$E, number = 9, overlap = 0.25), 
  sub = list("Figure 4.6",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Equivalence Ratio"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.7

attach(ethanol)
Compression.Ratio <- C
ans <- xyplot(NOx ~ E | Compression.Ratio,
  prepanel=function(x,y) prepanel.loess(x, y, span = 3/4, degree = 2),
  panel = function(x, y) {
    panel.grid(h=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 3/4, degree = 2)
  },
  layout=c(3,2),
  xlab = "Equivalence Ratio",
  ylab = "NOx (micrograms/J)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(ethanol$C), 
  sub = list("Figure 4.7",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Compression Ratio"),
  position=c(0,.2,1,.5))
invisible()

Figure 4.8

splom(~rubber[, 1:3], 
  sub = list("Figure 4.8",cex=.8),
  varnames = c("Hardness\n(degrees Shore)", 
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"))

Figure 4.9

splom(~rubber[, 1:3], 
  panel = function(x, y) {
       i <- rubber[,1] < 62
       panel.splom(x[!i], y[!i])
       panel.splom(x[i], y[i], pch = "+")
  },
  sub = list("Figure 4.9",cex=.8),
  varnames = c("Hardness\n(degrees Shore)", 
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"))

Figure 4.10

splom(~rubber[, 1:3], 
  panel = function(x, y) {
       i <- rubber[,1] > 54 & rubber[,1] < 72
       panel.splom(x[!i], y[!i])
       panel.splom(x[i], y[i], pch = "+")
  },
  sub = list("Figure 4.10",cex=.8),
  varnames = c("Hardness\n(degrees Shore)", 
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"))

Figure 4.11

splom(~rubber[, 1:3], 
  panel = function(x, y) {
       i <- rubber[,2] < 170
       panel.splom(x[!i], y[!i])
       panel.splom(x[i], y[i], pch = "+")
  },
  sub = list("Figure 4.11",cex=.8),
  varnames = c("Hardness\n(degrees Shore)", 
    "Tensile Strength\n(kg/square cm)",
    "Abrasion Loss\n(g/hp-hour)"))

Figure 4.12

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C",family="s")
CE.marginal <- list(C=seq(min(C),max(C),length=2),
  E=seq(min(E),max(E),length=16))
CE.grid <- expand.grid(CE.marginal)
eth.fit <- predict(eth.m,CE.grid)
EQ.Ratio <- CE.grid$E
ans <- xyplot(eth.fit ~ CE.grid$C | EQ.Ratio,
  aspect=2,
  layout = c(8, 2),
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab = "Compression Ratio",
  ylab = "NOx (micrograms/J)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

ans <- plot(shingle(seq(min(ethanol$E), max(ethanol$E), length=16)), 
  sub = list("Figure 4.12",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Equivalence Ratio")
print(ans, position=c(0,.1,1,.45))
invisible()

Figure 4.13

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C",family="s")
C.marginal <- seq(min(C),max(C),length=16)
E.marginal <- seq(min(E),max(E),length=50)
CE.marginal <- list(C=C.marginal,E=E.marginal)
CE.grid <- expand.grid(CE.marginal)
eth.fit <- predict(eth.m,CE.grid)
Compression.Ratio <- CE.grid$C
ans <- xyplot(eth.fit ~ CE.grid$E | Compression.Ratio,
  panel = function(x, y) {
    panel.grid(h=2)
    panel.xyplot(x, y, type = "l")
  },
  aspect="xy",
  layout = c(4, 4),
  xlab = "Equivalence Ratio",
  ylab = "NOx (micrograms/J)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

ans <- plot(shingle(seq(min(ethanol$C), max(ethanol$C), length=16)), 
  sub = list("Figure 4.13",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Compression Ratio")
print(ans, position=c(0,.1,1,.45))
invisible()

Figure 4.14

attach(ethanol)
c.values <- seq(min(C), max(C), length = 16)
e.values <- seq(min(E), max(E), length = 16)
ans <- xyplot(E ~ C,
  panel = substitute(function(x, y){
    segments(rep(min(x),16),e.values,rep(max(x),16),e.values)
    segments(c.values,rep(min(y),16),c.values,rep(max(y),16))
    panel.xyplot(x, y, col=0, pch=16)  # cover grid lines
    panel.xyplot(x, y)
  }),
  aspect = 1,
  sub = list("Figure 4.14",cex=.8),
  xlab = "Compression Ratio",
  ylab = "Equivalence Ratio")
detach()
ans

Figure 4.14

xyplot(hardness ~ tensile.strength,
  data = rubber,
  aspect=1,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.abline(h=c(55,83),v=c(144,237))
  },
  ylab= "Hardness (degrees Shore)",
  sub = list("Figure 4.15",cex=.8),
  xlab= "Tensile Strength (kg/square cm)")

Figure 4.16

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
grid <- expand.grid(hardness = seq(55, 83, length = 4),
  tensile.strength = c(seq(144,180,length=50),c(181, 190)))
center <- grid$tensile.strength-180
grid$ts.low <- center*(center<0)
rub.fit <- predict(rub.lm,grid)
Hardness <- grid$hardness
ans <- xyplot(rub.fit ~ grid$tensile.strength | Hardness,
  aspect="xy",
  layout = c(4, 1),
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab = "Tensile Strength (kg/square cm)", 
  ylab = "Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

print(plot(shingle(seq(55, 83, length = 4)), 
  sub = list("Figure 4.16",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Hardness (degrees Shore)"), 
  position=c(0,.1,1,.45))
invisible()

Figure 4.17

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
grid<-expand.grid(tensile.strength = seq(144,180,length = 3), 
  hardness = seq(55, 83,length=50))
#this mimics the computation of ts.low for this cropped version of tensile strength
center <- grid$tensile.strength-180
grid$ts.low <- center*(center<0)
rub.fit <- predict(rub.lm,grid)
Tensile.strength <- grid$tensile.strength
ans <- xyplot(rub.fit ~ grid$hardness | Tensile.strength,
  aspect="xy",
  layout = c(3, 1), 
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab = "Hardness (degrees Shore)",
  ylab = "Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

print(plot(shingle(seq(144,180,length = 3)), 
  sub = list("Figure 4.17",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Tensile Strength (kg/square cm)"),
  position=c(0,.1,1,.45))
invisible()

Figure 4.18

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
ans <- xyplot(residuals(rub.lm) ~ tensile.strength,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=1)
    panel.abline(h=0)
  },
  aspect=1,
  sub = list("Figure 4.18",cex=.8),
  xlab = "Tensile Strength (kg/square cm)",
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
ans

Figure 4.19

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
ans <- xyplot(residuals(rub.lm) ~ hardness,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=1)
    panel.abline(h=0)
  },
  aspect=1,
  sub = list("Figure 4.19",cex=.8),
  xlab = "Hardness (degrees Shore)",
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
ans

Figure 4.20

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
Hardness <- equal.count(hardness,6,3/4)
ans <- xyplot(residuals(rub.lm) ~ tensile.strength | Hardness,
  panel = function(x, y, ...) {
    panel.grid(v=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1, degree = 1)
    panel.abline(h=0)
  },
  aspect = 2,
  layout = c(3, 2),
  xlab = "Tensile Strength (kg/square cm)", 
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber$hardness,6,3/4), 
  sub = list("Figure 4.20",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Hardness (degrees Shore)"),    position=c(0,.1,1,.425))
invisible()

Figure 4.21

attach(rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness+ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
Tensile.Strength <- equal.count(tensile.strength,6,3/4)
ans <- xyplot(residuals(rub.lm) ~ hardness | Tensile.Strength,
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1, degree = 1)
    panel.abline(h=0)
  },
  aspect = 2,
  layout = c(3, 2),
  xlab = "Hardness (degrees Shore)", 
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber$tensile.strength,6,3/4), 
  sub = list("Figure 4.21",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Tensile Strength (kg/square cm)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.22

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low, weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
Hardness <- equal.count(hardness,6,3/4)
ans <- xyplot(residuals(rub.lm) ~ tensile.strength | Hardness,
  panel = function(x, y, ...) {
    panel.grid(v=2, h=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1, degree = 1)
    panel.abline(h=0)
  },
  aspect = 1,
  layout = c(3, 2),
  xlab = "Tensile Strength (kg/square cm)",
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber[rubber$tensile.strength>130,]$hardness,6,3/4),
  sub = list("Figure 4.22",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Hardness (degrees Shore)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.23

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
Tensile.Strength <- equal.count(tensile.strength,6,3/4)
ans <- xyplot(residuals(rub.lm) ~ hardness | Tensile.Strength,
  panel = function(x, y, ...) {
    panel.grid(v=2, h=2)
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1, degree = 1)
    panel.abline(h=0)
  },
  aspect = 1,
  layout = c(3, 2),
  xlab = "Hardness (degrees Shore)",
  ylab = "Residual Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(rubber[rubber$tensile.strength>130,]$tensile.strength,6,3/4), 
  sub = list("Figure 4.23",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Tensile Strength (kg/square cm)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.24

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low,weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
grid <- expand.grid(hardness = seq(55, 83, length = 4),
  tensile.strength = c(seq(144,180,length=50),181, 190))
center <- grid$tensile.strength-180
grid$ts.low <- center*(center<0)
rub.fit <- predict(rub.lm,grid)
Hardness <- grid$hardness
ans <- xyplot(rub.fit ~ grid$tensile.strength | Hardness,
  aspect="xy",
  layout = c(4, 1),
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab = "Tensile Strength (kg/square cm)",
  ylab = "Abrasion Loss (g/hp-hour)")
detach()  print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(seq(55, 83, length = 4)),     
  sub = list("Figure 4.24",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Hardness (degrees Shore)"),
  position=c(0,.1,1,.425))
invisible()

Figure

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low, weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
grid<-expand.grid(tensile.strength = seq(144,180,length = 3),   
    hardness = seq(55, 83,length=50))
center <- grid$tensile.strength-180
grid$ts.low <- center*(center<0)
rub.fit <- predict(rub.lm,grid)
Tensile.Strength <- grid$tensile.strength
ans <- xyplot(rub.fit ~ grid$hardness | Tensile.Strength,    aspect="xy",
  layout = c(3, 1),
  panel = function(x, y) {
    panel.grid(h=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab = "Hardness (degrees Shore)",
  ylab = "Abrasion Loss (g/hp-hour)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(seq(144,180,length = 3)), 
  sub = list("Figure 4.25",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Tensile Strength (kg/square cm)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.26

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low, weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
ans <- xyplot(sqrt(abs(rub.lm$res)) ~ rub.lm$fit,
  panel = function(x, y){
    panel.xyplot(x, y)
    panel.loess(x, y, span=2)
  },
  aspect=1,
  sub = list("Figure 4.26",cex=.8),
  xlab = "Fitted Abrasion Loss (g/hp-hour)",
  ylab = "Square Root Absolute Residual Abrasion Loss (g/hp-hour)")
detach()
ans

Figure 4.27

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low, weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
ans <- qqmath(~residuals(rub.lm), 
  prepanel = prepanel.qqmathline,
  panel = function(x, y) {
    panel.qqmathline(y, distribution = qnorm)
    panel.qqmath(x, y)
  },
  aspect=1,
  sub = list("Figure 4.27",cex=.8),
  xlab = "Unit Normal Quantile",
  ylab="Residual Abrasion Loss (g/hp-hour)")
detach()
ans

Figure 4.28

new.rubber <- rubber[rubber$tensile.strength>130,]
attach(new.rubber)
wt <- rep(1,length(abrasion.loss))
for(i in 1:10){
  rub.lm <- lm(abrasion.loss~hardness*ts.low, weights=wt)
  wt <- wt.bisquare(rub.lm$res/median(abs(rub.lm$res)),c=6)
ans <- rfs(rub.lm,
  sub = list("Figure 4.28",cex=.8),
  aspect=1.5,
  ylab="Residual Abrasion Loss (g/hp-hour)")
detach()
ans

Figure 4.29

xyplot(loess(NOx ~ C * E, span = 1/2, degree = 2,
    parametric = "C", drop.square = "C",family="s")$residuals ~ E,
  data = ethanol,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=1/2, degree=1)
    panel.abline(h=0)
  },
  aspect=1,
  sub = list("Figure 4.29",cex=.8),
  xlab="Equivalence Ratio",
  ylab = "Residual NOx (micrograms/J)")

Figure 4.30

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C",family="s")
ans <- xyplot(sqrt(abs(residuals(eth.m))) ~ fitted.values(eth.m),
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=2)
  },
  aspect=1,
  scales=list(cex=.9),
  sub = list("Figure 4.30",cex=.8),
  xlab = list("Fitted NOx (micrograms/J)",cex=.9),
  ylab = list("Square Root Absolute Residual NOx (square root micrograms/J)",cex=.9))
detach()
ans

Figure 4.31

qqmath(~loess(NOx ~ C * E, data = ethanol, span = 1/3, degree = 2,
    parametric = "C", drop.square = "C",family="s")$residuals,
  prepanel = prepanel.qqmathline,
  panel = function(x, y) {
    panel.qqmathline(y, distribution = qnorm)
    panel.qqmath(x, y)
  },
  aspect=1, 
  sub = list("Figure 4.31",cex=.8),
  xlab = "Unit Normal Quantile",
  ylab = "Residual NOx (micrograms/J)")

Figure 4.32

rfs(loess(NOx ~ C * E, data = ethanol, span = 1/3, degree = 2,
    parametric = "C", drop.square = "C",family="s"),
  sub = list("Figure 4.32",cex=.8),
  aspect=2,
  ylab = "NOx (micrograms/J)")

Figure 4.33

attach(ethanol)
eth.m <- loess(logb(NOx,2) ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C",family="s")
ans <- xyplot(sqrt(abs(eth.m$res)) ~ eth.m$fit,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=2)
  },
  aspect=1,
  scales=list(cex=.8),
  sub = list("Figure 4.33",cex=.8),
  xlab = list("Fitted Log NOx (log 2 micrograms/J)",cex=.8),
  ylab = list("Square Root Absolute Residual Log NOx\n(square root absolute log 2 micrograms/J)\n",cex=.8))
detach()
ans

# no fig 4.34

Figure 4.35

set.seed(19)
new.slope <- slope
new.slope$percent <- jitter(new.slope$percent,2)
splom(~new.slope,
  varnames =c("Absolute\nError\n(%)",
    "Jittered\nPercent\n(%)",
    "Distance\n(degrees)",
    "Resolution\n(degrees)"),
sub = list("Figure 4.35",cex=.8))

Figure 4.36

print(xyplot(error ~ distance | percent,
  data = slope, 
  aspect="xy",
  layout = c(6,2),
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y)
  },
  xlab="Distance (degrees)",
  ylab="Absolute Error (%)"),
  position=c(0,.375,1,1), more=T)

print(plot(shingle(slope$percent), 
  sub = list("Figure 4.36",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Percent (%)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.37

print(xyplot(error ~ resolution | percent,
  data = slope,
  aspect="xy",
  layout=c(6,2),
  panel = function(x, y) {
    panel.grid(h=2, v=2)
    panel.xyplot(x, y)
  },
  xlab="Resolution (degrees)",
  ylab="Absolute Error (%)"),
  position=c(0,.375,1,1), more=T)

print(plot(shingle(slope$percent),
  sub = list("Figure 4.37",cex=.8),  
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Percent (%)"),
  position=c(0,.25,1,.55))
invisible()

Figure 4.38

attach(slope)
Resolution <- equal.count(resolution,8,1/4)
ans <- xyplot(error ~ percent | Resolution,
  aspect=1,
  layout=c(4,2),
  panel = function(x, y) {
    panel.grid()
    panel.xyplot(x, y)
  },
  xlab="Percent (%)",
  ylab="Absolute Error (%)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

print(plot(equal.count(slope$resolution,8,1/4), 
  sub = list("Figure 4.38",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Resolution (degrees)"),
  position=c(0,.1,1,.45))
invisible()

Figure 4.39

attach(slope)
wt <- rep(1,length(error))
for(i in 1:10){
  slo.lm <- lm(error~percent+resolution,weights=wt)
  wt <- wt.bisquare(slo.lm$res/median(abs(slo.lm$res)),c=6)
ans <- xyplot(residuals(slo.lm) ~ percent,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=1/2)
    panel.abline(h=0)
  },
  aspect=1,
  sub = list("Figure 4.39",cex=.8),
  xlab="Percent (%)",
  ylab="Residual Absolute Error (%)")
detach()
ans

Figure 4.40

attach(slope)
wt <- rep(1,length(error))
for(i in 1:10){
  slo.lm <- lm(error~percent+resolution,weights=wt)
  wt <- wt.bisquare(slo.lm$res/median(abs(slo.lm$res)),c=6)
ans <- xyplot(residuals(slo.lm) ~ resolution,
  panel = function(x, y) {
    panel.xyplot(x, y)
    panel.loess(x, y, span=1/2)
    panel.abline(h=0)
  },
  aspect=1,
  sub = list("Figure 4.40",cex=.8),
  xlab="Resolution (degrees)",
  ylab="Residual Absolute Error (%)")
detach()
ans

Figure 4.41

attach(slope)
wt <- rep(1,length(error))  for(i in 1:10){
  slo.lm <- lm(error~percent+resolution,weights=wt)
  wt <- wt.bisquare(slo.lm$res/median(abs(slo.lm$res)),c=6)
ans <- rfs(slo.lm, 
  sub = list("Figure 4.41",cex=.8),
  aspect = 1.5,    ylab="Residual Absolute Error (%)")
detach()
ans

Figure 4.42

attach(slope)
grid <- expand.grid(d = seq(0, 45, length = 46), 
        p = seq(min(percent), max(percent), length = 11))
p <- grid$p
d <- grid$d
a <- 52.7 - 1.19 * asin(cos(d*pi/90) * (100 - p) / (100 + p)) * 180/pi - 0.48 * p
Percent <- p
ans <- xyplot(a ~ d | Percent,
  aspect="xy",
  layout=c(4,3),
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y, type = "l")
  },
  xlab="Distance (degrees)",
  ylab="Absolute Error (%)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(seq(min(slope$percent), 
  max(slope$percent), length = 11)), 
  sub = list("Figure 4.42",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Percent (%)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.43

attach(galaxy)
set.seed(19)
ans <- xyplot(jitter(north.south, .3) ~ jitter(east.west, .3),
  aspect = diff(range(north.south))/diff(range(east.west)),
  sub = list("Figure 4.43",cex=.8),
  xlab = "Jittered East-West Coordinate (arcsec)",
  ylab = "Jittered South-North Coordinate (arcsec)")
detach()
ans

Figure 4.44

attach(galaxy)
set.seed(19)
Velocity <- equal.count(velocity,15,1/4)
ans <- xyplot(jitter(north.south, .3) ~ jitter(east.west, .3) | Velocity,
  layout = c(5, 3), 
  panel = function(x, y) {
    panel.grid(v=2)
    panel.xyplot(x, y)
  },
  aspect=diff(range(north.south))/diff(range(east.west)),
  xlab = "Jittered East-West Coordinate (arcsec)",
  ylab = "Jittered South-North Coordinate (arcsec)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(equal.count(galaxy$velocity,15,1/4), 
  sub = list("Figure 4.44",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Velocity (km/sec)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.45

ans <- xyplot(velocity ~ radial.position | angle, 
  data = galaxy,
  prepanel=function(x,y) prepanel.loess(x, y, span = 1/2, degree = 2),
  panel = function(x, y) {
    panel.grid()
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1/2, degree = 2)
  },
  layout = c(4, 2),
  xlab = "Radial Position (arcsec)",
  ylab = "Velocity (km/sec)")
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(galaxy$angle), 
  sub = list("Figure 4.45",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Slit (degrees)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.46

attach(galaxy)
gal.m <- loess(velocity ~ east.west * north.south, 
         span = 0.25, degree = 2, normalize = F, family = "symmetric")
slitfit <- fitted.values(gal.m)
ans <- xyplot(velocity ~ radial.position | angle,
  prepanel = substitute(function(x, y, subscripts){
    k <- order(x)
    list(dx = diff(x[k]), dy = diff(slitfit[subscripts][k]))
  }),
  panel = substitute(function(x, y, subscripts) {
    add.line <- trellis.par.get("add.line")
    panel.grid()
    panel.xyplot(x,y)
     k <- order(x)
    lines(x[k], slitfit[subscripts][k],
      lwd = add.line$lwd, lty = add.line$lty, col = add.line$col)
  }),
  layout = c(4, 2),
  xlab = "Radial Position (arcsec)",
  ylab = "Velocity (km/sec)")
detach()
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(galaxy$angle), 
  sub = list("Figure 4.46",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Slit (degrees)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.47

ans <- xyplot(loess(velocity ~ east.west * north.south, 
  span = 0.25,
  degree = 2, normalize = F, 
  family = "symmetric")$residuals ~ radial.position | angle,  
  data = galaxy,
  panel = function(x, y) {
    panel.grid()
    panel.xyplot(x, y)
    panel.loess(x, y, span = 1, degree = 2)
    panel.abline(h=0)
  },
  aspect=1,
  layout = c(4, 2),
  xlab = "Radial Position (arcsec)",
  ylab = "Residual Velocity (km/sec)")
print(ans, position=c(0,.375,1,1), more=T)

print(plot(shingle(galaxy$angle), 
  sub = list("Figure 4.47",cex=.8),
  scales=list(x=list(cex=.7),y=list(cex=.7)),
  xlab = "Slit (degrees)"),
  position=c(0,.1,1,.425))
invisible()

Figure 4.48

qqmath(~loess(velocity ~ east.west * north.south, data = galaxy, 
    span = 0.25, degree = 2, normalize = F, 
    family = "symmetric")$residuals,
  prepanel = prepanel.qqmathline,
  panel = function(x, y) {
    panel.qqmathline(y, distribution = qnorm)
    panel.qqmath(x, y)
  },
  aspect=1,
  sub = list("Figure 4.48",cex=.8),
  xlab="Unit Normal Quantile",
  ylab="Residual Velocity (km/sec)")

Figure 4.49

rfs(loess(velocity ~ east.west * north.south, data = galaxy,
    span = 0.25, degree = 2, normalize = F, 
    family = "symmetric"),
  sub = list("Figure 4.49",cex=.8),
  aspect=1.5,
  ylab="Velocity (km/sec)")

Figure 4.50

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25,length=101),
  north.south = seq(-45, 45,length=181))
galaxy.grid <- expand.grid(galaxy.marginal)
gal.m <- loess(velocity ~ east.west * north.south, span = 0.25, 
  degree = 2, normalize = F, family = "symmetric")
galaxy.fit <- predict(gal.m, galaxy.grid)
x <- range(galaxy.marginal$east.west)
y <- range(galaxy.marginal$north.south)
ans <- contourplot(galaxy.fit ~ galaxy.grid$east.west * galaxy.grid$north.south,
  at = seq(1435, 1755, by = 5),
  panel = function(..., at) {
    ref <- trellis.par.get("reference.line")
    panel.contourplot(..., labels = F,
      at = at,
      lty = ref$lty, lwd = ref$lwd, col = ref$col)
    panel.contourplot(..., labels = F,
      at = seq(1460, 1740, by = 40))
    panel.contourplot(...,
      at = seq(1440, 1760, by = 40))
  },  
  aspect = diff(range(y))/diff(range(x)),
  sub = list("Figure 4.50",cex=.8),
  xlab = "East-West Coordinate (arcsec)",
  ylab = "South-North Coordinate (arcsec)")
detach()
ans

# no figs for 4.51 4.54

Figure 4.55

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25,length=101),
  north.south = seq(-45, 45,length=181))
galaxy.grid <- expand.grid(galaxy.marginal)
gal.m <- loess(velocity ~ east.west * north.south, span = 0.25, 
  degree = 2, normalize = F, family = "symmetric")
galaxy.fit <- predict(gal.m, galaxy.grid)
n.level <- 48
given <- seq(min(galaxy.fit),max(galaxy.fit),length=n.level+1)
Velocity <- shingle(as.vector(galaxy.fit), 
  cbind(given[-(n.level+1)],given[-1]))
ans <- xyplot(galaxy.grid$n ~ galaxy.grid$e | Velocity,
  layout = c(8, 6),
  panel = function(x, y) {
    panel.grid(v=2, h=2)
    panel.xyplot(x, y, pch = ".")
  },
  aspect=diff(range(galaxy.grid$n))/diff(range(galaxy.grid$e)),
  xlab = "East-West Coordinate (arcsec)",
  ylab = "South-North Coordinate (arcsec)")
detach()
print(ans, position=c(0,.3,1,1), more=T)

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25,length=101),
  north.south = seq(-45, 45,length=181))
galaxy.grid <- expand.grid(galaxy.marginal)
gal.m <- loess(velocity ~ east.west * north.south, span = 0.25, 
  degree = 2, normalize = F, family = "symmetric")
galaxy.fit <- predict(gal.m, galaxy.grid)
given <- seq(min(galaxy.fit),max(galaxy.fit),length=49)
ans <- plot(shingle(as.vector(galaxy.fit), 
  intervals = cbind(given[-length(given)], given[-1])),
  cex=.7,
  scales=list(x=list(cex=.7),y=list(at=seq(10,40,10),cex=.7)),
  sub = list("Figure 4.55",cex=.8),
  ylab="",
  xlab = "Velocity (km/sec)")
detach()
print(ans, position=c(0,.05,1,.35))
invisible()

Figure 4.56

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25, by = 2), 
  north.south = seq(-45, 45, by = 2))
galaxy.grid <- expand.grid(galaxy.marginal)
galaxy.fit <- predict(loess(velocity ~ east.west * north.south, 
  span = 0.25, degree = 2, normalize = F, family = "symmetric"),
  galaxy.grid)
ar <- diff(range(galaxy.grid$north.south))/diff(range(galaxy.grid$east.west))
ans <- wireframe(galaxy.fit ~ galaxy.grid$east.west * galaxy.grid$north.south,
  screen = list(z = 30, x = -60, y = 0),
  aspect = c(ar, 1.3),
  sub = list("Figure 4.56",cex=.8),
  xlab = "East-West",
  ylab = "South-North",
  zlab = "Velocity")
detach()
ans

Figure 4.57

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25, by = 2), 
  north.south = seq(-45, 45, by = 2))
galaxy.grid <- expand.grid(galaxy.marginal)
galaxy.fit <- predict(loess(velocity ~ east.west * north.south, 
  span = 0.25, degree = 2, normalize = F, family = "symmetric"),
  galaxy.grid)
ar <- diff(range(galaxy.grid$north.south))/diff(range(galaxy.grid$east.west))
ans <- wireframe(galaxy.fit ~ galaxy.grid$east.west * galaxy.grid$north.south,
  screen = list(z = 120, x = -60, y = 0),
  aspect = c(ar, 1.3),
  sub = list("Figure 4.57",cex=.8),
  xlab = "East-West",
  ylab = "South-North",
  zlab = "Velocity")
detach()
ans

Figure 4.58

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25, by = 2), 
  north.south = seq(-45, 45, by = 2))
galaxy.grid <- expand.grid(galaxy.marginal)
galaxy.fit <- predict(loess(velocity ~ east.west * north.south, 
  span = 0.25, degree = 2, normalize = F, family = "symmetric"),
  galaxy.grid)
ar <- diff(range(galaxy.grid$north.south))/diff(range(galaxy.grid$east.west))
ans <- wireframe(galaxy.fit ~ galaxy.grid$east.west * galaxy.grid$north.south,
  screen = list(z = 210, x = -60, y = 0),
  aspect = c(ar, 1.3),
  sub = list("Figure 4.58",cex=.8),
  xlab = "East-West",
  ylab = "South-North",
  zlab = "Velocity")
detach()
ans

Figure 4.59

attach(galaxy)
galaxy.marginal <- list(east.west = seq(-25, 25, by = 2), 
  north.south = seq(-45, 45, by = 2))
galaxy.grid <- expand.grid(galaxy.marginal)
galaxy.fit <- predict(loess(velocity ~ east.west * north.south, 
  span = 0.25, degree = 2, normalize = F, family = "symmetric"),
  galaxy.grid)
ar <- diff(range(galaxy.grid$north.south))/diff(range(galaxy.grid$east.west))
ans <- wireframe(galaxy.fit ~ galaxy.grid$east.west * galaxy.grid$north.south,
  screen = list(z = 300, x = -60, y = 0),
  aspect = c(ar, 1.3),
  sub = list("Figure 4.59",cex=.8),
  xlab = "East-West",
  ylab = "South-North",
  zlab = "Velocity")
detach()
ans

Figure 4.60

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C", family = "s")
eth.marginal <- list(C = seq(min(C), max(C), length = 6), 
  E = seq(min(E), max(E), length = 25))
eth.grid <- expand.grid(eth.marginal)
eth.fit <- predict(eth.m, eth.grid)
ans <- wireframe(eth.fit ~ eth.grid$C * eth.grid$E, 
  screen = list(z = 94, x = -95, y = 0),
  distance = .3,
  sub = list("Figure 4.60",cex=.8),
  xlab = "C",    ylab = "E",
  zlab = "NOx")
detach()
ans

Figure 4.61

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2, parametric = "C", 
  drop.square = "C", family="s")
eth.marginal <- list(C = seq(min(C), max(C), length = 25), 
  E = seq(min(E), max(E), length = 25))
eth.grid <- expand.grid(eth.marginal)
eth.fit <- predict(eth.m, eth.grid)
ans <- wireframe(eth.fit ~ eth.grid$C * eth.grid$E, 
  screen = list(z = 30, x = -60, y=0),
  distance = .3,
  sub = list("Figure 4.61",cex=.8),
  xlab = "C",
  ylab = "E",
  zlab = "NOx")
detach()
ans

Figure 4.62

attach(ethanol)
eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
  parametric = "C", drop.square = "C",family="s")
eth.marginal <- list(C = seq(min(C), max(C), length = 6), 
  E = seq(min(E), max(E), length = 25))
eth.grid <- expand.grid(eth.marginal)
eth.fit <- predict(eth.m, eth.grid)
ans <- wireframe(eth.fit ~ eth.grid$C * eth.grid$E, 
  screen = list(z = 90, x = -90, y = 0),
  perspective = F,
  distance =.3,
  sub = list("Figure 4.62",cex=.8),
  xlab = "",
  ylab = "E",
  zlab = "NOx")
detach()
ans

Figure 4.63

ar <- diff(range(galaxy$north.south))/diff(range(galaxy$east.west))
left <- cloud(velocity ~ east.west*north.south,
  data = galaxy,
  screen=list(z=123, x=-60, y=0),
  perspective=F,
  aspect = c(ar, 1.6),
  sub = list("Figure 4.63",cex=.8),
  scales=list(cex=.7),
  xlab=list("East-West",cex=.4),
  ylab=list("South-North",cex=.4),
  zlab=list("Velocity",cex=.4))
right <- update(left, screen=list(z=117,x =-60, y=0), sub="")
print(left, split = c(1,1,2,1), more = T)
print(right, split = c(2,1,2,1))
invisible()

Figure 4.64

xyplot(northing ~ easting, 
  data = soil,
  panel = function(x, y) points(x, y, pch=".", cex=.75),
  aspect=diff(range(soil$northing))/diff(range(soil$easting)),
  ylab= "Northing (km)",
  sub = list("Figure 4.64",cex=.8),
  xlab= "Easting (km)")

Figure 4.65

xyplot(northing ~ resistivity | track,
  data = soil,
  subset = is.ns,
  layout=c(4,2),
  panel=function(x, y) {
    panel.grid(v=2)
    points(x, y, pch=".", cex=1.25)
  },
  sub = list("Figure 4.65",cex=.8),
  xlab="Resistivity (ohm-cm)",
  ylab="Northing (km)")

Figure 4.66

xyplot(resistivity[!is.ns] ~ easting[!is.ns] | track[!is.ns],
  data = soil,
  layout=c(4,10),
  panel=function(x, y) {
    panel.grid(h=2)
    points(x, y, pch=".", cex=1.25)
  },
  sub = list("Figure 4.66",cex=.8),
  xlab="Easting (km)",
  ylab="Resistivity (ohm-cm)")

Figure 4.67

attach(soil)
soi.marginal <- list(easting = seq(.15, 1.410, by = .015),
  northing = seq(.150, 3.645, by = .015))
soi.grid <-  expand.grid(soi.marginal)
soi.m <- loess(resistivity~easting*northing, span = 0.25, degree = 2)
soi.fit <- predict(soi.m, soi.grid)
given <- seq(min(soi.fit), max(soi.fit), length = 29)
Resistivity <- shingle(as.vector(soi.fit), cbind(given[-29],given[-1]))
ans <- xyplot(soi.grid$n ~ soi.grid$e | Resistivity,
  layout = c(7, 4),
  panel = function(x, y) {
    panel.grid(v=2)
    points(x, y, pch = ".")
  },
  aspect=diff(range(soi.grid$n))/diff(range(soi.grid$e)),
  xlab = "Easting (km)",
  ylab = "Northing (km)")
detach()
print(ans, position=c(0,.25,1,1), more=T)


attach(soil)
soi.marginal <- list(easting = seq(.15, 1.410, by = .015),
  northing = seq(.150, 3.645, by = .015))
soi.grid <-  expand.grid(soi.marginal)
soi.m <- loess(resistivity~easting*northing, span = 0.25, degree = 2)
soi.fit <- predict(soi.m, soi.grid)
given <- seq(min(soi.fit), max(soi.fit), length = 29)
ans <- plot(shingle(as.vector(soi.fit), 
  intervals = cbind(given[-29],given[-1])), 
  cex=.7,
  sub = list("Figure 4.67", cex = 0.8),
  scales=list(x=list(cex=.7),
    y=list(cex=.7,at=seq(5,25,5))),
  xlab = "Resistivity (ohm-cm)")
detach()
print(ans, position=c(0,.05,1,.3))

invisible()

Figure 4.68

attach(soil)
soi.grid <- expand.grid(easting = seq(.15, 1.410, by = .015),
  northing = seq(.150, 3.645, by = .015))
soi.fit <- predict(loess(resistivity~easting*northing, span = 0.25, 
  degree = 2), soi.grid)
ans <- levelplot(soi.fit ~ soi.grid$easting * soi.grid$northing,
  cuts = 9,
  colorkey = list(labels = list(at = c(20, 40, 60, 80)), 
    width = 5),
  aspect=diff(range(soi.grid$n))/diff(range(soi.grid$e)),
  sub = list("Figure 4.68",cex=.8),
  xlab = "Easting (km)",
  ylab = "Northing (km)")
detach()
ans

How to cite this page

Report an error on this page

UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.

Splus Textbook Examples Visualizing Data: Chapter 4, Trivariate Data

Splus Textbook Examples
Visualizing Data: Chapter 4, Trivariate Data