The Art of R Programming

(WallPaper) #1
Note that we asked R to use blank labels for the figure as a whole and
for thex-axis. Otherwise, R would have gotten such labels fromd1, which
would have been specific to exam 1.
Also note that we needed to plot exam 1 first. The scores there were
less diverse, so the density estimate was narrower and taller. Had we plotted
exam 2, with its shorter curve, first, exam 1’s curve would have been too tall
for the plot window. Here, we first ran the two plots separately to see which
was taller, but let’s consider a more general situation.
Say we wish to write a broadly usable function that will plot several den-
sity estimates on the same graph. For this, we would need to automate the
process of determining which density estimate is tallest. To do so, we would
use the fact that the estimated density values are contained in theycompo-
nent of the return value from the call todensity(). We would then callmax()
on each density estimate and usewhich.max()to determine which density esti-
mate is the tallest.
The call toplot()both initiates the plot and draws the first curve. (With-
out specifyingtype="l", only the points would have been plotted.) The call to
lines()then adds the second curve.

12.1.5 Extended Example: More on the Polynomial Regression Example....


In Section 9.1.7, we defined a class"polyreg"that facilitates fitting polyno-
mial regression models. Our code there included an implementation of the
genericprint()function. Let’s now add one for the genericplot()function:

1 # polyfit(x,maxdeg) fits all polynomials up to degree maxdeg; y is
2 # vector for response variable, x for predictor; creates an object of
3 # class "polyreg", consisting of outputs from the various regression
4 # models, plus the original data
5 polyfit <- function(y,x,maxdeg) {
6 pwrs <- powers(x,maxdeg) # form powers of predictor variable
7 lmout <- list() # start to build class
8 class(lmout) <- "polyreg" # create a new class
9 for (i in 1:maxdeg) {
10 lmo <- lm(y ~ pwrs[,1:i])
11 # extend the lm class here, with the cross-validated predictions
12 lmo$fitted.xvvalues <- lvoneout(y,pwrs[,1:i,drop=F])
13 lmout[[i]] <- lmo
14 }
15 lmout$x <- x
16 lmout$y <- y
17 return(lmout)
18 }
19
20 # generic print() for an object fits of class "polyreg": print
21 # cross-validated mean-squared prediction errors
22 print.polyreg <- function(fits) {
23 maxdeg <- length(fits) - 2 # count lm() outputs only, not $x and $y

266 Chapter 12

Free download pdf