11 for (i in 1:maxdeg) {
12 lmo <- lm(y ~ pwrs[,1:i])
13 # extend the lm class here, with the cross-validated predictions
14 lmo$fitted.cvvalues <- lvoneout(y,pwrs[,1:i,drop=F])
15 lmout[[i]] <- lmo
16 }
17 lmout$x <- x
18 lmout$y <- y
19 return(lmout)
20 }
21
22 # print() for an object fits of class "polyreg": print
23 # cross-validated mean-squared prediction errors
24 print.polyreg <- function(fits) {
25 maxdeg <- length(fits) - 2
26 n <- length(fits$y)
27 tbl <- matrix(nrow=maxdeg,ncol=1)
28 colnames(tbl) <- "MSPE"
29 for (i in 1:maxdeg) {
30 fi <- fits[[i]]
31 errs <- fits$y - fi$fitted.cvvalues
32 spe <- crossprod(errs,errs) # sum of squared prediction errors
33 tbl[i,1] <- spe/n
34 }
35 cat("mean squared prediction errors, by degree\n")
36 print(tbl)
37 }
38
39 # forms matrix of powers of the vector x, through degree dg
40 powers <- function(x,dg) {
41 pw <- matrix(x,nrow=length(x))
42 prod <- x
43 for (i in 2:dg) {
44 prod <- prod*x
45 pw <- cbind(pw,prod)
46 }
47 return(pw)
48 }
49
50 # finds cross-validated predicted values; could be made much faster via
51 # matrix-update methods
52 lvoneout <- function(y,xmat) {
53 n <- length(y)
54 predy <- vector(length=n)
55 for (i in 1:n) {
56 # regress, leaving out ith observation
57 lmo <- lm(y[-i] ~ xmat[-i,])
58 betahat <- as.vector(lmo$coef)220 Chapter 9