book1.1 <-
function()
{
trellis.device(postscript, horizontal=T, file = "book1.1.ps", paper="legal")
print(
dotplot(variety ~ yield | site, 
data = barley,
groups = year,
layout = c(1,6),
aspect = .5,
main = list("Figure 1.1",cex=1.25),
scales=list(y=list(cex=.6)),
xlab = "Barley Yield (bushels/acre)",
key = list(points = Rows(trellis.par.get("superpose.symbol"), 1:2),
text = list(levels(barley$year)), columns = 2))
)
graphics.off()
}

book1.2 <-
function()
{
trellis.device(postscript, file = "book1.2.ps", color = "TRUE", paper="legal")
print(
	histogram(~ height | voice.part,
		  data = singer, 
		  nint = 17, 
		  layout = c(8, 1), 
		  main = list("Figure 1.2",cex=1),
		  xlab = "Height (inches)")
)
graphics.off()
}

book1.3 <-
function()
{
trellis.device(postscript, file = "book1.3.ps", color = "TRUE", paper="legal")
print(
	xyplot(babinet ~ concentration, 
	       data = polarization, 
	       aspect = 1,
		sub = list("Figure 1.3",cex=.8),
	       xlab = "Concentration (micrograms per cubic meter)",
	       ylab = "Babinet Point (degrees)")
)
graphics.off()    
}

book1.4 <-
function()
{
trellis.device(postscript, file = "book1.4.ps", color = "TRUE", paper="legal")
print(
	xyplot(NOx ~ E, 
	       data = ethanol, 
	       aspect = 1, 
		sub = list("Figure 1.4",cex=.8),
	       xlab= "E", 
	       ylab= "Oxides of Nitrogen\n(micrograms/J)")
)
graphics.off()
}

book1.5 <-
function()
{
trellis.device(postscript, file = "book1.5.ps", color = "TRUE", paper="legal")
print(
	xyplot(NOx ~ C, 
	       data = ethanol, 
	       aspect = 1,
		sub = list("Figure 1.5",cex=.8),
	       xlab= "C", 
	       ylab= "Oxides of Nitrogen\n(micrograms/J)")
)
graphics.off()
}

book1.6 <-
function()
{
trellis.device(postscript, file = "book1.6.ps", color = "TRUE", paper="legal")
print(
	dotplot( variety ~ yield | year * site, 
		data = barley, 
		aspect = .4, 
		sub = list("Figure 1.6",cex=.8),
		xlab = "Barley Yield (bushels/acre)")
)
graphics.off()
}

book2.1 <-
function()
{
trellis.device(postscript, file = "book2.1.ps", color = TRUE, paper = "legal")
graph <- qqmath(~ height | voice.part,
distribution=qunif,
data=singer,
panel = function(x, ...)
{
panel.grid()		
panel.qqmath(x,...)
},
layout=c(8,1), 
aspect=1,
sub = list("Figure 2.1",cex=1),
xlab = "f-value",
ylab="Height (inches)"
)
print(graph)
graphics.off()
}

book2.2 <- 
function(){
trellis.device(postscript, file = "book2.2.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ singer$height[singer$voice.part=="Tenor 1"],
		distribution = qunif,
                type="b",
		aspect = 1, 
		sub = list("Figure 2.2",cex=1),
		xlab = "f-value",
		ylab = "Tenor 1 Height (inches)")
)
graphics.off()
}

book2.3 <-
function()
{
voice.part <- ordered(singer$voice.part, 
	c("Soprano 1", "Soprano 2", "Alto 1", "Alto 2",
	"Tenor 1", "Tenor 2", "Bass 1", "Bass 2"))
trellis.device(postscript, file = "book2.3.ps", color = TRUE, paper = "legal")
print(
	qq(voice.part ~ height,
		subset=voice.part=="Tenor 1" | voice.part=="Bass 2",
		aspect=1,
		data = singer,
		sub = list("Figure 2.3",cex=1),
		xlab = "Tenor 1 Height (inches)",
		ylab = "Base 2 Height (inches)")
)
graphics.off()
}

book2.4 <-
function()
{
voice.part <- ordered(singer$voice.part,
	c("Soprano 1", "Soprano 2", "Alto 1", "Alto 2",
	"Tenor 1", "Tenor 2", "Bass 1", "Bass 2"))
bass.tenor.qq <- qq(voice.part ~ singer$height,
		subset=voice.part=="Bass 2" | voice.part=="Tenor 1")
trellis.device(postscript, file = "book2.4.ps", color = TRUE, paper = "legal")
print(
	tmd(bass.tenor.qq,
		aspect=1,
		ylab = "Difference (inches)",
		sub = list("Figure 2.4",cex=1),
		xlab = "Mean (inches)")
)
graphics.off()
}

book2.6 <-
function()
{
trellis.device(postscript, file = "book2.6.ps", color = TRUE, paper = "legal")
	oldpty <- par("pty")
	par(pty = "s")
	data <-
	c(0.9, 1.6, 2.26305, 2.55052, 2.61059, 2.69284, 2.78511, 2.80955, 
		2.94647, 2.96043, 3.05728, 3.15748, 3.18033, 3.20021, 
		3.20156, 3.24435, 3.33231, 3.34176, 3.3762, 3.39578, 3.4925,
		3.55195, 3.56207, 3.65149, 3.72746, 3.73338, 3.73869, 
		3.80469, 3.85224, 3.91386, 3.93034, 4.02351, 4.03947, 
	        4.05481, 4.10111, 4.26249, 4.28782, 4.37586, 4.48811, 
		4.6001, 4.65677, 4.66167, 4.73211, 4.80803, 4.9812, 5.17246,
		5.3156, 5.35086, 5.36848, 5.48167, 5.68, 5.98848, 6.2, 7.1, 
		7.4)
	boxplot(data, rep(NA, length(data)), ylab = "Data")
	usr <- par("usr")
	x <- usr[1] + (usr[2] - usr[1]) * 0.5
	at <- c(0.9, 1.6, 3.2, 3.8, 4.65, 6.2, 7.2)
	arrows(rep(x * 1.15, 7), at, rep(x, 7), at)
	mtext("Figure 2.6",1,1,cex=.8)
	text(rep(x * 1.2, 7), at, adj = 0,
		labels = c("outside value", "lower adjacent value", 
			"lower quartile", "median", "upper quartile", 
			"upper adjacent value", "outside values"))	
	par(pty = oldpty)
	invisible()
graphics.off()
}

	
book2.7 <-
function()
{
	data <- round(c(0.9, 1.6, 2.263047,
		2.550518, 2.610592, 2.69284, 2.785113, 
		2.809547, 2.946467, 2.96044, 3.057283, 
		3.15748, 3.180327, 3.200206, 
		3.20156, 3.244347, 3.332312, 
		3.341763, 3.3762, 3.395778, 3.492497, 
		3.551945, 3.562066, 3.65149, 
		3.7274632, 3.73338, 3.738686, 3.80469, 
		3.85224, 3.91386, 3.93034, 
		4.02351, 4.039466, 4.05481, 4.101108, 4.262486, 
		4.28782, 4.375864, 4.48811, 4.6001, 
		4.656775, 4.661673, 4.73211, 
		4.80803, 4.9812, 5.172464, 
		5.3156, 5.35086, 5.36848, 
		5.48167, 5.68, 5.98848, 6.2, 
		7.1, 7.4),5)
	uq <- quantile(data,.75)
	lq <- quantile(data,.25)
	r <- 1.5*(uq-lq)
	h <- c(lq-r,1.6,lq,uq,6.2,uq+r)
	writing <- c("lower quartile - 1.5 r",
		"lower adjacent value",
		"lower quartile",
		"upper quartile",
		 "upper adjacent value",
		 "upper quartile + 1.5 r")
trellis.device(postscript, file = "book2.7.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ data,
		distribution = qunif,
		panel = substitute(function(x, ...) {
			reference.line <- trellis.par.get("reference.line")
			panel.abline(h = h, lwd = reference.line$lwd, lty =
reference.line$lty, col = reference.line$col)
			panel.qqmath(x, ..., col = 0, pch = 16)
			panel.qqmath(x, ...)
			panel.text(rep(0,3), h[4:6], writing[4:6], adj=0)
			panel.text(rep(1,3), h[1:3], writing[1:3], adj=1)
		}),
		aspect = 1, 
		sub = list("Figure 2.7",cex=.8),
		xlab = "f-value", 
		ylab = "Data")
)
graphics.off()
}

book2.8 <-
function(){
trellis.device(postscript, file = "book2.8.ps", color = TRUE, paper = "legal")
print(
	bwplot(voice.part ~ height,
		data=singer,
		aspect=1,
		sub = list("Figure 2.8",cex=1),
		xlab="Height (inches)")
)
graphics.off()
}

book2.9 <-
function()
{
data <- sort(singer$height[singer$voice.part=="Alto 1"])
trellis.device(postscript, file = "book2.9.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ data, 
		distribution = qunif,type="b",
		panel = function(x, ...) {
			panel.grid()
			panel.qqmath(x, ..., )
		},
		aspect = 1, 
		ylim = range(data, qnorm(ppoints(data), mean(data),
sqrt(var(data)))),
		sub = list("Figure 2.9",cex=1),
		xlab = "f-value",
		ylab = "Alto 1 Height (inches)")
)
graphics.off()
}

book2.10 <-
function()
{
data <- sort(singer$height[singer$voice.part=="Alto 1"])
x <- ppoints(data)
y <- qnorm(x, mean(data), sqrt(var(data)))
trellis.device(postscript, file = "book2.10.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ x, 
		panel = function(x, y){
			panel.grid()
			panel.xyplot(x, y, type = "l")
		},
		ylim = range(data, y),
		aspect = 1, 
		sub = list("Figure 2.10",cex=1),
		xlab = "f-value", 
		ylab = "Normal Quantile Function")
)
graphics.off()
}


book2.11 <-
function(){
trellis.device(postscript, file = "book2.11.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ height | voice.part,
		data=singer,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		layout=c(8,1),
		aspect=1, 
		sub = list("Figure 2.11",cex=1),
		xlab = "Unit Normal Quantile",
		ylab="Height (inches)")
)
graphics.off()
}

book2.12 <-
function(){
trellis.device(postscript, file = "book2.12.ps", color = TRUE, paper = "legal")
print(
	dotplot(tapply(singer$height,singer$voice.part,mean), 
		aspect=1,
		sub = list("Figure 2.12",cex=1),
		xlab="Mean Height (inches)")
)
graphics.off()
}

book2.13 <-
function(){
trellis.device(postscript, file = "book2.13.ps", color = TRUE, paper = "legal")
print(
	bwplot(voice.part ~ oneway(height~voice.part, spread = 1)$residuals,
		data = singer,
		aspect=0.75,
		panel = function(x,y){
			panel.abline(v=0)
			panel.bwplot(x,y)
		},
		sub = list("Figure 2.13",cex=1),
		xlab = "Residual Height (inches)")
)
graphics.off()
}

book2.14 <-
function()
{
res.height <- oneway(height ~ voice.part, data = singer, spread = 1)$residuals
trellis.device(postscript, file = "book2.14.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ res.height | singer$voice.part, 
		distribution = substitute(function(p) quantile(res.height, p)),
		panel=function(x,...){
			panel.grid()
			panel.abline(0, 1)
			panel.qqmath(x, ...)
		},
		aspect=1,
		layout=c(8,1),
		sub = list("Figure 2.14",cex=1),
		xlab = "Pooled Residual Height (inches)",
		ylab = "Residual Height (inches)")
)
graphics.off()
}

book2.15 <-
function(){
trellis.device(postscript, file = "book2.15.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ oneway(height ~ voice.part, spread = 1)$residuals, 
		data = singer,
		distribution = qunif,
		aspect = 1,
		sub = list("Figure 2.15",cex=1),
		xlab = "f-value",
		ylab = "Residual Height (inches)")
)
graphics.off()
}

book2.16 <-
function(){
trellis.device(postscript, file = "book2.16.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ oneway(height~voice.part, spread = 1)$residuals,
		data = singer,
		prepanel = prepanel.qqmathline, 
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		sub = list("Figure 2.16",cex=1),
		xlab = "Unit Normal Quantile",
		ylab="Residual Height (inches)")
)
graphics.off()
}

book2.17 <-
function(){
trellis.device(postscript, file = "book2.17.ps", color = TRUE, paper = "legal")
print(
	rfs(oneway(height~voice.part, data = singer, spread = 1), 
		aspect=1, 
		sub = list("Figure 2.17",cex=1),
		ylab = "Height (inches)")
)
graphics.off()
}

book2.19 <-
function()
{
trellis.device(postscript, file = "book2.19.ps", color = TRUE, paper = "legal")
plot<-	qqmath(~ time | nv.vv,
		data=fusion.time,
		distribution = qunif,
	        panel = function(x, y, ...) {
			panel.grid()
	                panel.qqmath(x, y, ...)
	        },
		aspect=1,
		layout=c(2,1),
		sub = list("Figure 2.19",cex=.8),
		xlab = "f-value",
		ylab="Time (seconds)")
print(plot)
graphics.off()
}


book2.20 <-
function()
{
data <- 5+qnorm(ppoints(25))
trellis.device(postscript, file = "book2.20.ps", color = TRUE, paper = "legal")
print(
	qqmath(~data,
		distribution = qunif,type="b",
		panel = function(x, ...) {
			reference.line <- trellis.par.get("reference.line")
			m <- median(x)
			panel.segments(c(.1, .9), m, c(.1, .9), quantile(x, c(.1,.9)), lwd = reference.line$lwd, lty =reference.line$lty, col = reference.line$col)
			panel.qqmath(x, col = 0, pch = 16)
			panel.qqmath(x,... )
			panel.abline(h = m)
			panel.text(.05, 4.25, "d(0.1)", srt = 90, adj = 0)
			panel.text(.85, 5.25, "d(0.9)", srt = 90, adj = 0)
		},
		aspect = 1, 
		sub = list("Figure 2.20",cex=.8),
		xlab = "f-value",
	        ylab = "Data")
)
graphics.off()
}


book2.21 <-
function()
{
data <- 2 ^ (5 + qnorm(ppoints(25)))
trellis.device(postscript, file = "book2.21.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ data,
		distribution = qunif, type="b",
		panel = function(x, ...) {
			reference.line <- trellis.par.get("reference.line")
			m <- median(x)
			panel.segments(c(.1, .9), m, c(.1, .9), quantile(x, c(.1, .9)),
				lwd = reference.line$lwd, lty = reference.line$lty, col = reference.line$col)
			panel.qqmath(x, ..., col = 0, pch = 16)
			panel.qqmath(x, ...)
			panel.abline(h = m)
	                panel.text(.05, 15, "d(0.1)", srt = 90, adj = 0)
	                panel.text(.85, 40, "d(0.9)", srt = 90, adj = 0)
		},
		aspect = 1, 
		sub = list("Figure 2.21",cex=.8),
		xlab = "f-value",
	        ylab = "Data")
)
graphics.off()
}

book2.22 <-
function()
{
trellis.device(postscript, file = "book2.22.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ time | nv.vv,
		data=fusion.time,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		layout=c(2,1),
		sub = list("Figure 2.22",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab="Time (seconds)")
)
graphics.off()
}

book2.23 <-
function()
{
trellis.device(postscript, file = "book2.23.ps", color = TRUE, paper = "legal")
print(
	bwplot(nv.vv ~ time,
		data=fusion.time,
		aspect = .5,
		sub = list("Figure 2.23",cex=.8),
		xlab="Time (seconds)")
)
graphics.off()
}

book2.24 <-
function(){
trellis.device(postscript, file = "book2.24.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ log(time,2) | nv.vv,
		data=fusion.time,
		prepanel=prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		layout=c(2,1),
		sub = list("Figure 2.24",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab="Log Time (log 2 seconds)")
)
graphics.off()
}

book2.25 <-
function()
{
fusion.time.m <- oneway(time ~ nv.vv, data=fusion.time,location=median, spread=1)
trellis.device(postscript, file = "book2.25.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(fusion.time.m)))~jitter(fitted.values(fusion.time.m),amount=.2,factor=3),
		aspect=1,
		panel=substitute(function(x,y){
			panel.xyplot(x,y)
			srmads <- sqrt(tapply(abs(residuals(fusion.time.m)),
				fusion.time$nv.vv, median))
			panel.lines(fusion.time.m$location,srmads)
		}),
		sub = list("Figure 2.25",cex=.8),
		xlab="Jittered Median Time (sec)",
		ylab="Square Root Absolute Residual Time (square root sec)")
)
graphics.off()
}

book2.26 <-
function()
{
fusion.time.m <- oneway(log(time,2) ~ nv.vv,data=fusion.time, location = median, spread=1)
trellis.device(postscript, file = "book2.26.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(fusion.time.m))) ~ jitter(fitted.values(fusion.time.m),factor=3,amount=.2),
		aspect=1,
		panel=substitute(function(x,y){
			panel.xyplot(x,y)
			srmads <- tapply(abs(residuals(fusion.time.m)),
				fusion.time$nv.vv,median)
			panel.lines(fusion.time.m$location,srmads)
		}),
		sub = list("Figure 2.26",cex=.8),
		xlab="Jittered Median Log Time (log 2 sec)",
		ylab="Square Root Absolute Residual Log Time (square root absolute log 2 sec)")
)
graphics.off()
}

book2.27 <-
function()
{
trellis.device(postscript, file = "book2.27.ps", color = TRUE, paper = "legal")
print(
	qq(nv.vv ~ time,
		data = fusion.time,
		aspect = 1,
		sub = list("Figure 2.27",cex=.8),
		xlab="NV Time (seconds)",
		ylab="VV Time (seconds)")
)
graphics.off()
}

book2.28 <-
function()
{
trellis.device(postscript, file = "book2.28.ps", color = TRUE, paper = "legal")
print(
	qq(nv.vv ~ log(time, 2),
		data = fusion.time,
		aspect = 1,
		sub = list("Figure 2.28",cex=.8),
		xlab = "Log NV Time (log 2 seconds)",
		ylab = "Log VV Time (log 2 seconds)")
)
graphics.off()
}

book2.29 <-
function()
{
trellis.device(postscript, file = "book2.29.ps", color = TRUE, paper = "legal")
print(
	tmd(qq(nv.vv ~ log(time, 2), data = fusion.time),
		aspect = 1,
		sub = list("Figure 2.29",cex=.8),
		xlab = "Mean (log 2 seconds)",
		ylab = "Difference (log 2 seconds)")
)
graphics.off()
}

book2.30 <-
function()
{
trellis.device(postscript, file = "book2.30.ps", color = TRUE, paper = "legal")
res <- oneway(log(time,2)~nv.vv, data = fusion.time, spread = 1)$residuals
print(
	qqmath(~ res | fusion.time$nv.vv,
		distribution = substitute(function(p){quantile(res, p)}),
		panel=function(x,...){
			panel.grid()
			panel.abline(0, 1)
			panel.qqmath(x, ...)
		},
		aspect=1,
		layout=c(2,1),
		sub = list("Figure 2.30",cex=.8),
		xlab = "Pooled Residual Log Time (log 2 seconds)",
		ylab = "Residual Log Time (log 2 seconds)")
)
graphics.off()
}

book2.31 <-
function()
{
trellis.device(postscript, file = "book2.31.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ oneway(log(time,2)~nv.vv, data = fusion.time, spread = 1)$residuals,
		prepanel = prepanel.qqmathline, 
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect = 1,
		sub = list("Figure 2.31",cex=.8),
		xlab = "Unit Normal Quantile", 
		ylab = "Residual Log Time (log 2 seconds)")
)
graphics.off()
}

book2.32 <-
function()
{
trellis.device(postscript, file = "book2.32.ps", color = TRUE, paper = "legal")
print(
	rfs(oneway(log(time, 2)~nv.vv, data = fusion.time, spread = 1),
		aspect=1,
		sub = list("Figure 2.32",cex=.8),
	        ylab = "Log Time (log 2 seconds)")
)
graphics.off()
}

book2.33 <-
function()
{
vvtime <- fusion.time$time[fusion.time$nv.vv=="VV"]
lambda.factor <- factor(rep(c("-1","-1/2", "-1/4","0","1/4", "1/2", "1"), rep(length(vvtime),7)))
lambda.factor <- ordered(lambda.factor, c("-1","-1/2", "-1/4","0","1/4", "1/2", "1"))

#method 1
powfun <- function(x, lambda) x^lambda/lambda
transformed <- c(cbind(outer(vvtime,c(-1,-1/2,-1/4), powfun),log(vvtime),
	(outer(vvtime,c(1/4,1/2,1),powfun))))

#method 2
powfun <- function(x, lambda) if(lambda==0) log(x) else x^lambda/lambda
lambda <- c(-1,-1/2,-1/4,0,1/4,1/2,1)
tmp <- expand.grid(vvtime, lambda)
names(tmp) <- c("x", "lambda")
transformed <- apply(tmp, 1, function(r) do.call(powfun, as.list(r)))

trellis.device(postscript, file = "book2.33.ps", color = TRUE, paper = "legal")
print(
	ans <- qqmath(~ transformed | lambda.factor,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid(h = 0)
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		scale = list(y = "free"),
		layout=c(4,2), 
		sub = list("Figure 2.33 Revised",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "VV Time")
)
graphics.off()
}

book2.34 <-
function()
{
trellis.device(postscript, file = "book2.34.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ mean.length | dimension,
		distribution = qunif,
		data=food.web,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmath(x, ...)
	        },
		layout=c(3,1),
		aspect=1,
		sub = list("Figure 2.34",cex=.8),
		xlab = "f-value",
		ylab="Chain Length")
)
graphics.off()
}

book2.35 <-
function()
{
foo.m <- oneway(mean.length~dimension, data = food.web, location = median, spread=1)
set.seed(19)
trellis.device(postscript, file = "book2.35.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(foo.m))) ~ jitter(fitted.values(foo.m),factor=3,amount=.03),
		aspect=1,
	        panel = substitute(function(x,...){
			panel.xyplot(x,...)
			srmads <- tapply(abs(residuals(foo.m)),
				food.web$dimension, median)
			panel.lines(foo.m$location,srmads)
		}),
		sub = list("Figure 2.35",cex=.8),
		xlab="Jittered Median Chain Length",
		ylab="Square Root Absolute Residual Chain Length")
)
graphics.off()
}


book2.36 <-
function()
{
trellis.device(postscript, file = "book2.36.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ mean.length | dimension,
		data=food.web,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		layout=c(3,1),
		aspect=1, 
		sub = list("Figure 2.36",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab="Chain Length")
)
graphics.off()
}

book2.37 <-
function()
{
	foo.m <- oneway(log(mean.length, 2) ~ dimension, data = food.web, location = median, spread = 1)
	set.seed(19)
trellis.device(postscript, file = "book2.37.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(foo.m))) ~ jitter(fitted.values(foo.m),amount=.03, factor = 3),
	        panel = substitute(function(x, ...) {
			panel.xyplot(x, ...)
			add.line <- trellis.par.get("add.line")
			panel.lines(foo.m$location, tapply(y, food.web$dimension, median), 
				lty = add.line$lty, lwd = add.line$lwd, col = add.line$col)
		}),
		aspect = 1,
		sub = list("Figure 2.37",cex=.8),
		xlab = "Jittered Median Log 2 Chain Length",
		ylab = "Square Root Absolute Residual Log 2 Chain Length")
)
graphics.off()
}

book2.38 <-
function()
{
trellis.device(postscript, file = "book2.38.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ log(mean.length,2) | dimension,
		data=food.web,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		layout=c(3,1),
		aspect=1, 
		sub = list("Figure 2.38",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab="Log 2 Chain Length")
)
graphics.off()
}

book2.39 <-
function()
{
foo.m <- oneway(1/mean.length ~ dimension, data = food.web, location = median, spread = 1)
set.seed(19)
trellis.device(postscript, file = "book2.39.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(foo.m))) ~ jitter(fitted.values(foo.m), amount=.004, factor = 3),
	        panel = substitute(function(x,...) {
			panel.xyplot(x,...)
	                add.line <- trellis.par.get("add.line")
			panel.lines(foo.m$location, tapply(y, food.web$dimension, median),
			lty = add.line$lty, lwd = add.line$lwd, col = add.line$col)
		}),
		aspect = 1,
		sub = list("Figure 2.39",cex=.8),
		xlab = "Jittered Median Link Fraction",
		ylab = "Square Root Absolute Residual Link Fraction")
)
graphics.off()
}

book2.40 <-
function()
{
trellis.device(postscript, file = "book2.40.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ (1/mean.length) | dimension,
	       data = food.web,
	       panel = function(x, ...){
			panel.grid()
			panel.qqmath(x, ...)
			panel.qqmathline(x, distribution = qnorm)
	       },
	       layout = c(3, 1),
	       aspect = 1,
		sub = list("Figure 2.40",cex=.8),
	       xlab = "Unit Normal Quantile", 
	       ylab = "Link Fraction")
)
graphics.off()
}

book2.41 <-
function()
{
res <- oneway((1/mean.length)~dimension, data = food.web, spread = 1)$residuals
trellis.device(postscript, file = "book2.41.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ res | food.web$dimension,
		distribution = substitute(function(p) quantile(res, p)),
		panel=function(x,...){
			panel.grid()
			panel.abline(0, 1)
			panel.qqmath(x, ...)
		},
		layout=c(3,1),
		aspect=1,
		sub = list("Figure 2.41",cex=.8),
		xlab = "Pooled Residual Link Fraction",
		ylab = "Residual Link Fraction")
)
graphics.off()
}

book2.42 <-
function()
{
trellis.device(postscript, file = "book2.42.ps", color = TRUE, paper = "legal")
print(
	rfs(oneway((1/mean.length)~dimension, data = food.web, spread = 1),
		sub = list("Figure 2.42",cex=.8),
		aspect=1,
		ylab = "Link Fraction")
)
graphics.off()
}

book2.43 <-
function()
{
trellis.device(postscript, file = "book2.43.ps", color = TRUE, paper = "legal")
print(
	bwplot(factor(number.runs) ~ log(empty.space,2),
		data=bin.packing,
		aspect=1,
		sub = list("Figure 2.43",cex=.8),
		xlab="Log 2 Empty Space")
)
graphics.off()
}

book2.44 <-
function()
{
trellis.device(postscript, file = "book2.44.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ log(empty.space,2) | factor(number.runs),
		data = bin.packing,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		layout = c(4, 3),
		aspect = 1,
		sub = list("Figure 2.44",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Log 2 Empty Space")
)
graphics.off()
}

book2.45 <-
function()
{
res <- oneway(log(empty.space,2) ~ number.runs, 
		data = bin.packing, 
		location = median,
# This next line is the way it should be.
#		spread = function(x) median(abs(x-median(x)))
# This next line is the way it is now.
		spread = function(x)(quantile(x,.75)-quantile(x,.25))/1.33
	      )$scaled.residuals
trellis.device(postscript, file = "book2.45.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ res | factor(bin.packing$number.runs),
		prepanel = prepanel.qqmathline, 
	        panel = function(x, ...) {
			panel.grid()
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
        	},
		aspect=1,
		layout = c(4,3),
		sub = list("Figure 2.45",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Spread-Standardized Residual Log 2 Empty Space")
)
graphics.off()
}

book2.46 <-
function()
{

data <- log(empty.space,2)
bin.m <- oneway(data~number.runs,
		location = median,
		spread=function(x){diff(quantile(x,c(.25,.75)))/1.35}
)
res <- bin.m$scaled.res[number.runs>1000]
gr <- factor(number.runs[number.runs>1000])
trellis.device(postscript, file = "book2.46.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ res | gr,
		distribution = substitute(function(p) quantile(res, p)),
		panel = function(x, ...){
			panel.grid()
			panel.abline(0, 1)
			panel.qqmath(x, ...)
		},
		aspect = 1,
		layout = c(4, 2),
		sub = list("Figure 2.46",cex=.8),
		xlab = "Pooled Spread-Standardized Residual Log 2 Empty Space",
		ylab = "Spread-Standardized Residual Log 2 Empty Space")
)
graphics.off()

}

book2.47 <-
function()
{
bin.m <- oneway(log(empty.space,2) ~ number.runs, 
		data = bin.packing, 
		location = median, 
		spread = function(x){diff(quantile(x,c(.25,.75)))/1.35}
)
trellis.device(postscript, file = "book2.47.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ bin.m$scaled.res[bin.packing$number.runs > 1000],
		prepanel = prepanel.qqmathline, 
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect = 1, 
		sub = list("Figure 2.47",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Spread-Standardized Residual Log 2 Empty Space")
)
graphics.off()
}

book2.48 <-
function()
{
	
	data <- log(empty.space,2)
	mq <- tapply(data, number.runs, median)
	nw <- log(sort(unique(number.runs)),2)
trellis.device(postscript, file = "book2.48.ps", color = TRUE, paper = "legal")
print(
	xyplot(mq ~ nw,
	        panel = function(x, y){
			panel.xyplot(x, y)
			panel.abline(y[11] - x[11]/3, 1/3)
		},
		aspect = 1,
		sub = list("Figure 2.48",cex=.8),
		xlab = "Log 2 Number of Weights",
		ylab = "Median Log 2 Empty Space")
)
graphics.off()

}

book2.49 <-
function()
{

bin.packing.m <- oneway(log(empty.space, 2) ~ number.runs, 
			location = median, 
			spread = 1
)
srmads <- tapply(abs(residuals(bin.packing.m)), number.runs, median)
trellis.device(postscript, file = "book2.49.ps", color = TRUE, paper = "legal")
print(
	xyplot(log(srmads, 2) ~ log(sort(unique(number.runs)), 2),
		aspect = 1,
		sub = list("Figure 2.49",cex=.8),
		xlab = "Log 2 Number of Weights",
		ylab = "Log 2 Mad of Log 2 Empty Space")
)
graphics.off()

}

book2.50 <-
function()
{

bin.packing.m <- oneway(log(empty.space,2) ~ number.runs, 
			location = median, 
			spread = 1
)
srmads <- tapply(abs(residuals(bin.packing.m)), number.runs, median)
trellis.device(postscript, file = "book2.50.ps", color = TRUE, paper = "legal")
print(
	xyplot(log(srmads/min(srmads), 2) ~ bin.packing.m$location,
		aspect = 1,
		sub = list("Figure 2.50",cex=.8),
		xlab = "Median Log 2 Empty Space",
		ylab = "Log 2 Relative Spread")
)
graphics.off()

}

book2.51 <-
function()
{

bin.packing.m <- oneway(empty.space~number.runs,
			location = median,
			spread=1
)
srmads <- tapply(abs(residuals(bin.packing.m)),number.runs,median)
log.srmads <- log(srmads/min(srmads),2)
trellis.device(postscript, file = "book2.51.ps", color = TRUE, paper = "legal")
print(
xyplot(log.srmads ~ bin.packing.m$location,
		aspect=1,
		sub = list("Figure 2.51",cex=.8),
		xlab="Median Empty Space",
		ylab="Log 2 Relative Spread")
)
graphics.off()

}

book2.52 <-
function()
{
res <- oneway(log(empty.space,2) ~ number.runs, 
		data = bin.packing,
		location = median,
		spread = function(x){median(abs(x-median(x)))}
	)$scaled.residuals/1.68
trellis.device(postscript, file = "book2.52.ps", color = TRUE, paper = "legal")
print(
	qqmath(~ res | factor(bin.packing$number.runs), 
		prepanel = prepanel.qqmathline, 
		panel = function(x, ...) {
			panel.grid()
			panel.qqmathline(x, distribution = qnorm)
			panel.qqmath(x, ..., err=-1)   # no warnings for out of bounds
		},
		ylim = c(-2.75, 2.75),
		aspect= 1,
		layout = c(4, 3), 
		sub = list("Figure 2.52",cex=.8),
		xlab = "Unit Normal Quantile", 
		ylab = "Spread-Standardized Log 2 Empty Space")
)
graphics.off()
}



book3.1 <- 
function()
{
trellis.device(postscript, file = "book3.1.ps", color = TRUE, paper = "legal")
print(
	xyplot(cp.ratio ~ area,
		data = ganglion,
		aspect=1,
		sub = list("Figure 3.1",cex=.8),
		xlab = "Area (square mm)",
		ylab = "CP Ratio")
)
graphics.off()
}

book3.10 <- 
function()
{
X <- seq(5, 10, length=150)
z <- seq(0, 1, length=50)
Y <- 5 + c(-(1-z)^.5, sin(z*2*pi), .5*z^2)
base.y <- 5 + c(-(1-z)^.5, sin(z*2*pi), .5*z^2)
true.y <- loess.smooth(X, base.y, degree = 2, family = "g", span = .5,
evaluation = 150)$y
set.seed(20)
Y <- true.y+rnorm(150, 0, .2)
fit1 <- loess.smooth(X, Y, degree = 1, family = "g", span = .1, evaluation = 150)
fit2 <- loess.smooth(X, Y, degree = 1, family = "g", span = .3, evaluation = 150)
fit3 <- loess.smooth(X, Y, degree = 2, family = "g", span = .3, evaluation = 150)
fits <- c(fit1$y, fit2$y, fit3$y)
which <- factor(rep(c("alpha .1, lambda 1", "alpha .3, lambda 1", "alpha .3, lambda 2"), rep(length(X),3)))
X <- rep(X,3)
Y <- rep(Y,3)
trellis.device(postscript, file = "book3.10.ps", color = TRUE, paper = "legal")
print(
	xyplot(Y ~ X | which, 
		prepanel = substitute(function(x, y, subscripts) {
			list(dx = diff(x), dy=diff(fits[subscripts]))}),
		panel = substitute(function(x, y, subscripts){
			add.line <- trellis.par.get("add.line")
			panel.xyplot(x, y)
			panel.lines(x, fits[subscripts], lwd = add.line$lwd,
				lty = add.line$lty, col = add.line$col)
		}),
		layout = c(1, 3),
		sub = list("Figure 3.10",cex=.8),
		xlab = "x", 	
		ylab = "y")
)
graphics.off()
}


book3.11 <- 
function()
{
x <- seq(5, 10, length=150)
z <- seq(0, 1, length=50)
base.y <- 5 + c(-(1-z)^.5, sin(z*2*pi), .5*z^2)
true.y <- loess.smooth(x,base.y,degree = 2, family = "g", span = .5, evaluation = 150)$y
set.seed(20)
y <- true.y+rnorm(150, 0, .2)
fit <- loess.smooth(x, y, degree = 2, family = "g", span = .5, evaluation = 150)
trellis.device(postscript, file = "book3.11.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ fit$x, 
		prepanel = function(x, y) 
			prepanel.loess(x, y, degree = 2, family = "g", span = .3, evaluation = 300), 
		panel = function(x, y){
	        	panel.xyplot(x, y)
			panel.loess(x, y, degree = 2, family = "g", span = .3, evaluation = 300)
		        panel.abline(v = c(5.3, 5.7))
		},
		sub = list("Figure 3.11",cex=.8),
		xlab = "x", 	
		ylab = "y")
)
graphics.off()
}

book3.12 <- 
function()
{
trellis.device(postscript, file = "book3.12.ps", color = TRUE, paper = "legal")
print(
	xyplot(lm(cp.ratio~area)$residuals ~ area,
		data = ganglion,
		prepanel = function(x, y)
			prepanel.loess(x, y, span=1, family = "g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x, y, span=1, family = "g")
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.12",cex=.8),
		xlab = "Area (square mm)", 
		ylab="Residual CP Ratio")
)
graphics.off()
}

book3.13 <- 
function()
{
trellis.device(postscript, file = "book3.13.ps", color = TRUE, paper = "legal")
print(
	xyplot(lm(cp.ratio~area+I(area^2))$residuals ~ area,
		data = ganglion,
		prepanel = function(x, y)
			prepanel.loess(x, y, span=1, family = "g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x, y, span=1, family = "g")
			panel.abline(h=0)
		},	
		aspect=1,
		sub = list("Figure 3.13",cex=.8),
		xlab = "Area (square mm)",
		ylab="Residual CP Ratio")
)
graphics.off()
}

book3.14 <- 
function()
{
gan.lm <- lm(cp.ratio~area+I(area^2), data = ganglion)
trellis.device(postscript, file = "book3.14.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(gan.lm))) ~ gan.lm$fit,
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x, y, span=2, evaluation=100, family = "g")
		},
		aspect=1,
		sub = list("Figure 3.14",cex=.8),
		xlab = "Fitted CP Ratio",
		ylab = "Square Root Absolute Residual CP Ratio")
)
graphics.off()
}

book3.15 <- 
function()
{
trellis.device(postscript, file = "book3.15.ps", color = TRUE, paper = "legal")
print(
	xyplot(log(cp.ratio,2) ~ area,
		data = ganglion,
		prepanel=function(x,y)
			prepanel.loess(x,y,family="g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x,y, family = "g")
		},
		aspect="xy",
		sub = list("Figure 3.15",cex=.8),
		xlab = "Area (square mm)",
		ylab = "Log Base 2 CP Ratio")
)
graphics.off()
}

book3.16 <- 
function()
{
trellis.device(postscript, file = "book3.16.ps", color = TRUE, paper = "legal")
print(
	xyplot(log(cp.ratio,2) ~ area, 
		data = ganglion,
		prepanel = prepanel.lmline,
		panel = substitute(function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(lm(log(cp.ratio,2)~area, data = ganglion))
		}),
		aspect=1,
		sub = list("Figure 3.16",cex=.8),
		xlab = "Area (square mm)",
		ylab = "Log Base 2 CP Ratio")
)
graphics.off()
}

book3.17 <- 
function()
{
gan.lm <- lm(log(cp.ratio,2) ~ area, data = ganglion)
trellis.device(postscript, file = "book3.17.ps", color = TRUE, paper = "legal")
print(
	xyplot(gan.lm$res ~ ganglion$area,
		prepanel = function(x, y)
			prepanel.loess(x, y,span=1, evaluation=100, family = "g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x, y, span=1, evaluation=100, family = "g")
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.17",cex=.8),
		xlab = "Area (square mm)",
		ylab="Residual Log Base 2 CP Ratio")
)
graphics.off()
}

book3.18 <- 
function()
{
gan.lm <- lm(log(cp.ratio,2) ~ area, data = ganglion)
trellis.device(postscript, file = "book3.18.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(gan.lm$res)) ~ gan.lm$fit,
		prepanel = function(x, y)
			prepanel.loess(x, y,span=2, evaluation=100, family = "g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x, y, span=2, evaluation=100, family = "g")
		},
		aspect=1,
		sub = list("Figure 3.18",cex=.8),
		xlab = "Fitted Log Base 2 CP Ratio",
		ylab = "Square Root Absolute Residual Log Base 2 CP Ratio")
)
graphics.off()
}

book3.19 <- 
function()
{
trellis.device(postscript, file = "book3.19.ps", color = TRUE, paper = "legal")
print(
	rfs(lm(log(cp.ratio,2) ~ area, data = ganglion),
		sub = list("Figure 3.19",cex=.8),
		aspect=1,
		ylab="Log Base 2 CP Ratio")
)
graphics.off()
}

book3.2 <- 
function()
{
trellis.device(postscript, file = "book3.2.ps", color = TRUE, paper = "legal")
print(
	xyplot(cp.ratio ~ area, 
		data = ganglion,
		prepanel=function(x,y)
			prepanel.loess(x,y,family="g"),
		panel = function(x,y){
			panel.xyplot(x,y)
			panel.loess(x,y, family = "g")
		},
		sub = list("Figure 3.2",cex=.8),
		xlab = "Area (square mm)",
		ylab = "CP Ratio",
		aspect="xy")
)
graphics.off()
}


book3.20 <- 
function()
{
trellis.device(postscript, file = "book3.20.ps", color = TRUE, paper = "legal")
print(
	qqmath(~lm(log(cp.ratio,2) ~ area, data = ganglion)$residuals,
		prepanel = prepanel.qqmathline,
		panel = function(x,...){
			panel.qqmathline(x, distribution = qnorm)
			panel.qqmath(x,...)
		},
		aspect = 1,
		sub = list("Figure 3.20",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual Log Base 2 CP Ratio")
)
graphics.off()
}

book3.21 <- 
function()
{
trellis.device(postscript, file = "book3.21.ps", color = TRUE, paper = "legal")
print(
	xyplot((thorium - carbon) ~ carbon,
		data = dating,
		aspect=1, 
		sub = list("Figure 3.21",cex=.8),
		xlab = "Carbon Age (kyr BP)",
		ylab = "Thorium Age - Carbon Age (kyr BP)")
)
graphics.off()
}

book3.22 <- 
function()
{
difference <- dating$thorium - dating$carbon
trellis.device(postscript, file = "book3.22.ps", color = TRUE, paper = "legal")
print(
	xyplot(difference ~ carbon,data=dating, 
		###prepanel = prepanel.lmline,
		panel = function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(lm(difference~carbon,data=dating))
		},
		aspect = 1,
		sub = list("Figure 3.22",cex=.8),
		xlab = "Carbon Age (kyr BP)",
		ylab = "Thorium Age - Carbon Age (kyr BP)")
)
graphics.off()
}

book3.23 <- 
function()
{
trellis.device(postscript, file = "book3.23.ps", color = TRUE, paper = "legal")
print(
	xyplot(lm((thorium - carbon)~carbon, data = dating)$residuals ~ dating$carbon,
		panel = function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.23",cex=.8),
		xlab = "Carbon Age (kyr BP)", 
		ylab = "Residual Age Difference (kyr BP)")
)
graphics.off()
}

book3.25 <- 
function()
{

difference <- thorium - carbon
wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
	w[U > 1.] <- 0.
	w
}
wt <- rep(1,length(difference))
for(i in 1:10){
        dat.lm <- lm(difference~carbon,weights=wt)
        wt <- wt.bisquare(dat.lm$res/median(abs(dat.lm$res)),c=6)
}
ylim <- range(fitted.values(dat.lm),difference)
trellis.device(postscript, file = "book3.25.ps", color = TRUE, paper = "legal")
print(
	xyplot(difference~carbon,
	#	ylim = ylim,
		panel = substitute(function(x, y) {
			panel.xyplot(x, y)
			panel.abline(dat.lm)
		}),
		aspect=1,
		sub = list("Figure 3.25",cex=.8),
		xlab = "Carbon Age (kyr BP)",
		ylab = "Thorium Age - Carbon Age (kyr BP)")
)
graphics.off()

}

book3.26 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

difference <- thorium - carbon
wt <- rep(1,length(difference))
for(i in 1:10){
	dat.lm <- lm(difference~carbon,weights=wt)
	wt <- wt.bisquare(dat.lm$res/median(abs(dat.lm$res)),c=6)
}
trellis.device(postscript, file = "book3.26.ps", color = TRUE, paper = "legal")
print(
	xyplot(residuals(dat.lm)~carbon,
	        panel = function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.26",cex=.8),
		xlab = "Carbon Age (kyr BP)",
		ylab = "Residual Age Difference (kyr BP)")	
)
graphics.off()

}

book3.27 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}																      

difference <- thorium - carbon
wt <- rep(1,length(difference))
for(i in 1:10){
	dat.lm <- lm(difference~carbon,weights=wt)
	wt <- wt.bisquare(dat.lm$res/median(abs(dat.lm$res)),c=6)
}
trellis.device(postscript, file = "book3.27.ps", color = TRUE, paper = "legal")
print(
	qqmath(~residuals(dat.lm),
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		sub = list("Figure 3.27",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual Age Difference (kyr BP)",
		aspect=1)
)
graphics.off()

}

book3.28 <- 
function()
{
trellis.device(postscript, file = "book3.28.ps", color = TRUE, paper = "legal")
print(
	xyplot(babinet~concentration,
		data = polarization,
		aspect=1,
		sub = list("Figure 3.28",cex=.8),
		xlab = "Concentration (micrograms per cubic meter)",
		ylab = "Babinet Point (degrees)")
)
graphics.off()
}

book3.29 <- 
function()
{
lambda <- factor(rep(c("-1","-1/2", "-1/3","0","1/3", "1/2", "1"), rep(length(concentration),7)))
lambda <- ordered(lambda, c("-1","-1/2", "-1/3","0","1/3", "1/2", "1"))
powfun <- function(x, lambda=1){
   if(lambda==0){
	log(x)
   }
   else{
        x^lambda/lambda
   }
}
transformed <- c(powfun(concentration, lambda=-1), powfun(concentration, lambda=-1/2), powfun(concentration, lambda=-1/3),
		 powfun(concentration, lambda=0), powfun(concentration, lambda=1/3), powfun(concentration, lambda=1/2),
		 powfun(concentration, lambda=1))
trellis.device(postscript, file = "book3.29.ps", color = TRUE, paper = "legal")
print(
	qqmath(~transformed|lambda,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
			panel.grid(h = 0)
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1, 
		scale = list(y = "free"),
		layout=c(7,1),
		sub = list("Figure 3.29 Revised",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Concentration")
)
graphics.off()

}


book3.3 <- 
function()
{
x <- seq(-1, 1, length = 100)
y <- c(x[1:50]^2, x[51:100])
trellis.device(postscript, file = "book3.3.ps", color = TRUE, paper = "legal")
	ans1 <- xyplot(y ~ x,
		aspect = "xy", 
		scale = list(draw = F),
		type = "l",
		xlab = "",
		ylab = "")
	print(ans1, position = c(0, 0.2, 0.75, 1), more = T)
	ans1 <- xyplot(y ~ x,
		aspect = "xy", 
		scale = list(draw = F),
		xlim=c(-1.2,1.2),
		type = "l",
		xlab = "",
		ylab = "")
	print(update(ans1, aspect = 5), position = c(0.75, 0.2, 1, 1), 
		more = T)
	ans1 <- xyplot(y ~ x,
		aspect = "xy", 
		scale = list(draw = F),
		ylim=c(-.2,1.2),
		type = "l",
		sub = list("Figure 3.3",cex=.8),
		xlab = "",
		ylab = "")
	print(update(ans1, aspect = 0.05), position = c(0, 0, 1, 0.2))
	invisible()
graphics.off()
}

book3.30 <- 
function()
{
trellis.device(postscript, file = "book3.30.ps", color = TRUE, paper = "legal")
print(
	xyplot(babinet~concentration^(1/3),
		data=polarization,
		aspect=1,
		sub = list("Figure 3.30",cex=.8),
		xlab = "Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Babinet Point (degrees)")
)
graphics.off()
}

book3.31 <- 
function()
{
set.seed(19)
trellis.device(postscript, file = "book3.31.ps", color = TRUE, paper = "legal")
print(
	xyplot(jitter(babinet)~jitter(concentration^(1/3)),
		data=polarization,
		aspect=1,
		sub = list("Figure 3.31",cex=.8),
		xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Jittered Babinet Point (degrees)")
)
graphics.off()
}

book3.32 <- 
function()
{

add.line <- trellis.par.get("add.line")
set.seed(19)
curve <- loess.smooth(concentration^(1/3), babinet, span = 3/4)
y <- jitter(babinet)
x <- jitter(concentration^(1/3))
trellis.device(postscript, file = "book3.32.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ x,
	       prepanel = function(x, y)
			list(dx = diff(curve$x), dy = diff(curve$y)),
	       panel = function(x, y){
			panel.xyplot(x, y)
			panel.lines(curve, lwd = add.line$lwd, lty = add.line$lty, col = add.line$col)
	       },
		sub = list("Figure 3.32",cex=.8),
	       xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
	       ylab = "Jittered Babinet Point (degrees)")
)
graphics.off()
}

book3.33 <- 
function()
{

conc <- concentration^(1/3)
set.seed(19)
trellis.device(postscript, file = "book3.33.ps", color = TRUE, paper = "legal")
print(
	xyplot(residuals(loess(babinet~conc,span=3/4,family="s",degree=1)) ~ jitter(conc),
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.loess(conc,y,span = 1/3)
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.33",cex=.8),
		xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Residual Babinet Point (degrees)")
)
graphics.off()

}

book3.34 <- 
function()
{

add.line <- trellis.par.get("add.line")
set.seed(19)
curve <- loess.smooth(concentration^(1/3), babinet, span = 1/6, evaluation=200)
y <- jitter(babinet)
x <- jitter(concentration^(1/3))
trellis.device(postscript, file = "book3.34.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ x,
	       prepanel = function(x, y)
			list(dx = diff(curve$x), dy = diff(curve$y)),
	       panel = substitute(function(x, y){
			panel.xyplot(x, y)
			panel.lines(curve, lwd = add.line$lwd, lty = add.line$lty, col = add.line$col)
	       }),
	       sub = list("Figure 3.34",cex=.8),
	       xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
	       ylab = "Jittered Babinet Point (degrees)")
)
graphics.off()

}

book3.35 <- 
function()
{
conc <- polarization$concentration^(1/3)
set.seed(19)
trellis.device(postscript, file = "book3.35.ps", color = TRUE, paper = "legal")
print(
	xyplot(
	residuals(loess(polarization$babinet ~ conc, span=1/6,family="s",degree=1)) ~ jitter(conc),
		panel = substitute(function(x, y) {
			panel.xyplot(x, y)
			panel.loess(conc, y, span = 1/3)
			panel.abline(h=0)
		}),
		aspect=1,
		sub = list("Figure 3.35",cex=.8),
		xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Residual Babinet Point (degrees)")
)
graphics.off()
}

book3.36 <- 
function()
{

set.seed(19)
curve <- loess.smooth(concentration^(1/3), babinet, span = 1/3)
y <- jitter(babinet)
x <- jitter(concentration^(1/3))
trellis.device(postscript, file = "book3.36.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ x,
	       prepanel = function(x, y)
			list(dx = diff(curve$x), dy = diff(curve$y)),
	       panel = substitute(function(x, y){
			add.line <- trellis.par.get("add.line")
			panel.xyplot(x, y)
			panel.lines(curve, lwd = add.line$lwd, lty = add.line$lty, col = add.line$col)
	       }),
		sub = list("Figure 3.36",cex=.8),
		aspect=1.1,
	       xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
	       ylab = "Babinet Point (degrees)")
)
graphics.off()

}

book3.37 <- 
function()
{
conc <- polarization$concentration^(1/3)
set.seed(19)
trellis.device(postscript, file = "book3.37.ps", color = TRUE, paper = "legal")
print(
	xyplot(
	residuals(loess(polarization$babinet ~ conc, span=1/3,family="s",degree=1))~jitter(conc),
		panel = substitute(function(x, y) {
			panel.xyplot(x, y)
			panel.loess(conc, y, span = 1/3)
			panel.abline(h=0)
		}),
		aspect=1,
		sub = list("Figure 3.37",cex=.8),
		xlab = "Jittered Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Residual Babinet Point (degrees)")
)
graphics.off()
}

book3.38 <- 
function()
{
cubiccon=polarization$concentration^(1/3)
trellis.device(postscript, file = "book3.38.ps", color = TRUE, paper = "legal")
print(
      qqmath(~residuals(loess(polarization$babinet~cubiccon,  span=1/3, family="s", degree=1)),
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect = 1,
		sub = list("Figure 3.38",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual Babinet Point (degrees)")
)
graphics.off()
}

book3.39 <- 
function()
{
cubiccon=polarization$concentration^(1/3)
trellis.device(postscript, file = "book3.39.ps", color = TRUE, paper = "legal")
print(
	rfs(loess(polarization$babinet~cubiccon,  span=1/3, family="s", degree=1),
		sub = list("Figure 3.39",cex=.8),
		aspect=1,
		ylab="Babinet Point (degrees)")
)
graphics.off()
}

		
book3.4 <- 
function()
{
fit <- lm(cp.ratio~area, data = ganglion)
trellis.device(postscript, file = "book3.4.ps", color = TRUE, paper = "legal")
print(
	xyplot(cp.ratio~area,
		data = ganglion,
		prepanel = prepanel.lmline,
		panel = substitute(function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(fit)
		}),
		aspect=1,
		sub = list("Figure 3.4",cex=.8),
		xlab = "Area (square mm)", ylab="CP Ratio")
)
graphics.off()
}

book3.40 <- 
function()
{
trellis.device(postscript, file = "book3.40.ps", color = TRUE, paper = "legal")
print(
	xyplot(babinet ~ concentration^(1/3),
		data = polarization,
		prepanel = function(x, y)
			prepanel.loess(x, y, span = 1/3),
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.abline(v = c(3.936497, 4.188859))
		},
		sub = list("Figure 3.40",cex=.8),
		aspect=1.1,
		xlab = "Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Babinet Point (degrees)")
)
graphics.off()
}

book3.41 <- 
function()
{
trellis.device(postscript, file = "book3.41.ps", color = TRUE, paper = "legal")
print(
	xyplot(babinet ~ concentration^(1/3),
		data = polarization,
		prepanel = function(x, y)
			prepanel.loess(x, y, span = 1/3),
		panel = function(x, y) {
			panel.xyplot(x, y)	
			panel.loess(x, y, span = 1/3)
			panel.abline(v = c(3.936497, 4.188859))
		},
		sub = list("Figure 3.41",cex=.8),
		aspect=1.1,
		xlab = "Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Babinet Point (degrees)")
)
graphics.off()
}

book3.42 <- 
function()
{
cubiccon=polarization$concentration^(1/3)
trellis.device(postscript, file = "book3.42.ps", color = TRUE, paper = "legal")
print(
	xyplot(residuals(loess(polarization$babinet ~ cubiccon, span = 1/3, 
			family = "s", degree = 1)) ~ cubiccon,
		panel = function(x, y) {
	                panel.xyplot(x, y)
	                panel.abline(h = 0, v = c(3.936497, 4.188859))
		},
		sub = list("Figure 3.42",cex=.8),
		aspect=1,
		xlab = "Cube Root Concentration (cube root micrograms per meter)",
		ylab = "Residual Babinet Point (degrees)")
)
graphics.off()		
}

book3.43 <- 
function()
{
trellis.device(postscript, file = "book3.43.ps", color = TRUE, paper = "legal")
print(
	plot(equal.count(polarization$concentration^(1/3), 15, .5),
		aspect = 1,
		sub = list("Figure 3.43",cex=.8),
		xlab = "Cube Root Concentration (cube root micrograms per meter)")
)
graphics.off()
}

book3.44 <- 
function()
{

cube.root <- concentration^(1/3)
pol.m <- loess(babinet~cube.root,span=1/3,family="s",degree=1)
trellis.device(postscript, file = "book3.44.ps", color = TRUE, paper = "legal")
print(
	bwplot(equal.count(cube.root,15,1/2) ~ residuals(pol.m), 
		panel = function(x, y) {
			panel.bwplot(x, y)
			panel.abline(v=0)
		},
		aspect=1,
		sub = list("Figure 3.44",cex=.8),
		xlab="Residual Babinet Point (degrees)")
)
graphics.off()

}

book3.45 <- 
function()
{
trellis.device(postscript, file = "book3.45.ps", color = TRUE, paper = "legal")
print(
	xyplot(facet~temperature,
		data = fly,
		aspect=1,
		sub = list("Figure 3.45",cex=.8),
		xlab="Temperature (Degrees Celsius)",
		ylab="Facet Number")
)
graphics.off()
}

book3.46 <- 
function()
{
set.seed(19)
trellis.device(postscript, file = "book3.46.ps", color = TRUE, paper = "legal")
print(
	xyplot(jitter(facet,1)~jitter(temperature,2),
		data = fly,
		aspect=1,
		sub = list("Figure 3.46",cex=.8),
		xlab="Jittered Temperature (Degrees Celsius)",
		ylab="Jittered Facet Number")
)
graphics.off()
}

book3.47 <- 
function()
{
trellis.device(postscript, file = "book3.47.ps", color = TRUE, paper = "legal")
print(
	bwplot(factor(temperature) ~ facet,
		data=fly,
		aspect=1,
		sub = list("Figure 3.47",cex=.8),
		xlab = "Facet Number",
		ylab = "Temperature (Degrees Celsius)")
)
graphics.off()
}

book3.48 <- 
function()
{
fly.m <- oneway(facet~temperature,data=fly,spread=1)
Temperature <- factor(fly$temperature)
trellis.device(postscript, file = "book3.48.ps", color = TRUE, paper = "legal")
print(
	qqmath(~residuals(fly.m)|Temperature,
		prepanel = prepanel.qqmathline,
		panel = function(x, ...) {
			panel.grid()
			panel.qqmathline(x, distribution = qnorm)
			panel.qqmath(x, ...)
		},
		aspect = 1,
		sub = list("Figure 3.48",cex=.8),
		xlab="Unit Normal Quantile",
		ylab="Residual Facet Number")
)
graphics.off()
}

book3.49 <- 
function()
{

fly.lm <- lm(facet~temperature)
ylim <- range(fitted.values(fly.lm))
fly.means <- tapply(facet,temperature,mean)
fly.temps <- unique(temperature)
trellis.device(postscript, file = "book3.49.ps", color = TRUE, paper = "legal")
print(
	xyplot(fly.means~fly.temps,
	        ylim = ylim,
		panel = substitute(function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(fly.lm)
		}),
		aspect=1,
		sub = list("Figure 3.49",cex=.8),
		xlab="Temperature (Degrees Celsius)",
		ylab="Facet Number")
)
graphics.off()

}

book3.5 <- 
function()
{

add.line <- trellis.par.get("add.line")
gan.lm <- lm(cp.ratio~area+I(area^2))
gan.x <- seq(min(area),max(area),length=50)
gan.fit <- gan.lm$coef[1]+gan.lm$coef[2]*gan.x+gan.lm$coef[3]*gan.x^2
trellis.device(postscript, file = "book3.5.ps", color = TRUE, paper = "legal")
print(
	xyplot(cp.ratio~area,
		prepanel = substitute(function(x, y) 
			list(dx = diff(gan.x), dy = diff(gan.fit))),
		panel = substitute(function(x, y) {
	                panel.xyplot(x, y)
			panel.lines(gan.x, gan.fit, lwd = add.line$lwd, lty =add.line$lty, col = add.line$col)
		}),
		sub = list("Figure 3.5",cex=.8),
		xlab = "Area (square mm)",
		ylab = "CP Ratio")
)
graphics.off()

}

book3.50 <- 
function()
{

fly.lm <- lm(facet~temperature)
fly.means <- tapply(facet,temperature,mean)
fly.temps <- unique(temperature)
fly.res <- fly.means-predict(fly.lm,data.frame(temperature=fly.temps))
trellis.device(postscript, file = "book3.50.ps", color = TRUE, paper = "legal")
print(
	xyplot(fly.res~fly.temps,
		panel = function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(h=0)
		},
		aspect=1,
		sub = list("Figure 3.50",cex=.8),
		xlab="Temperature",
		ylab="Residual Facet Number")
)
graphics.off()

}


book3.51 <- 
function()
{
trellis.device(postscript, file = "book3.51.ps", color = TRUE, paper = "legal")
on.exit(par(oldpar))
oldpar <- par("pty", "xaxt", "yaxt")
par(pty = "s", xaxt = "n", yaxt = "n")

x <- c(seq(22, 1, length = 11), seq(22, 1, length = 11))
y <- c(rep(1, 11), rep(2, 11))
symbols(x, y, circles = diameter, inches = 1/6, ylim = c(0.5, 2.2))
y <- c(rep(0.9, 11), rep(1.8, 11))
mtext("Figure 3.51",1,1,cex=.8)
text(x, y, labels=dimnames(playfair)[[1]], srt = 45, adj = 1)

invisible()
graphics.off()
}

book3.52 <- 
function()
{
# The first two lines is not necessary if Playfair is a matrix;
# in which case dotplot(Playfair[,"population"], ...) will work.
	population <- playfair$population	
	names(population) <- row.names(playfair)
trellis.device(postscript, file = "book3.52.ps", color = TRUE, paper = "legal")
print(
	dotplot(population,
		aspect = 1,
		xlim = range(0, population),
		sub = list("Figure 3.52",cex=.8),
		xlab = "Population (Thousands)")
)
graphics.off()
}

book3.53 <- 
function()
{
trellis.device(postscript, file = "book3.53.ps", color = TRUE, paper = "legal")
print(
	xyplot(diameter ~ sqrt(population),
		data = playfair,
		aspect=1,
		ylab = "Diameter (mm)",
		sub = list("Figure 3.53",cex=.8),
		xlab = "Square Root Population (square root thousands)")
)
graphics.off()
}

book3.54 <- 
function()
{
pla.lm <- lm(diameter ~ sqrt(population) - 1, data = playfair)
trellis.device(postscript, file = "book3.54.ps", color = TRUE, paper = "legal")
print(
	xyplot(diameter ~ sqrt(population),
		data = playfair,
		ylim = range(playfair$diameter, fitted.values(pla.lm)),
		panel = substitute(function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(pla.lm)
		}),
		aspect=1,
		ylab = "Diameter (mm)",
		sub = list("Figure 3.54",cex=.8),
		xlab = "Square Root Population (square root thousands)")
)
graphics.off()
}

book3.55 <- 
function()
{
trellis.device(postscript, file = "book3.55.ps", color = TRUE, paper = "legal")
print(
	xyplot(
	residuals(lm(diameter ~ sqrt(population) - 1, data = playfair)) ~ sqrt(playfair$population),
		panel = function(x, y) {
	                panel.xyplot(x, y)
			panel.abline(h=0)
		},
		aspect = "xy",
		sub = list("Figure 3.55",cex=.8),
		xlab = "Square Root Population (square root thousands)",
		ylab = "Residual Diameter (mm)")
)
graphics.off()
}

book3.56 <- 
function()
{
set.seed(19)
trellis.device(postscript, file = "book3.56.ps", color = TRUE, paper = "legal")
print(
	xyplot(jitter(wind)~jitter(temperature),
		data = environmental,
		aspect=1,
		sub = list("Figure 3.56",cex=.8),
		xlab="Jittered Temperature (Degrees Fahrenheit)",
		ylab="Jittered Wind Speed (mph)")
)
graphics.off()
}

book3.57 <- 
function()
{
trellis.device(postscript, file = "book3.57.ps", color = TRUE, paper = "legal")
print(
	qqmath(~environmental$temperature,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect = 1,
		sub = list("Figure 3.57",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Temperature (Degrees Fahrenheit)")
)
graphics.off()
}

book3.58 <- 
function()
{
trellis.device(postscript, file = "book3.58.ps", color = TRUE, paper = "legal")
print(
	qqmath(~environmental$wind,
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect = 1,
		sub = list("Figure 3.58",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Wind Speed (mph)")
)
graphics.off()
}

book3.59 <- 
function()
{
trellis.device(postscript, file = "book3.59.ps", color = TRUE, paper = "legal")
print(
	xyplot(wind~temperature,
		data = environmental,
		panel=function(x,y){
			add.line <- trellis.par.get("add.line")
			panel.xyplot(x,y)
			panel.lines(loess.smooth(x, y, span=3/4, family="g"), 
				lwd = add.line$lwd, lty = add.line$lty,
				col = add.line$col)
			other.fit <- loess.smooth(y, x, span=3/4, family="g")
			panel.lines(other.fit$y,other.fit$x, 
				lwd = add.line$lwd, lty = add.line$lty,
				col = add.line$col)
		},
		aspect = 1,
		sub = list("Figure 3.59",cex=.8),
		xlab = "Temperature (Degrees Fahrenheit)",
		ylab = "Wind Speed (mph)")
)
graphics.off()
}

book3.6 <- 
function()
{
x <- seq(5, 10, length=150)
z <- seq(0, 1, length=50)
y <- 5 + c(-(1-z)^.5, sin(z*2*pi), .5*z^2)
fit <- loess.smooth(x, y, degree = 2, family = "g", span = .5, evaluation = 150)
set.seed(20)
y <- fit$y + rnorm(150, 0, .2)
x <- fit$x
trellis.device(postscript, file = "book3.6.ps", color = TRUE, paper = "legal")
print(
	xyplot(y ~ x, 
		prepanel = function(x, y) 
			prepanel.loess(x, y, degree = 2, family = "g", span = .3, evaluation = 300), 
		panel = function(x, y){
	        	panel.xyplot(x, y)
			panel.loess(x, y, degree = 2, family = "g", span = .3, evaluation = 300)
		},
		sub = list("Figure 3.6",cex=.8),
		xlab = "x", 	
		ylab = "y")
)
graphics.off()
}

book3.60 <- 
function()
{
trellis.device(postscript, file = "book3.60.ps", color = TRUE, paper = "legal")
print( 
	xyplot(stamford~yonkers,
		data = ozone,
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.abline(0, 1)
		},
		aspect=1,
		scale = list(limits = range(ozone)),
		sub = list("Figure 3.60",cex=.8),
		xlab="Yonkers Concentration (ppb)",
		ylab="Stamford Concentration (ppb)")
)
graphics.off()
}

book3.61 <- 
function()
{
diff=ozone$stamford-ozone$yonkers
mean=(ozone$yonkers+ozone$stamford)/2
trellis.device(postscript, file = "book3.61.ps", color = TRUE, paper = "legal")
print(
	tmd(xyplot(diff~mean	),
		prepanel = function(x, y)
			prepanel.loess(x, y,span=1,degree=1,family="g"),
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.loess(x, y, span=1,degree=1,family="g")
			panel.abline(h = 0)
		},
		aspect = "xy",
		ylab="Difference (ppb)",
		sub = list("Figure 3.61",cex=.8),
		xlab="Mean (ppb)")
)
graphics.off()
}

book3.62 <- 
function()
{
trellis.device(postscript, file = "book3.62.ps", color = TRUE, paper = "legal")
print(
	xyplot(log(stamford, 2) ~ log(yonkers, 2),
		data = ozone,
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.abline(0, 1)
		},
		scale = list(limits = range(log(ozone, 2))),
		aspect = 1,
		sub = list("Figure 3.62",cex=.8),
		xlab="Log Yonkers Concentration (log base 2 ppb)",
		ylab="Log Stamford Concentration (log base 2 ppb)")
)
graphics.off()
}

book3.63 <- 
function()
{
diff=log(ozone$stamford,2)-log(ozone$yonkers,2)
mean=(log(ozone$stamford,2)+log(ozone$yonkers,2))/2
trellis.device(postscript, file = "book3.63.ps", color = TRUE, paper = "legal")
print(
	tmd(xyplot(diff~mean),
		prepanel = function(x, y)
			prepanel.loess(x, y,span=1,degree=1,family="g"),
		panel = function(x, y) {
			panel.xyplot(x, y)
			panel.loess(x, y, span=1,degree=1,family="g")
			panel.abline(h = 0)
		},
		aspect=1.5,
		scale = list(y = list(labels = c("-1", "0", "1"),
			at = c(-1, 0, 1))),
		ylab="Log Difference (log base 2 ppb)",
		sub = list("Figure 3.63",cex=.8),
		xlab="Log Mean (log base 2 ppb)")
)
graphics.off()
}

book3.64 <- 
function()
{
trellis.device(postscript, file = "book3.64.ps", color = TRUE, paper = "legal")
print(
	xyplot(incidence ~ year,
		data = melanoma,
		panel = function(x, y)
			panel.xyplot(x, y, type="o", pch = 16, cex = .75),
		aspect = "xy",
		ylim = c(0, 5),
		sub = list("Figure 3.64",cex=.8),
		xlab = "Year",
		ylab = "Incidence")
)
graphics.off()
}

book3.65 <- 
function()
{

mel.low <- stl2new(ts(incidence, start = 1936, frequency=1, end = 1972), 
		s.window = 30,
		s.degree = 1)
trellis.device(postscript, file = "book3.65.ps", color = TRUE, paper = "legal")
print(
	xyplot(mel.low$fc.30~year,
		aspect="xy",
		type="l",
		ylab = "Incidence",
		sub = list("Figure 3.65",cex=.8),
		xlab = "Year")
)
graphics.off()

}

book3.66 <- 
function()
{

tsnow=ts(incidence, start = 1936, frequency=1, end = 1972)
mel.m <- stl(tsnow, s.window = 30, s.degree = 1)
trellis.device(postscript, file = "book3.66.ps", color = TRUE, paper = "legal")
print(          
	xyplot(mel.m$time.series[,"remainder"] ~ year,
		panel = function(x, y) {
			panel.xyplot(x, y, type = "b", pch = 16, cex=.6)
			panel.abline(h = 0)
		},
		aspect = "xy",
		ylab = "Incidence",
		sub = list("Figure 3.66",cex=.8),
		xlab = "Year")
)
graphics.off()

}

book3.67 <- 
function()
{

mel.mlow <- stl(ts(incidence, start = 1936, frequency=1, end = 1972), 
		l.window = 30,
		l.degree = 1)
mel.mo <- stl(mel.mlow$remainder,
		l.window = 9,
		l.degree = 2)
trellis.device(postscript, file = "book3.67.ps", color = TRUE, paper = "legal")
print(
	xyplot(mel.mo$fc.9~year,
		panel = function(x, y) {
			panel.xyplot(x, y, type = "o", pch = 16, cex = .6)
			panel.abline(h = 0)
		},
		aspect = "xy",
		ylab = "Incidence",
		sub = list("Figure 3.67",cex=.8),
		xlab = "Year")
)
graphics.off()

}

book3.68 <- 
function()
{

mel.mboth <- stl(ts(incidence, start = 1936, frequency=1, end = 1972),
		l.window = c(30,9),
		l.degree = c(1,2))
ylim <- c(-1,1)*2*max(abs(mel.mboth$remainder))
trellis.device(postscript, file = "book3.68.ps", color = TRUE, paper = "legal")
print(
	ans <- xyplot(mel.mboth$remainder~year,
		panel = function(x, y) {
			panel.xyplot(x, y, type="o", pch = 16, cex=.75)
			panel.abline(h=0)
		},
		aspect="xy",
		las = 0, # vertical y-axis labels
		ylim=ylim,
		ylab = "Residual\nIncidence",
		sub = list("Figure 3.68",cex=.8),
		xlab = "Year")
)
graphics.off()

}

book3.69 <- 
function()
{

mel.mboth <- stl(ts(incidence, start = 1936, frequency=1, end = 1972),
		l.window = c(30,9),
		l.degree = c(1,2))
trellis.device(postscript, file = "book3.69.ps", color = TRUE, paper = "legal")
print(
	qqmath(~mel.mboth$remainder,
		prepanel = prepanel.qqmathline,
	        panel = function(x, y) {
	                panel.qqmathline(y, distribution = qnorm)
	                panel.qqmath(x, y)
	        },
		aspect=1,
		sub = list("Figure 3.69",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual Incidence")
)
graphics.off()

}

book3.70 <- 
function()
{

mel.mboth <- stl(ts(incidence, start = 1936, frequency=1, end = 1972),
		l.window = c(30,9),
		l.degree = c(1,2))
n <- length(incidence)
mel.mboth$fc.30 <- mel.mboth$fc.30 - mean(mel.mboth$fc.30)
mel.components <- c(mel.mboth$fc.30,mel.mboth$fc.9, mel.mboth$remainder)
mel.year <-  rep(1936:1972,3)
mel.names <-
ordered(rep(c("Trend","Oscillatory","Residuals"),c(n,n,n)),
	 c("Trend","Oscillatory","Residuals"))
trellis.device(postscript, file = "book3.70.ps", color = TRUE, paper = "legal")
print(
	xyplot(mel.components~mel.year|mel.names,
	        panel = function(x, y) {
	                panel.grid(h = 7)
	                panel.xyplot(x, y, type = "o", pch = 16, cex = .7)
	        },
		layout=c(3,1),
		ylim = c(-1,1)*max(abs(mel.components)),
		sub = list("Figure 3.70",cex=.8),
		xlab="Year",
		ylab="Incidence")
)
graphics.off()

}


book3.71 <- 
function()
{
trellis.device(postscript, file = "book3.71.ps", color = TRUE, paper = "legal")
	print(xyplot(stl(ts(incidence, start = 1936, frequency=1, end = 1972),
		        fc.window = c(30,9), fc.degree = c(1,2))$fc.9 ~ year,
			data = melanoma,
			panel=function(x, y) {
				panel.grid(v = 5,h = 0)
				panel.xyplot(x, y, type = "l")
			},
			aspect = "xy",
			ylab = "Incidence",
			xlab = "Year"),	
			position = c(0, .4, 1, 1), more = T)
	print(xyplot(sunspot ~ year,
			data = melanoma,
			panel=function(x, y) {
				panel.grid(v = 5, h = 0)
				panel.xyplot(x, y, type = "l")
			},
			aspect = "xy",
			ylab = "Sunspot Number",
			sub = list("Figure 3.71",cex=.8),
			xlab = "Year"),
			position = c(0, .2, 1, .6))
	invisible()
graphics.off()
}

book3.72 <- 
function()
{
trellis.device(postscript, file = "book3.72.ps", color = TRUE, paper = "legal")
print(
	xyplot(carbon.dioxide~time(carbon.dioxide),
		aspect="xy",
		type="l",
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		ylab = "CO2 (ppm)",
		sub = list("Figure 3.72",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.73 <- 
function()
{
trellis.device(postscript, file = "book3.73.ps", color = TRUE, paper = "legal")
print(
	xyplot(carbon.dioxide ~ time(carbon.dioxide),
		aspect = 1,
		type = "l",
		ylab = "Carbon Dioxide (ppm)",
		sub = list("Figure 3.73",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.74 <- 
function()
{
car.sfit <- as.data.frame(stl(carbon.dioxide, s.window = 25, s.degree=1)$time.series)$seasonal
ylim <- c(-1,1)*max(abs(car.sfit))
car.time <- time(carbon.dioxide)
trellis.device(postscript, file = "book3.74.ps", color = TRUE, paper = "legal")
print(
	xyplot(car.sfit ~ car.time | equal.count(car.time, 7, overlap=0),
	        panel = function(x, y) {
	                panel.grid(v = 0)
	                panel.xyplot(x, y, type="l")
	        },
		aspect = "xy",
		strip = F,
		ylim = ylim,
		scale = list(x = "free"),
		layout = c(1, 7),
		sub = list("Figure 3.74",cex=.8),
		xlab = "Year",
		ylab = "Carbon Dioxide (ppm)")
)
graphics.off()
}

book3.75 <- 
function()
{
car.sfit <- as.data.frame(stl(carbon.dioxide, s.window = 25, s.degree=1)$time.series)$seasonal
car.time <- time(carbon.dioxide)-1959
car.subseries <- factor(cycle(carbon.dioxide),label=month.abb)
trellis.device(postscript, file = "book3.75.ps", color = TRUE, paper = "legal")
print(
	xyplot(car.sfit~car.time|car.subseries,
		layout=c(12,1),
		panel=function(x,y) {
			panel.xyplot(x,y,type="l")
			panel.abline(h=mean(y))
		},
#		aspect="xy",
		scales=list(cex=.7, tick.number=3),
		sub = list("Figure 3.75",cex=.8),
		xlab="",
		ylab="Carbon Dioxide (ppm)")
)
graphics.off()
}

book3.76 <- 
function()
{
deseasonal <- as.data.frame(stl(carbon.dioxide,s.window = 25, s.degree=1)$time.series)$remainder + 
		as.data.frame(stl(carbon.dioxide,s.window = 25, s.degree=1)$time.series)$trend
trellis.device(postscript, file = "book3.76.ps", color = TRUE, paper = "legal")
print(
	xyplot(deseasonal ~ time(carbon.dioxide),
		aspect = 1,
		type = "l",
		ylab = "Carbon Dioxide (ppm)",
		sub = list("Figure 3.76",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.77 <- 
function()
{
library(stl2)
trellis.device(postscript, file = "book3.77.ps", color = TRUE, paper = "legal")
print(
	xyplot(stl2(carbon.dioxide,s.window = 25, s.degree=1, fc.window=101, 
			fc.degree=1)$fc$fc.101 ~ time(carbon.dioxide),
		aspect="xy",
		type="l",
		ylab = "Carbon Dioxide (ppm)",
		sub = list("Figure 3.77",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.78 <- 
function()
{
library(stl2)
trellis.device(postscript, file = "book3.78.ps", color = TRUE, paper = "legal")
print(
	xyplot(stl2(carbon.dioxide,s.window = 25, s.degree=1, fc.window=101, 
			fc.degree=1)$fc$remainder ~ time(carbon.dioxide),
		panel = function(x, y) {
			panel.xyplot(x, y, type = "l")
			panel.abline(h = 0)
		},
		aspect = 1/3,
		ylab = "Carbon Dioxide (ppm)",
		sub = list("Figure 3.78",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.79 <- 
function()
{
library(stl2)
trellis.device(postscript, file = "book3.79.ps", color = TRUE, paper = "legal")
print(
	xyplot(stl2(carbon.dioxide,s.window = 25, s.degree=1, fc.window=c(101,35), 
			fc.degree=c(1,2))$fc$fc.35 ~ time(carbon.dioxide),
		panel = function(x, y) {
			panel.xyplot(x, y, type="l")
			panel.abline(h = 0)
		},
		aspect = "xy",
		ylab = "CO2 (ppm)",
		sub = list("Figure 3.79",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.80 <- 
function()
{
library(stl2)
trellis.device(postscript, file = "book3.80.ps", color = TRUE, paper = "legal")
print(
	xyplot(stl2(carbon.dioxide,s.window = 25, s.degree=1, fc.window=c(101,35), 
			fc.degree=c(1,2))$data$remainder ~ time(carbon.dioxide),
		panel = function(x, y, ...) {
			panel.xyplot(x, y, type = "l")
			panel.abline(h = 0)
		},
		aspect = .1,
		ylab = "Residual CO2 (ppm)",
		sub = list("Figure 3.80",cex=.8),
		xlab = "Year")
)
graphics.off()
}

book3.81 <- 
function()
{
trellis.device(postscript, file = "book3.81.ps", color = TRUE, paper = "legal")
print(
	qqmath( ~ stl2(carbon.dioxide,s.window = 25, s.degree=1,
	                fc.window=c(101,35),fc.degree=c(1,2))$data$remainder,
		prepanel = prepanel.qqmathline,
	        panel = function(x,...) {
	                panel.qqmathline(x,...)
	                panel.qqmath(x,...)
	        },
		aspect = 1,
		sub = list("Figure 3.81",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual Carbon Dioxide (ppm)")
)
graphics.off()
}

book3.82 <- 
function()
{
car.fit <- stl2(carbon.dioxide,s.window = 25, s.degree=1,
		fc.window=c(101,35), fc.degree=c(1,2))
co2 <- c(car.fit$fc$fc.35,car.fit$data$seasonal,car.fit$data$remainder)
n <- length(carbon.dioxide)
co2.names <- factor(rep(c("3.5-Year","Seasonal","Residual"),rep(n,3)))
co2.names <- ordered(co2.names, c("3.5-Year","Seasonal","Residual"))
co2.year <- rep(time(carbon.dioxide),3)
trellis.device(postscript, file = "book3.82.ps", color = TRUE, paper = "legal")
print(
	xyplot(co2 ~ co2.year | co2.names,
	        panel = function(x, y) {
	                panel.grid(h = 5, v = 0)
	                panel.xyplot(x, y, type = "l")},
		layout = c(1, 3),
		ylim = c(-1,1) *1.05 * max(abs(co2)),
		sub = list("Figure 3.82",cex=.8),
		xlab = "Year",
		ylab = "Carbon Dioxide (ppm)")
)
graphics.off()
}

book3.9 <- 
function()
{
X <- seq(5, 10, length=150)
z <- seq(0, 1, length=50)
base.y <- 5 + c(-(1-z)^.5, sin(z*2*pi), .5*z^2)
true.y <- loess.smooth(X, base.y, degree = 2, family = "g", span = .5, evaluation = 150)$y
set.seed(20)
Y <- true.y+rnorm(150, 0, .2)
fit1 <- loess.smooth(X, Y, degree = 2, family = "g", span = .1, evaluation = 150)
fit2 <- loess.smooth(X, Y, degree = 2, family = "g", span = .3, evaluation = 150)
fit3 <- loess.smooth(X, Y, degree = 2, family = "g", span = .6, evaluation = 150)
alpha <- factor(rep(c(.1,.3,.6),rep(length(X),3)))
fits <- c(fit1$y, fit2$y, fit3$y)
X <- rep(X,3)
Y <- rep(Y,3)
trellis.device(postscript, file = "book3.9.ps", color = TRUE, paper = "legal")
print(
	xyplot(Y ~ X | alpha, 
		prepanel = substitute(function(x, y, subscripts) {
			list(dx = diff(x), dy=diff(fits[subscripts]))}),
		panel = substitute(function(x, y, subscripts){
			add.line <- trellis.par.get("add.line")
			panel.xyplot(x, y)
			panel.lines(x, fits[subscripts], lwd = add.line$lwd,
				lty = add.line$lty, col = add.line$col)
		}),
		strip = function(...) strip.default(..., strip.names = T),
		layout = c(1, 3),
		sub = list("Figure 3.9",cex=.8),
		xlab = "x", 	
		ylab = "y")
)
graphics.off()
}








book4.1 <- 
function()
{
trellis.device(postscript, file = "book4.1.ps", color = TRUE, paper = "legal")
print(
	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))
)
graphics.off()
}
	
book4.10 <- 
function()
{
trellis.device(postscript, file = "book4.10.ps", color = TRUE, paper = "legal")
print(
	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)"))
)
graphics.off()
}

book4.11 <- 
function()
{
trellis.device(postscript, file = "book4.11.ps", color = TRUE, paper = "legal")
print(
	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)"))
)
graphics.off()
}

book4.12 <- 
function()
{

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
trellis.device(postscript, file = "book4.12.ps", color = TRUE, paper = "legal")
	ans <- xyplot(eth.fit ~ CE.grid$C | EQ.Ratio,
#		aspect=2,
		strip = FALSE,
		layout = c(8, 2),
		panel = function(x, y) {
			panel.grid(v=2)
			panel.xyplot(x, y, type = "l")
		},
		xlab = "Compression Ratio",
                sub = list("Figure 4.12",cex=.8),
		ylab = "NOx (micrograms/J)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	ans <- plot(shingle(seq(min(ethanol$E), max(ethanol$E), length=16)), 
		scales=list(x=list(cex=.7),y=list(cex=.5)),
		xlab = "Equivalence Ratio")
	print(ans, position=c(0,.6,1,1))
	invisible()
graphics.off()
}

book4.13 <- 
function()
{

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
trellis.device(postscript, file = "book4.13.ps", color = TRUE, paper = "legal")
	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",
		strip = FALSE, 
		layout = c(8, 2),
		xlab = "Equivalence Ratio",                
		sub = list("Figure 4.13",cex=.8),
		ylab = "NOx (micrograms/J)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	ans <- plot(shingle(seq(min(ethanol$C), max(ethanol$C), length=16)), 
		scales=list(x=list(cex=.7),y=list(cex=.5)),
		xlab = "Compression Ratio")
	print(ans, position=c(0,.6,1,1))
	invisible()
graphics.off()
}

book4.14 <- 
function()
{

c.values <- seq(min(C), max(C), length = 16)
e.values <- seq(min(E), max(E), length = 16)
trellis.device(postscript, file = "book4.14.ps", color = TRUE, paper = "legal")
print(
	xyplot(E ~ C,
		panel = substitute(function(x, y){
			panel.segments(rep(min(x),16),e.values,rep(max(x),16),e.values)
			panel.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")
)
graphics.off()

}

book4.15 <- 
function()
{
trellis.device(postscript, file = "book4.15.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()
}

book4.16 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
    	U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
	w[U > 1.] <- 0.
	w
}		  

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
trellis.device(postscript, file = "book4.16.ps", color = TRUE, paper = "legal")
	ans <- xyplot(rub.fit ~ grid$tensile.strength | Hardness,
		aspect="xy",
		layout = c(4, 1),
                sub = list("Figure 4.16",cex=.8),
		strip = FALSE,
		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)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	print(plot(shingle(seq(55, 83, length = 4)), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Hardness (degrees Shore)"), 
		position=c(0,.6,1,1))
	invisible()
graphics.off()
}

book4.17 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
    U <- abs(u/c)
    w <- ((1. + U) * (1. - U))^2.
    w[U > 1.] <- 0.
    w
}	  

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
trellis.device(postscript, file = "book4.17.ps", color = TRUE, paper = "legal")
	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")
		},
		strip = FALSE,
                sub = list("Figure 4.17",cex=.8),
		xlab = "Hardness (degrees Shore)",
		ylab = "Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)
	print(plot(shingle(seq(144,180,length = 3)), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Tensile Strength (kg/square cm)"),
		position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.18 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
    U <- abs(u/c)
    w <- ((1. + U) * (1. - U))^2.
    w[U > 1.] <- 0.
    w
}	  

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)
}
trellis.device(postscript, file = "book4.18.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()

}

book4.19 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
    U <- abs(u/c)
    w <- ((1. + U) * (1. - U))^2.
    w[U > 1.] <- 0.
    w 
}

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)
}
trellis.device(postscript, file = "book4.19.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()

}

book4.2 <- 
function()
{
trellis.device(postscript, file = "book4.2.ps", color = TRUE, paper = "legal")
print(
	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)"))
)
graphics.off()
}

book4.20 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
	w[U > 1.] <- 0.
        w
}

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)
trellis.device(postscript, file = "book4.20.ps", color = TRUE, paper = "legal")
	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)
		},
                sub = list("Figure 4.20",cex=.8),
		aspect = 2,
		strip = FALSE,
		layout = c(6, 1),
		xlab = "Tensile Strength (kg/square cm)", 
		ylab = "Residual Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)

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

book4.21 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}		   

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)
trellis.device(postscript, file = "book4.21.ps", color = TRUE, paper = "legal")
	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)
		},
		strip = FALSE,
		aspect = 2,
		layout = c(6, 1),
                sub = list("Figure 4.21",cex=.8),
		xlab = "Hardness (degrees Shore)", 
		ylab = "Residual Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,0.6), more=T)

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

book4.22 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w	  
}	   
new.rubber <- rubber[rubber$tensile.strength>130,]

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)
trellis.device(postscript, file = "book4.22.ps", color = TRUE, paper = "legal")
	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,
		strip = FALSE,
                sub = list("Figure 4.22",cex=.8),
		layout = c(6, 1),
	        xlab = "Tensile Strength (kg/square cm)",
	        ylab = "Residual Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,0.6), more=T)

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

book4.23 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}		   
new.rubber <- rubber[rubber$tensile.strength>130,]

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)
trellis.device(postscript, file = "book4.23.ps", color = TRUE, paper = "legal")
	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,
		strip = FALSE, 
		layout = c(6, 1),
                sub = list("Figure 4.23",cex=.8),
		xlab = "Hardness (degrees Shore)",
		ylab = "Residual Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,0.6), more=T)

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

book4.24 <- 
function()
{
 wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}	   
new.rubber <- rubber[rubber$tensile.strength>130,]

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
trellis.device(postscript, file = "book4.24.ps", color = TRUE, paper = "legal")
	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")
	        },
		strip = FALSE, 
                sub = list("Figure 4.24",cex=.8),
		xlab = "Tensile Strength (kg/square cm)",
		ylab = "Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)
        print(plot(shingle(seq(55, 83, length = 4)), 		
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Hardness (degrees Shore)"),
		position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.25 <- 
function()
{
 wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}	   
new.rubber <- rubber[rubber$tensile.strength>130,]

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
trellis.device(postscript, file = "book4.25.ps", color = TRUE, paper = "legal")
	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")
	        },
		strip = FALSE,
                scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Hardness (degrees Shore)",
		ylab = "Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)
	print(plot(shingle(seq(144,180,length = 3)), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Tensile Strength (kg/square cm)"),
                position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.26 <- 
function()
{
 wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}		   
new.rubber <- rubber[rubber$tensile.strength>130,]

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)
}
trellis.device(postscript, file = "book4.26.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()

}

book4.27 <- 
function()
{
 wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}	   
new.rubber <- rubber[rubber$tensile.strength>130,]

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)
}
trellis.device(postscript, file = "book4.27.ps", color = TRUE, paper = "legal")
print(
	qqmath(~residuals(rub.lm), 
		prepanel = prepanel.qqmathline,
	        panel = function(x, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		sub = list("Figure 4.27",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab="Residual Abrasion Loss (g/hp-hour)")
)
graphics.off()

}

book4.28 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}
new.rubber <- rubber[rubber$tensile.strength>130,]

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)
}
trellis.device(postscript, file = "book4.28.ps", color = TRUE, paper = "legal")
print(
	ans <- rfs(rub.lm,
		sub = list("Figure 4.28",cex=.8),
		aspect=1.5,
		ylab="Residual Abrasion Loss (g/hp-hour)")
)
graphics.off()

}

book4.29 <- 
function()
{
trellis.device(postscript, file = "book4.29.ps", color = TRUE, paper = "legal")
print(          
	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)")
)
graphics.off()
}

book4.3 <- 
function()
{

Hardness <- equal.count(hardness,6,3/4)
trellis.device(postscript, file = "book4.3.ps", color = TRUE, paper = "legal")
	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(6,1),
		strip = FALSE, 
                sub = list("Figure 4.3",cex=.8),
		xlab = "Tensile Strength (kg/square cm)",
		ylab = "Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)

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

book4.30 <- 
function()
{

eth.m <- loess(NOx ~ C * E, span = 1/3, degree = 2,
	parametric = "C", drop.square = "C",family="s")
trellis.device(postscript, file = "book4.30.ps", color = TRUE, paper = "legal")
print(
	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))
)
graphics.off()

}

book4.31 <- 
function()
{
trellis.device(postscript, file = "book4.31.ps", color = TRUE, paper = "legal")
print(
	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, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1, 
		sub = list("Figure 4.31",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Residual NOx (micrograms/J)")
)
graphics.off()
}

book4.32 <- 
function()
{
trellis.device(postscript, file = "book4.32.ps", color = TRUE, paper = "legal")
print(
	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=1,
		ylab = "NOx (micrograms/J)")
)
graphics.off()
}

book4.33 <- 
function()
{

eth.m <- loess(log(NOx,2) ~ C * E, span = 1/3, degree = 2,
	parametric = "C", drop.square = "C",family="s")
trellis.device(postscript, file = "book4.33.ps", color = TRUE, paper = "legal")
print(
	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))
)
graphics.off()

}

book4.35 <- 
function()
{
set.seed(19)
new.slope <- slope
new.slope$percent <- jitter(new.slope$percent,2)
trellis.device(postscript, file = "book4.35.ps", color = TRUE, paper = "legal")
print(
	splom(~new.slope,
		varnames =c("Absolute\nError\n(%)",
			"Jittered\nPercent\n(%)",
			"Distance\n(degrees)",
			"Resolution\n(degrees)"),
	sub = list("Figure 4.35",cex=.8))
)
graphics.off()
}

book4.36 <- 
function()
{
trellis.device(postscript, file = "book4.36.ps", color = TRUE, paper = "legal")
	print(xyplot(error ~ distance | percent,
		data = slope, 
		aspect="xy",
		layout = c(11,1),
		panel = function(x, y) {
	                panel.grid(v=2)
	                panel.xyplot(x, y)
	        },
		strip = FALSE,
                sub = list("Figure 4.36",cex=.8),
		xlab="Distance (degrees)",
		ylab="Absolute Error (%)"),
		position=c(0,0,1,.6), more=T)

	print(plot(shingle(slope$percent), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Percent (%)"),
                position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.37 <- 
function()
{
trellis.device(postscript, file = "book4.37.ps", color = TRUE, paper = "legal")
	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)
	        },
		strip = FALSE,
                sub = list("Figure 4.37",cex=.8),
		xlab="Resolution (degrees)",
		ylab="Absolute Error (%)"),
                position=c(0,0,1,.6), more=T)

        print(plot(shingle(slope$percent),
		scales=list(x=list(cex=.7),y=list(cex=.7)),
                xlab = "Percent (%)"),
                position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.38 <- 
function()
{

Resolution <- equal.count(resolution,8,1/4)
trellis.device(postscript, file = "book4.38.ps", color = TRUE, paper = "legal")
	ans <- xyplot(error ~ percent | Resolution,
		aspect=1,
		layout=c(8,1),
		panel = function(x, y) {
	                panel.grid()
	                panel.xyplot(x, y)
	        },
		strip = FALSE,
                sub = list("Figure 4.38",cex=.8),
		xlab="Percent (%)",
		ylab="Absolute Error (%)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	print(plot(equal.count(slope$resolution,8,1/4), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Resolution (degrees)"),
		position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.39 <- 
function()
{
 wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}

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)
}
trellis.device(postscript, file = "book4.39.ps", color = TRUE, paper = "legal")
print(
	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 (%)")
)
graphics.off()

}

book4.4 <- 
function()
{

Tensile.strength <- equal.count(tensile.strength,6,3/4)
trellis.device(postscript, file = "book4.4.ps", color = TRUE, paper = "legal")
	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(6, 1),
		strip = FALSE,
                sub = list("Figure 4.4",cex=.8),
		xlab = "Hardness (degrees Shore)",
		ylab = "Abrasion Loss (g/hp-hour)")
	
	print(ans, position=c(0,0,1,.6), more=T)

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

book4.40 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}

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)
}
trellis.device(postscript, file = "book4.40.ps", color = TRUE, paper = "legal")
print(
	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 (%)")
)
graphics.off()

}

book4.41 <- 
function()
{
wt.bisquare = function(u, c = 4.685) {
     U <- abs(u/c)
     w <- ((1. + U) * (1. - U))^2.
     w[U > 1.] <- 0.
     w
}

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)
}
trellis.device(postscript, file = "book4.41.ps", color = TRUE, paper = "legal")
print(
	rfs(slo.lm, 
		sub = list("Figure 4.41",cex=.8),
		aspect = 1,
		ylab="Residual Absolute Error (%)")
)
graphics.off()

}

book4.42 <- 
function()
{

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
trellis.device(postscript, file = "book4.42.ps", color = TRUE, paper = "legal")
	ans <- xyplot(a ~ d | Percent,
		aspect="xy",
		layout=c(11,1),
		panel = function(x, y) {
	                panel.grid(v=2)
	                panel.xyplot(x, y, type = "l")
	        },
		strip = FALSE,
                sub = list("Figure 4.42",cex=.8),
		xlab="Distance (degrees)",
		ylab="Absolute Error (%)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	print(plot(shingle(seq(min(slope$percent), 
		max(slope$percent), length = 11)), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Percent (%)"),
		position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.43 <- 
function()
{

set.seed(19)
trellis.device(postscript, file = "book4.43.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()

}

book4.44 <- 
function()
{

set.seed(19)
Velocity <- equal.count(velocity,15,1/4)
trellis.device(postscript, file = "book4.44.ps", color = TRUE, paper = "legal")
	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)
		},
		strip = FALSE, 
		aspect=diff(range(north.south))/diff(range(east.west)),
		sub = list("Figure 4.44",cex=.8),
		xlab = "Jittered East-West Coordinate (arcsec)",
		ylab = "Jittered South-North Coordinate (arcsec)")
	
	print(ans, position=c(0,0,1,.6), more=T)

	print(plot(equal.count(galaxy$velocity,15,1/4), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Velocity (km/sec)"),
		position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.45 <- 
function()
{
trellis.device(postscript, file = "book4.45.ps", color = TRUE, paper = "legal")
	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(7, 1),
		strip = FALSE,
		xlab = "Radial Position (arcsec)",
		ylab = "Velocity (km/sec)")
        print(ans, position=c(0,0,1,.6), 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,.6,1,1))
        invisible()
graphics.off()
}


book4.46 <- 
function()
{

gal.m <- loess(velocity ~ east.west * north.south, 
	span = 0.25, degree = 2, normalize = F, family = "symmetric")
slitfit <- fitted.values(gal.m)
trellis.device(postscript, file = "book4.46.ps", color = TRUE, paper = "legal")
	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)
			panel.lines(x[k], slitfit[subscripts][k],
				lwd = add.line$lwd, lty = add.line$lty, col = add.line$col)
		}),
		layout = c(7, 1),
		strip = FALSE, 
                sub = list("Figure 4.46",cex=.8),
		xlab = "Radial Position (arcsec)",
		ylab = "Velocity (km/sec)")
	
        print(ans, position=c(0,0,1,.6), more=T)

	print(plot(shingle(galaxy$angle), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
		xlab = "Slit (degrees)"),
                position=c(0,.6,1,1))
        invisible()
graphics.off()
}


book4.47 <- 
function()
{
trellis.device(postscript, file = "book4.47.ps", color = TRUE, paper = "legal")
	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,
		strip = FALSE,
                scales=list(x=list(cex=.7),y=list(cex=.7)),
		layout = c(7, 1),
		xlab = "Radial Position (arcsec)",
		ylab = "Residual Velocity (km/sec)")
        print(ans, position=c(0,0,1,.6), more=T)

	print(plot(shingle(galaxy$angle), 
		scales=list(x=list(cex=.7),y=list(cex=.7)),
        	xlab = "Slit (degrees)"),
                position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.48 <- 
function()
{
trellis.device(postscript, file = "book4.48.ps", color = TRUE, paper = "legal")
print(
	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, ...) {
	                panel.qqmathline(x, distribution = qnorm)
	                panel.qqmath(x, ...)
	        },
		aspect=1,
		sub = list("Figure 4.48",cex=.8),
		xlab="Unit Normal Quantile",
		ylab="Residual Velocity (km/sec)")
)
graphics.off()
}

book4.49 <- 
function()
{
trellis.device(postscript, file = "book4.49.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()
}

book4.5 <- 
function()
{
trellis.device(postscript, file = "book4.5.ps", color = TRUE, paper = "legal")
print(
	splom(~ethanol, 
		sub = list("Figure 4.5",cex=.8),
		aspect = 1,
		varnames=c("NOx\n(micrograms/J)",
			"Compression\nRatio",
			"Equivalence\nRatio"))
)
graphics.off()
}

book4.50 <- 
function()
{

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)
trellis.device(postscript, file = "book4.50.ps", color = TRUE, paper = "legal")
print(
	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(...,  at = at,
				lty = ref$lty, lwd = ref$lwd, col = ref$col)
			panel.contourplot(..., 
				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)")
)
graphics.off()

}

book4.55 <- 
function()
{

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]))
trellis.device(postscript, file = "book4.55.ps", color = TRUE, paper = "legal")
	ans <- xyplot(galaxy.grid$n ~ galaxy.grid$e | Velocity,
		layout = c(16, 3),
		panel = function(x, y) {
			panel.grid(v=2, h=2)
			panel.xyplot(x, y, pch = ".")
		},
		strip = FALSE,
		aspect=diff(range(galaxy.grid$n))/diff(range(galaxy.grid$e)),
		xlab = "East-West Coordinate (arcsec)",
                sub = list("Figure 4.55",cex=.8),
		ylab = "South-North Coordinate (arcsec)")
	
        print(ans, position=c(0,0,1,.6), more=T)
	
	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)),
		ylab="",
		xlab = "Velocity (km/sec)")
	
        print(ans, position=c(0,.6,1,1))
        invisible()
graphics.off()
}

book4.56 <- 
function()
{

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))
trellis.device(postscript, file = "book4.56.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.57 <- 
function()
{

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))
trellis.device(postscript, file = "book4.57.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.58 <- 
function()
{

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))
trellis.device(postscript, file = "book4.58.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.59 <- 
function()
{

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))
trellis.device(postscript, file = "book4.59.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.6 <- 
function()
{

Equivalence.Ratio <- equal.count(E, number = 9, overlap = 0.25)
trellis.device(postscript, file = "book4.6.ps", color = TRUE, paper = "legal")
	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(9, 1),
		sub = list("Figure 4.6",cex=.8), 
		strip = FALSE,
		xlab = "Compression Ratio",
		ylab = "NOx (micrograms/J)")
	
	print(ans, position=c(0,0,1,.6), 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,.6,1,1))
        invisible()
graphics.off()
}

book4.60 <- 
function()
{

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)
trellis.device(postscript, file = "book4.60.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.61 <- 
function()
{

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)
trellis.device(postscript, file = "book4.61.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.62 <- 
function()
{

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)
trellis.device(postscript, file = "book4.62.ps", color = TRUE, paper = "legal")
print(
	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")
)
graphics.off()

}

book4.63 <- 
function()
{
ar <- diff(range(galaxy$north.south))/diff(range(galaxy$east.west))
trellis.device(postscript, file = "book4.63.ps", color = TRUE, paper = "legal")
	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=.5),
		ylab=list("South-North",cex=.5),
		zlab=list("Velocity",cex=.5))
	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()
graphics.off()
}

book4.64 <- 
function()
{
trellis.device(postscript, file = "book4.64.ps", color = TRUE, paper = "legal")
print(
	xyplot(northing ~ easting, 
		data = soil,
		panel = function(x, y) panel.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)")
)
graphics.off()
}

book4.65 <- 
function()
{
trellis.device(postscript, file = "book4.65.ps", color = TRUE, paper = "legal")
print(
	xyplot(northing ~ resistivity | track,
		data = soil,
		subset = is.ns,
		layout=c(4,2),
		panel=function(x, y) {
	                panel.grid(v=2)
	                panel.points(x, y, pch=".", cex=1.25)
	        },
		sub = list("Figure 4.65",cex=.8),
		xlab="Resistivity (ohm-cm)",
		ylab="Northing (km)")
)
graphics.off()
}

book4.66 <- 
function()
{
trellis.device(postscript, file = "book4.66.ps", color = TRUE, paper = "legal")
print(
	xyplot(resistivity[!is.ns] ~ easting[!is.ns] | track[!is.ns],
		data = soil,
		layout=c(5,8),
		panel=function(x, y) {
	                panel.grid(h=2)
	                panel.points(x, y, pch=".", cex=1.25)
	        },
		sub = list("Figure 4.66",cex=.8),
		xlab="Easting (km)",
		ylab="Resistivity (ohm-cm)")
)
graphics.off()
}

book4.67 <- 
function()
{

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]))
trellis.device(postscript, file = "book4.67.ps", color = TRUE, paper = "legal")
	ans <- xyplot(soi.grid$n ~ soi.grid$e | Resistivity,
		layout = c(14, 2),
		strip = FALSE, 
		panel = function(x, y) {
	                panel.grid(v=2)
	                panel.points(x, y, pch = ".")
	        },
		aspect=diff(range(soi.grid$n))/diff(range(soi.grid$e)),
		xlab = "Easting (km)",
                sub = list("Figure 4.67", cex = 0.8),
		ylab = "Northing (km)")
	print(ans, position=c(0,0,1,.6), more=T)
	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,
		scales=list(x=list(cex=.7),
			y=list(cex=.7,at=seq(5,25,5))),
		xlab = "Resistivity (ohm-cm)")
	print(ans, position=c(0,.6,1,1))
	
	invisible()
graphics.off()
}

book4.68 <- 
function()
{

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)
trellis.device(postscript, file = "book4.68.ps", color = TRUE, paper = "legal")
print(
	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)")
)
graphics.off()

}

book4.7 <- 
function()
{

Compression.Ratio <- C
trellis.device(postscript, file = "book4.7.ps", color = TRUE, paper = "legal")
	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(5,1),
		strip = FALSE,
		xlab = "Equivalence Ratio",
		ylab = "NOx (micrograms/J)")
	
        print(ans, position=c(0,0,1,.6), 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,.6,1,1))
        invisible()
graphics.off()
}

book4.8 <- 
function()
{
trellis.device(postscript, file = "book4.8.ps", color = TRUE, paper = "legal")
print(
	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)"))
)
graphics.off()
}

book4.9 <- 
function()
{
trellis.device(postscript, file = "book4.9.ps", color = TRUE, paper = "legal")
print(          
	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)"))
)
graphics.off()
}

book5.1 <- 
function()
{
env <- environmental
env$ozone <- env$ozone ^ (1/3)
trellis.device(postscript, file = "book5.1.ps", color = TRUE, paper = "legal")
print(
	splom(~env,
		sub = list("Figure 5.1",cex=.8),
		varnames = c("Cube Root\nOzone\n(cube root ppb)",
			"Solar\nRadiation\n(langleys)",
			"Temperature\n(Degrees Fahrenheit)",
			"Wind Speed\n(mph)"))
)
graphics.off()
}

book5.10 <- 
function()
{

data(environmental)
env.m <- loess((environmental$ozone^(1/3))~radiation*temperature*wind,
	parametric=c("radiation","wind"), span = 1, degree = 2)
r.marginal <- seq(min(radiation),max(radiation),length=5)
t.marginal <- seq(61,92,length=5)
w.marginal <- seq(4,16,length=50)
twr.marginal <- list(temperature=t.marginal,wind=w.marginal,radiation=r.marginal)
twr.grid <- expand.grid(twr.marginal)
env.fit <- predict(env.m,twr.grid)
ind <- (twr.grid$wind>=4)&
	(twr.grid$wind<=16)&
	(twr.grid$temperature<(-2*twr.grid$wind+112))&
	(twr.grid$temperature>(-2*twr.grid$wind+87))&
	(twr.grid$temperature>=61)&
	(twr.grid$temperature<=92)
env.fit[!ind] <- NA
Radiation <- twr.grid$radiation
Temp <- twr.grid$temperature
trellis.device(postscript, file = "book5.10.ps", color = TRUE, paper = "legal")
print(
	xyplot(env.fit ~ twr.grid$wind | Radiation * Temp,
	        panel=function(x, y) {
			panel.grid(h = 2, v = 2)
			panel.xyplot(x, y, type = "l")
	        },
		aspect = 2,
		layout = c(13,2), 
		sub = list("Figure 5.10",cex=.8),
		xlab = "Wind Speed (mph)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.11 <- 
function()
{
	
	data(environmental)
	env.m <- loess((environmental$ozone^(1/3))~wind*radiation*temperature,
		parametric=c("radiation","wind"),
		span = 1,
		degree = 2)
	r.marginal <- seq(min(radiation),max(radiation),length=5)
	t.marginal <- seq(61,92,length=50)
	w.marginal <- seq(4,16,length=5)
	#r.marginal <- seq(50,275,length=50)
	twr.marginal <-
list(temperature=t.marginal,wind=w.marginal,radiation=r.marginal)
	twr.grid <- expand.grid(twr.marginal)
	env.fit <- predict(env.m,twr.grid)
	ind <- (twr.grid$wind>=4)&
		(twr.grid$wind<=16)&
		(twr.grid$temperature<(-2*twr.grid$wind+112))&
		(twr.grid$temperature>(-2*twr.grid$wind+87))&
		(twr.grid$temperature>=61)&
		(twr.grid$temperature<=92)
	env.fit[!ind] <- NA
	Radiation <- twr.grid$radiation
	Wind <- twr.grid$wind
trellis.device(postscript, file = "book5.11.ps", color = TRUE, paper = "legal")
print(
	xyplot(env.fit ~ twr.grid$temperature | Radiation * Wind,
	        panel = function(x, y) {
	          panel.grid(h = 2, v = 2)
	          panel.xyplot(x, y, type = "l")
	        },
		aspect = 2,
                layout = c(13,2),  
		sub = list("Figure 5.11",cex=.8),
		xlab = "Temperature (Degrees Fahrenheit)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.12 <- 
function()
{
	
	data(environmental)
	env.m <- loess(environmental$ozone ~ radiation * temperature * wind,
		parametric = c("radiation", "wind"), span = 1, degree = 2)
trellis.device(postscript, file = "book5.12.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(env.m))) ~ fitted.values(env.m),
	        panel = function(x, y){
			panel.xyplot(x, y)
			panel.loess(x, y, span = 1, family = "g")
		},
		aspect = 1,
		las = 0, # tick labels parallel to axis
		sub = list("Figure 5.12",cex=.8),
		xlab = "Fitted Ozone (ppb)",
		ylab = "Square Root Absolute Residual Ozone (square root ppb)")
)
graphics.off()

}

book5.13 <- 
function()
{
	
	data(environmental)
	env.m <- loess(log(environmental$ozone,2) ~ radiation * temperature * wind,
		parametric = c("radiation", "wind"), span = 1, degree = 2)
trellis.device(postscript, file = "book5.13.ps", color = TRUE, paper = "legal")
print(
	xyplot(sqrt(abs(residuals(env.m))) ~ fitted.values(env.m),
		panel = function(x, y){
			panel.xyplot(x, y)
			panel.loess(x, y, span = 1, family = "g")
		},
		aspect = 1,
		las = 0, # tick labels parallel to axis
		sub = list("Figure 5.13",cex=.8),
		xlab = "Fitted Log Ozone (log 2 ppb)",
		ylab = "Square Root Absolute Residual\nLog Ozone (square root absolute log 2 ppb)")
)
graphics.off()

}

book5.14 <- 
function()
{
	
	data(environmental)
	env.m <- loess((environmental$ozone^(1/3)) ~ radiation * temperature * wind,
		parametric = c("radiation", "wind"), span = 1, degree = 2, surface = "d")
trellis.device(postscript, file = "book5.14.ps", color = TRUE, paper = "legal")
print(
	ans <- xyplot(sqrt(abs(residuals(env.m))) ~ fitted.values(env.m),
		panel = function(x, y){
			panel.xyplot(x, y)
			panel.loess(x, y, span = 1, family = "g")
		},
		aspect = 1,
		las = 0, # tick labels parallel to axis
		ylim = c(0, 1.2),
		sub = list("Figure 5.14",cex=.8),
		xlab = "Fitted Cube Root Ozone (cube root ppb)",
		ylab = "Square Root Absolute Residual\nCube Root Ozone (sixth root ppb)")
)
graphics.off()

}

book5.15 <- 
function()
{
data(environmental)
trellis.device(postscript, file = "book5.15.ps", color = TRUE, paper = "legal")
print(
	qqmath(~loess((environmental$ozone^(1/3))~radiation*temperature*wind,
			data=environmental, parametric=c("radiation","wind"), 
			span = 1, degree = 2)$residuals,
		prepanel = prepanel.qqmathline, 
	        panel = function(x, ...){
			panel.qqmath(x,...)
			panel.qqmathline(x, distribution = qnorm)
		},
		aspect=1,
		sub = list("Figure 5.15",cex=.8),
		xlab = "Unit Normal Quantile",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()
}

book5.16 <- 
function()
{
data(environmental)
trellis.device(postscript, file = "book5.16.ps", color = TRUE, paper = "legal")
print(
	rfs(loess((environmental$ozone^(1/3))~radiation*temperature*wind,
			data = environmental, parametric=c("radiation","wind"),
			span=1, degree=2),
	        panel = function(x, ...){
		        panel.grid()
		        panel.qqmath(x,...)
		},
		sub = list("Figure 5.16",cex=.8),
		aspect=1.5,
		ylab="Cube Root Ozone (cube root ppb)")
)
graphics.off()
}

book5.17 <- 
function()
{
weight <- unlist(log(hamster,2))
organ <- factor(unlist(col(as.matrix(hamster))),labels=names(hamster))
organ <- reorder(organ,weight,median)
trellis.device(postscript, file = "book5.17.ps", color = TRUE, paper = "legal")
print(
	bwplot(organ ~ weight,
		aspect = 0.3,
		sub = list("Figure 5.17",cex=.8),
		xlab="Log Organ Weight (log 2 grams)")
)
graphics.off()
}

book5.18 <- 
function()
{
trellis.device(postscript, file = "book5.18.ps", color = TRUE, paper = "legal")
print(
	splom(~log(hamster,2),
	sub = list("Figure 5.18",cex=.8), varnames = names(hamster))
)
graphics.off()
}

book5.19 <- 
function()
{
trellis.device(postscript, file = "book5.19.ps", color = TRUE, paper = "legal")
print(          
	splom(~log(hamster, 2), 
		varnames = names(hamster),
		sub = list("Figure 5.19",cex=.8),
		panel = function(x, y){
			plot.symbol <- trellis.par.get("plot.symbol")
			panel.xyplot(x[-38], y[-38])
			panel.points(x[38], y[38], pch = "+", cex = 1.5, font = plot.symbol$font, col = plot.symbol$col)
		})
)
graphics.off()
}

book5.2 <- 
function()
{

data(environmental)
Temperature <- equal.count(temperature, 4, 1/2)
Wind <- equal.count(wind, 4, 1/2)
trellis.device(postscript, file = "book5.2.ps", color = TRUE, paper = "legal")
print(
	xyplot((environmental$ozone^(1/3)) ~ radiation | Temperature * Wind,
		prepanel = function(x, y) 
			prepanel.loess(x, y, span = 1),
		panel = function(x, y){
			panel.grid(h = 2, v = 2)
			panel.xyplot(x, y, cex = 1)
			panel.loess(x, y, span = 1)
		},
		layout = c(8,2),
		aspect = 2,
		sub = list("Figure 5.2",cex=.8),
		xlab = "Solar Radiation (langleys)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.20 <- 
function()
{
new.iris <- as.matrix(iris[,1:4])
set.seed(19)
for(i in 1:4){
	new.iris[,i] <- jitter(new.iris[,i])
}
variety <- iris[,5]
n <- length(levels(variety))
trellis.device(postscript, file = "book5.20.ps", color = TRUE, paper = "legal")
print(
	splom(~new.iris,
		varnames = c("Sepal Length\n(cm)", "Sepal Width\n(cm)", "Petal Length\n(cm)", "Petal Width\n(cm)"),
		panel = panel.superpose,
		groups = variety,
		sub = list("Figure 5.20",cex=.8),
		key = list(points = Rows(trellis.par.get("superpose.symbol"), 1:n),
			text = list(paste("Iris", levels(variety))), 
			columns = n))
)
graphics.off()
}

book5.21 <- 
function()
{
data(iris)
set.seed(19)
petal.length <- iris[,3]
petal.width <- iris[,4]	
variety <- iris[,5]
n <- length(levels(variety))
mea <- (log(petal.length,2)+log(petal.width,2))/2
dif <- jitter(log(petal.length,2)-log(petal.width,2), 2)
trellis.device(postscript, file = "book5.21.ps", color = TRUE, paper = "legal")
print(
	xyplot(dif ~ mea,
		panel = function(...){
			panel.superpose(...)
			panel.abline(v = c(0.4, 1.46))
		},
		groups = variety,
		aspect = 1,
		sub = list("Figure 5.21",cex=.8),
		xlab = "Size (log 2 cm)", 
		ylab = "Jittered Elongation (log 2 ratio)",
	        key = list(points = Rows(trellis.par.get("superpose.symbol"), 1:n),
			text = list(paste("Iris", levels(variety))),
			columns = n))
)
graphics.off()
}

book5.3 <- 
function()
{

data(environmental)
Temperature <- equal.count(temperature, 4, 1/2)
Wind <- equal.count(wind, 4, 1/2)
trellis.device(postscript, file = "book5.3.ps", color = TRUE, paper = "legal")
print(
	xyplot((environmental$ozone^(1/3)) ~ radiation | Temperature * Wind,
		prepanel = function(x, y) 
			prepanel.loess(x, y, span = 1),
		panel = function(x, y){
			panel.grid(h = 2, v = 2)
			panel.xyplot(x, y, cex = .5)
			panel.loess(x, y, span = 1)
		},
		layout = c(8,2),
		aspect = 2,
		sub = list("Figure 5.3",cex=.8),
		xlab = "Solar Radiation (langleys)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.4 <- 
function()
{

data(environmental)
Radiation <- equal.count(radiation, 4, 1/2)
Temperature <- equal.count(temperature, 4, 1/2)
trellis.device(postscript, file = "book5.4.ps", color = TRUE, paper = "legal")
print(
	xyplot((environmental$ozone^(1/3)) ~ wind | Radiation * Temperature,
		prepanel = function(x, y) 
			prepanel.loess(x, y, span = 1),
	        panel = function(x, y){
	                panel.grid(h = 2, v = 2)
	                panel.xyplot(x, y, cex = .5)
	                panel.loess(x, y, span = 1)
	        },
		aspect = 2,
		layout = c(8,2),
		sub = list("Figure 5.4",cex=.8),
		xlab = "Wind Speed (mph)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.5 <- 
function()
{

Radiation <- equal.count(radiation, 4, 1/2)
Wind <- equal.count(wind, 4, 1/2)
trellis.device(postscript, file = "book5.5.ps", color = TRUE, paper = "legal")
print(
	xyplot((environmental$ozone^(1/3)) ~ temperature | Radiation * Wind,
		prepanel = function(x, y)
			prepanel.loess(x, y, span = 1),
	        panel = function(x, y) {
	                panel.grid(h = 2, v = 2)
	                panel.xyplot(x, y, cex = .5)
	                panel.loess(x, y, span = 1)
	        },
		layout = c(8,2),
		aspect = 2,
		sub = list("Figure 5.5",cex=.8),
		xlab = "Temperature (degrees F)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.6 <- 
function()
{

env.tm <- loess((environmental$ozone^(1/3))~wind,span=2/3,degree=1,surface="d")
ws <- seq(min(wind),35,length=50)
env.wfit <- predict(env.tm, ws)
xlim <- range(ws)
ylim <- range(env.wfit, environmental$ozone^(1/3))
trellis.device(postscript, file = "book5.6.ps", color = TRUE, paper = "legal")
print(
	xyplot((environmental$ozone^(1/3)) ~ wind, 
		prepanel = substitute(function(x, y)
			list(dx = diff(ws), dy = diff(env.wfit))),
		panel = substitute(function(x,y){
			add.line <- trellis.par.get("add.line")
			panel.xyplot(x,y)
			panel.lines(ws, env.wfit, lwd = add.line$lwd, 
				lty = add.line$lty, col = add.line$col)
		}),
		xlim = xlim, 
		ylim = ylim,
		aspect = 1,
		sub = list("Figure 5.6",cex=.8),
		xlab="Wind Speed (mph)",
		ylab="Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book5.7 <- 
function()
{

trellis.device(postscript, file = "book5.7.ps", color = TRUE, paper = "legal")
	left <- xyplot(environmental$r ~ environmental$t,
		aspect=1,
		panel=function(x,y){
			panel.xyplot(x,y)
			panel.abline(h=range(y))
		},
		sub = list("Figure 5.7",cex=.8),
		xlab = "Temperature (degrees F)",
		ylab = "Radiation (langleys)")
	right <- xyplot(environmental$t ~ environmental$w,
		aspect=1,
		panel=function(x,y){
			panel.xyplot(x,y)
			panel.abline(v=c(4,16),h=c(61,92))
			panel.abline(112,-2)
			panel.abline(87,-2)
		},
		sub = list("Figure 5.7",cex=.8),
		xlab = "Wind Speed (mph)",
		ylab = "Temperature (degrees F)")
	print(left, split = c(1,1,2,1), more = T)
	print(right, split = c(2,1,2,1))
	invisible()
graphics.off()
}

book5.8 <- 
function()
{

ind <- (wind > 4)&
	(wind < 16)&
	(temperature < (-2 * wind + 112))&
	(temperature > (-2 * wind + 87))&
	(temperature > 61)&
	(temperature < 92)
trellis.device(postscript, file = "book5.8.ps", color = TRUE, paper = "legal")
print(
	splom(~environmental[ind, 2:4],
		 sub = list("Figure 5.8",cex=.8),
		varnames = c("Solar\nRadiation\n(langleys)",
			"Temperature\n(Degrees Fahrenheit)",
			"Wind Speed\n(mph)"))
)
graphics.off()
}

book5.9 <- 
function()
{

env.m <- loess((environmental$ozone^(1/3))~radiation*temperature*wind,
	parametric=c("radiation","wind"),span = 1, degree = 2)
r.marginal <- seq(min(radiation),max(radiation),length=50)
t.marginal <- seq(61,92,length=5)
w.marginal <- seq(4,16,length=5)
twr.marginal <- list(temperature=t.marginal,wind=w.marginal,radiation=r.marginal)
twr.grid <- expand.grid(twr.marginal)
env.fit <- predict(env.m,twr.grid)
ind <- (twr.grid$wind>=4)&
	(twr.grid$wind<=16)&
	(twr.grid$temperature<(-2*twr.grid$wind+112))&
	(twr.grid$temperature>(-2*twr.grid$wind+87))&
	(twr.grid$temperature>=61)&
	(twr.grid$temperature<=92)
env.fit[!ind] <- NA
Temp <- twr.grid$temperature
Wind <- twr.grid$wind
trellis.device(postscript, file = "book5.9.ps", color = TRUE, paper = "legal")
print(
	xyplot(env.fit ~ twr.grid$radiation | Temp * Wind,
		panel = function(x, y)
			if(!all(is.na(y))) {
				panel.grid(h = 2, v = 2)
				panel.xyplot(x, y, type = "l")
			},
		aspect = 2,
		layout = c(13,2),
		sub = list("Figure 5.9",cex=.8),
		xlab = "Solar Radiation (langleys)",
		ylab = "Cube Root Ozone (cube root ppb)")
)
graphics.off()

}

book6.1 <- 
function()
{

logcount <- log10(count)
new.country <- reorder(country, logcount, median)
new.livestock <- reorder(livestock.type, logcount, median)
trellis.device(postscript, file = "book6.1.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.country ~ logcount | new.livestock,
		sub = list("Figure 6.1",cex=.8),
		xlab = "Log 10 Number of Livestock",
		layout=c(3,2))
)
graphics.off()

}

book6.10 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

logcount <- log(count,2)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
liv.res <- residuals(liv.lm)
liv.res[(livestock.type=="Pigs")&(country=="Turkey")] <- NA
trellis.device(postscript, file = "book6.10.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.livestock ~ liv.res | new.country,
		panel = function(x, y) {
			panel.dotplot(x, y)
			panel.abline(v=0)
		},
		aspect=1,
		layout=c(9,3),
		sub = list("Figure 6.10",cex=.8),
		xlab = "Residual Log 2 Number of Livestock")
)
graphics.off()

}

book6.11 <- 
function()
{
### requires map library -- with the world database (available from statlib)
### To install library maps, R version later than R-2.9.1 is required.
### To install library maps manually, download maps at:   http://cran.r-project.org/web/packages/maps/index.html
### run: > R CMD INSTALL maps
       wt.bisquare<-function(u, c = 4.685) {
	       U <- abs(u/c)
	       w <- ((1. + U) * (1. - U))^2.
               w[U > 1.] <- 0.
               w
	}
	
	library(maps)
	xlim <- c(-10,32)
	ylim <- c(35,70)
	 map.ar <- diff(ylim)/(diff(xlim) * cos(abs(mean(ylim))*pi/180))
	"map.xy"<-
	list(x = c(20.0565758, 9.55326843, 21.5958538, 25.9420471, 14.9860134, 
		8.19508266, -8.19369793, 14.7143764, 9.19108582, 19.241663, 
		4.75434399, -7.74096918 , 25.0365906, 17.159111, 5.38816404, 
		12.4507322, 19.241663, 12.5412788, -3.39477348
		, 19.0605717, 24.9460449, -1.67440414, 2.67179179, 8.28562832, 
		29.835516, 28.2962379), 
		y = c(40.9712944, 61.4110756, 39.8017159, 62.6363487, 
		62.7477379, 46.9862709, 40.0801888, 47.7102966, 56.1758194, 
		47.0419655, 50.7734795, 53.2797165, 42.7535095, 49.4368172, 
		51.8873634, 52.4443054, 44.1458664, 42.8648987, 40.1358795, 
		52.1101379, 45.9280815, 52.6113892, 46.8748817, 50.439312, 
		38.799221, 54.4492989))
	country.names <- levels(country)
	names(map.xy$y) <- names(map.xy$x) <-country.names
	europe <- map("world", xlim = xlim, ylim = ylim, plot = F)
	logcount <- log(count,2)
	wt <- rep(1,length(logcount))
	for(i in 1:10){
	        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
	        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
	}
	residual.sheep <- residuals(liv.lm)[livestock.type=="Sheep"]
	names(residual.sheep) <- as.character(country[livestock.type=="Sheep"])
	residual.sheep <- residual.sheep[country.names]
	residual.sheep <- equal.count(residual.sheep, number = 3, overlap = 1/4)
	trellis.device(postscript, file = "book6.11.ps", color = TRUE, paper = "legal")
	ans <- xyplot(map.xy$y ~ map.xy$x | residual.sheep,
		panel = substitute(function(x, y){
			do.call("panel.lines", c(list(x = europe), 
			trellis.par.get("reference.line")))
			do.call("panel.xyplot", c(list(x = x, y = y), 
			trellis.par.get("dot.symbol")))
		}),
		layout = c(3, 1),
	        aspect = map.ar, 
		xlim = c(-10, 32),
		ylim = c(35, 70),
		xlab = "", 
		ylab = "")
	print(ans, position=c(0,0,1,.6), more=T)
	logcount <- log(count,2)
	wt <- rep(1,length(logcount))
	for(i in 1:10){
	        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
	        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
	}
	residual.sheep <- residuals(liv.lm)[livestock.type=="Sheep"]
	ans <- plot(equal.count(residual.sheep, number = 3, overlap = 1/4),
		xlab = "Residual Log 2 Sheep Count\n\nFigure 6.11")
	print(ans, position=c(0,.6,1,1))
	
	invisible()
	graphics.off()
}

book6.12 <- 
function()
{
trellis.device(postscript, file = "book6.12.ps", color = TRUE, paper = "legal")
print(                  
	dotplot(algorithm ~ log(time,2) | machine * input, 
		data = run.time,
		aspect = "xy",
		sub = list("Figure 6.12",cex=.8),
		xlab = "Log Run Time (log 2 seconds)")
)
graphics.off()
}

book6.13 <- 
function()
{
Machine <- levels(run.time$machine)
trellis.device(postscript, file = "book6.13.ps", color = TRUE, paper = "legal")
print(
	dotplot(algorithm ~ log(time, 2) | input, 
		data = run.time,
		panel = function(x, y, subscripts, ...) {
			do.call("panel.abline", c(list(h = y), trellis.par.get("dot.line")))
			panel.superpose(x, y, subscripts, ...)
		},
		groups = machine,
		aspect = "xy",
		layout = c(2, 3),
		sub = list("Figure 6.13",cex=.8),
		xlab = "Log Run Time (log 2 seconds)",
		key = list(points = Rows(trellis.par.get("superpose.symbol"), 1:length(Machine)),
			text = list(Machine),
			columns = length(Machine)))
)
graphics.off()
}

book6.14 <- 
function()
{
Machine <- levels(run.time$machine)
trellis.device(postscript, file = "book6.14.ps", color = TRUE, paper = "legal")
print(
	dotplot(input ~ log(time, 2) | algorithm,
		data = run.time,
		panel = function(x, y, subscripts, ...) {
			do.call("panel.abline", c(list(h = y), trellis.par.get("dot.line")))
			panel.superpose(x, y, subscripts, ...)
		},
		groups = machine,
		aspect = .5,	
		layout = c(1,3),
		sub = list("Figure 6.14",cex=.8),
		xlab = "Log Run Time (log 2 seconds)",
		key = list(points = Rows(trellis.par.get("superpose.symbol"), 1:length(Machine)),
			text = list(Machine),
			columns = length(Machine)))
)
graphics.off()
}

book6.15 <- 
function()
{
	differences <- run.time
	differences$time <- log(differences$time,2)
	new.vax <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="vax")]
	new.mips <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="mips")]
	differences$time <- differences$time - c(rep(new.vax,3),rep(new.mips,3))
	differences <- differences[differences$algorithm!="new",]
	differences$algorithm <- ordered(differences$algorithm,
		names(sort(tapply(differences$time,differences$algorithm,mean))))
	differences$machine <- ordered(differences$machine,
		names(sort(tapply(differences$time,differences$machine,mean))))
	differences$input <- ordered(differences$input,
		names(sort(tapply(differences$time,differences$input,mean))))
	Algorithm <- levels(differences$algorithm)
	trellis.device(postscript, file = "book6.15.ps", color = TRUE, paper = "legal")
	print(
	dotplot(input ~ time | machine,
		data = differences, 
		panel = function(x, y, subscripts, ...){
			do.call("panel.abline", c(list(h = y), trellis.par.get("dot.line")))
			panel.superpose(x, y, subscripts, ...)
		},
		groups = algorithm,
		aspect = .5,
		layout = c(2, 1),
		xlim = c(0, 4),
		sub = list("Figure 6.15",cex=.8),
		xlab = "Improvement (log 2 seconds)",
		key = list(y = 1.3,
			points = Rows(trellis.par.get("superpose.symbol"), 1:length(Algorithm)),
			text = list(Algorithm),
			columns = length(Algorithm)))
	)
	graphics.off()
}



book6.16 <- 
function()
{
	differences <- run.time
	differences$time <- log(differences$time,2)
	new.vax <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="vax")]
	new.mips <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="mips")]
	differences$time <- differences$time - c(rep(new.vax,3),rep(new.mips,3))
	differences <- differences[differences$algorithm!="new",]
	differences$algorithm <- ordered(differences$algorithm,
		names(sort(tapply(differences$time,differences$algorithm,mean))))
	differences$machine <- ordered(differences$machine,
		names(sort(tapply(differences$time,differences$machine,mean))))
	differences$input <- ordered(differences$input,
		names(sort(tapply(differences$time,differences$input,mean))))
	
	runtime.lm <- lm(time~(input+algorithm)*machine)
        trellis.device(postscript, file = "book6.16.ps", color = TRUE, paper = "legal")
	print(
	    dotplot(input ~ residuals(runtime.lm) | machine * algorithm, 
		panel = function(x, y) {
			panel.dotplot(x, y)
			panel.abline(v = 0)
		},
		aspect = 2/3,
		layout = c(2, 2),
		sub = list("Figure 6.16",cex=.8),
		xlab = "Residual Improvement (log 2 seconds)")
	)
	graphics.off()
	
}

book6.17 <- 
function()
{
	differences <- run.time
	differences$time <- log(differences$time,2)
	new.vax <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="vax")]
	new.mips <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="mips")]
	differences$time <- differences$time - c(rep(new.vax,3),rep(new.mips,3))
	differences <- differences[differences$algorithm!="new",]
	differences$algorithm <- ordered(differences$algorithm,
		names(sort(tapply(differences$time,differences$algorithm,mean))))
	differences$machine <- ordered(differences$machine,
		names(sort(tapply(differences$time,differences$machine,mean))))
	differences$input <- ordered(differences$input,
		names(sort(tapply(differences$time,differences$input,mean))))
        trellis.device(postscript, file = "book6.17.ps", color = TRUE, paper = "legal")
        print(
	  rfs(lm(time~(input+algorithm)*machine, data = differences),
		sub = list("Figure 6.17",cex=.8),
		aspect = 2, 
		ylab = "Log Run Time (log 2 seconds)")
	)
	graphics.off()
}

book6.18 <- 
function()
{
	differences <- run.time
	differences$time <- log(differences$time,2)
	new.vax <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="vax")]
	new.mips <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="mips")]
	differences$time <- differences$time - c(rep(new.vax,3),rep(new.mips,3))
	differences <- differences[differences$algorithm!="new",]
	differences$algorithm <- ordered(differences$algorithm,
		names(sort(tapply(differences$time,differences$algorithm,mean))))
	differences$machine <- ordered(differences$machine,
		names(sort(tapply(differences$time,differences$machine,mean))))
	differences$input <- ordered(differences$input,
		names(sort(tapply(differences$time,differences$input,mean))))
	
	runtime.lm <- lm(time~(input+algorithm)*machine)
	which <- algorithm == "7th"
        trellis.device(postscript, file = "book6.18.ps", color = TRUE, paper = "legal")
        print(
	ans <- dotplot(input[which] ~ fitted(runtime.lm)[which] | machine[which],
		aspect = 1/4,
		layout = c(1, 2),
		sub = list("Figure 6.18",cex=.8),
		xlab = "Fitted Improvement (log 2 seconds)")
	)
	graphics.off()
	
}

book6.19 <- 
function()
{
	differences <- run.time
	differences$time <- log(differences$time,2)
	new.vax <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="vax")]
	new.mips <- differences$time[(differences$algorithm=="new")&
		(differences$machine=="mips")]
	differences$time <- differences$time - c(rep(new.vax,3),rep(new.mips,3))
	differences <- differences[differences$algorithm!="new",]
	differences$algorithm <- ordered(differences$algorithm,
		names(sort(tapply(differences$time,differences$algorithm,mean))))
	differences$machine <- ordered(differences$machine,
		names(sort(tapply(differences$time,differences$machine,mean))))
	differences$input <- ordered(differences$input,
		names(sort(tapply(differences$time,differences$input,mean))))
	
	runtime.lm.fitted <- 2^fitted(lm(time~(input+algorithm)*machine))
	which <- algorithm == "7th"
        trellis.device(postscript, file = "book6.19.ps", color = TRUE, paper = "legal")
        print(
	   dotplot(input ~ runtime.lm.fitted | machine,
		subset = algorithm == "7th",
		aspect = 1/4,
		layout = c(1, 2),
		sub = list("Figure 6.19",cex=.8),
		xlab = "Fitted Improvement Factor")
	)
	graphics.off()
	
}

book6.2 <- 
function()
{

logcount <- log10(count)
new.country <- reorder(country, logcount, median)
new.livestock <- reorder(livestock.type, logcount, median)
trellis.device(postscript, file = "book6.2.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.livestock ~ logcount | new.country,
		sub = list("Figure 6.2",cex=.8),
		xlab = "Log 10 Number of Livestock",
		layout=c(4,7),
		aspect=1/2)
)
graphics.off()

}

book6.20 <- 
function()
{
trellis.device(postscript, file = "book6.20.ps", color = TRUE, paper = "legal")
print(
	dotplot(variety ~ yield | year * site,
		data = barley,
		aspect=.4,
		sub = list("Figure 6.20",cex=.8),
		scales=list(y=list(cex=.4)),
		xlab = "Barley Yield (bushels/acre)")
)
graphics.off()		
}

book6.21 <- 
function()
{

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
trellis.device(postscript, file = "book6.21.ps", color = TRUE, paper = "legal")
print(
	dotplot(variety ~ new.yield | year * site, data=barley,
		sub = list("Figure 6.21",cex=.8),
		 scales=list(y=list(cex=.4)),
		xlab = "Barley Yield (bushels/acre)",
		aspect=.4)
)
graphics.off()

}

book6.22 <- 
function()
{

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
which <- year=="1931"
new.yield.diff <- new.yield[which]-new.yield[!which]
trellis.device(postscript, file = "book6.22.ps", color = TRUE, paper = "legal")
print(
	dotplot(variety[which] ~ new.yield.diff | site[which], 
		sub = list("Figure 6.22",cex=.8),
		 scales=list(y=list(cex=.4)),
		 xlab = "Differences of Barley Yield (bushels/acre)",
		layout = c(1,6),
		aspect=.4)
)
graphics.off()

}

book6.23 <- 
function()
{

wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
wt <- rep(1,length(yield))
for(i in 1:10){
        barley.lm <- lm(new.yield~variety+year*site,weights=wt)
        wt <- wt.bisquare(barley.lm$res/median(abs(barley.lm$res)),c=6)
}
trellis.device(postscript, file = "book6.23.ps", color = TRUE, paper = "legal")
print(
	rfs(barley.lm, 
		sub = list("Figure 6.23",cex=.8),
		aspect=2,
		ylab = "Yield (bushels/acre)")
)
graphics.off()

}

book6.24 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
U <- abs(u/c)
w <- ((1. + U) * (1. - U))^2.
w[U > 1.] <- 0.
w
}

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31	
wt <- rep(1,length(yield))
for(i in 1:10){
        barley.lm <- lm(new.yield~variety+year*site,weights=wt)
        wt <- wt.bisquare(barley.lm$res/median(abs(barley.lm$res)),c=6)
}
barley.effects <- dummy.coef(barley.lm)
ys.effects <- c(barley.effects$"year:site" + outer(
	barley.effects$year,barley.effects$site,"+"))
ys.year <- ordered(rep(levels(year),6),levels(year))
ys.site <- ordered(rep(levels(site),rep(2,6)),levels(site))
n <- length(levels(ys.year))
trellis.device(postscript, file = "book6.24.ps", color = TRUE, paper = "legal")
print(
	dotplot(ys.site ~ ys.effects,
		panel = function(x, y, subscripts, ...) {
			do.call("panel.abline", c(list(h = y), trellis.par.get("dot.line")))
			panel.superpose(x, y, subscripts, ...)
		},
		groups = ys.year,
		aspect = 2/3,
		sub = list("Figure 6.24",cex=.8),
		xlab = "Site by Year Effects (bushels/acre)",
		key = list(
			points = Rows(trellis.par.get("superpose.symbol"), 1:n),
			text = list(levels(ys.year)),
			columns = n))
)
graphics.off()

}

book6.25 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
	U <- abs(u/c)
	w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
wt <- rep(1,length(yield))
for(i in 1:10){
        barley.lm <- lm(new.yield~variety+year*site,weights=wt)
        wt <- wt.bisquare(barley.lm$res/median(abs(barley.lm$res)),c=6)
}
barley.effects <- dummy.coef(barley.lm)
trellis.device(postscript, file = "book6.25.ps", color = TRUE, paper = "legal")
print(
	dotplot(sort(barley.effects$variety),
		sub = list("Figure 6.25",cex=.8),
		xlab = "Variety Effects (bushels/acre)", 
		aspect=2/3)
)
graphics.off()

}

book6.26 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
wt <- rep(1,length(yield))
for(i in 1:10){
        barley.lm <- lm(new.yield~variety+year*site,weights=wt)
        wt <- wt.bisquare(barley.lm$res/median(abs(barley.lm$res)),c=6)
}
trellis.device(postscript, file = "book6.26.ps", color = TRUE, paper = "legal")
print(
	dotplot(site ~ residuals(barley.lm) | year * variety,
		layout=c(4,5),
		aspect=2/3,
		panel = function(x, y) {
			panel.dotplot(x, y)
			panel.abline(v=0)
		},
		sub = list("Figure 6.26",cex=.8),
		xlab = "Residual Barley Yield (bushels/acre)")
)
graphics.off()

}

book6.27 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

morris31 <- yield[(site=="Morris")&(year=="1931")]
morris32 <- yield[(site=="Morris")&(year=="1932")]
new.yield <- yield
new.yield[(site=="Morris")&(year=="1931")] <- morris32
new.yield[(site=="Morris")&(year=="1932")] <- morris31
wt <- rep(1,length(yield))
for(i in 1:10){
        barley.lm <- lm(new.yield~variety+year*site,weights=wt)
        wt <- wt.bisquare(barley.lm$res/median(abs(barley.lm$res)),c=6)
}
trellis.device(postscript, file = "book6.27.ps", color = TRUE, paper = "legal")
print(
	dotplot(variety ~ residuals(barley.lm) | year * site,
		layout = c(2,6),
		aspect=2/3,
		panel = function(x, y) {
			panel.dotplot(x, y)
			panel.abline(v=0)
		},
		sub = list("Figure 6.27",cex=.8),
		xlab = "Residual Barley Yield (bushels/acre)")
)
graphics.off()

}

book6.3 <- 
function()
{

logcount <- log10(count)
new.country <- ordered(country, rev(sort(levels(country))))
new.livestock <- ordered(livestock.type,
	c("Sheep","Poultry","Pigs","Horses","Cattle"))
trellis.device(postscript, file = "book6.3.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.country ~ logcount | new.livestock,
		sub = list("Figure 6.3",cex=.8),
		xlab = "Log 10 Number of Livestock",
		layout=c(3,2))
)
graphics.off()

}

book6.4 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

logcount <- log10(count)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,sort)
all.effects <- c(liv.effects$livestock.type,NA,liv.effects$country)
trellis.device(postscript, file = "book6.4.ps", color = TRUE, paper = "legal")
print(
	ans <- dotplot(all.effects,
		aspect=1,
		sub = list("Figure 6.4",cex=.8),
		xlab = "Log Livestock Effects (log 10 count)")
)
graphics.off()

}

book6.5 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

logcount <- log10(count)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
trellis.device(postscript, file = "book6.5.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.country ~ logcount | new.livestock,
		sub = list("Figure 6.5",cex=.8),
		scales=list(y=list(cex=.4)),
		xlab = "Log 10 Number of Livestock",
		layout=c(3,2))
)
graphics.off()

}

book6.6 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

logcount <- log10(count)
wt <- rep(1,length(logcount))
for(i in 1:10){
	liv.lm <- lm(logcount~livestock.type+country,weights=wt)
	wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
trellis.device(postscript, file = "book6.6.ps", color = TRUE, paper = "legal")
print(
	ans <- dotplot(new.livestock ~ logcount | new.country,
		sub = list("Figure 6.6",cex=.8),
		xlab = "Log 10 Number of Livestock",
		layout=c(9,3),
		aspect=1)
)
graphics.off()

}

book6.7 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}       

logcount <- log10(count)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
trellis.device(postscript, file = "book6.7.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.country ~ fitted.values(liv.lm) | new.livestock,
		sub = list("Figure 6.7",cex=.8),
	        xlab = "Fitted Log 10 Number of Livestock",
	        layout=c(3,2))
)
graphics.off()

}

book6.8 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}

logcount <- log10(count)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
trellis.device(postscript, file = "book6.8.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.livestock ~ fitted.values(liv.lm) | new.country,
		sub = list("Figure 6.8",cex=.8),
		xlab = "Fitted Log 10 Number of Livestock",
		layout=c(9,3),
	        aspect=1)
)
graphics.off()

}

book6.9 <- 
function()
{
wt.bisquare<-function(u, c = 4.685) {
        U <- abs(u/c)
        w <- ((1. + U) * (1. - U))^2.
        w[U > 1.] <- 0.
        w
}       

logcount <- log(count,2)
wt <- rep(1,length(logcount))
for(i in 1:10){
        liv.lm <- lm(logcount~livestock.type+country,weights=wt)
        wt <- wt.bisquare(liv.lm$res/median(abs(liv.lm$res)),c=6)
}
liv.effects <- dummy.coef(liv.lm)
liv.effects <- sapply(liv.effects,"sort")
new.country <- ordered(country,names(liv.effects$country))
new.livestock <- ordered(livestock.type,names(liv.effects$livestock.type))
liv.res <- residuals(liv.lm)
liv.res[(livestock.type=="Pigs")&(country=="Turkey")] <- NA
trellis.device(postscript, file = "book6.9.ps", color = TRUE, paper = "legal")
print(
	dotplot(new.country ~ liv.res | new.livestock,
		panel = function(x, y) {
			panel.dotplot(x, y)
	                panel.abline(v=0)
		},
		layout=c(3,2),
		sub = list("Figure 6.9",cex=.8),
		xlab = "Residual Log 2 Number of Livestock")
)
graphics.off()

}