1 Setup

Libraries used in this exercise: sp for spatial objects, rpart for regression trees and random forests, rpart.plot for nicer plots of trees.

library(sp)
library(rpart)
library(rpart.plot)

Sample dataset, comes with sp; make a log10-transformed copy of the metal concentrations.

data(meuse)
str(meuse)

## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...

meuse$logZn <- log10(meuse$zinc)    # 锌
meuse$logCu <- log10(meuse$copper)  # 铜
meuse$logPb <- log10(meuse$lead)    # 铅
meuse$logCd <- log10(meuse$cadmium) # 镉

Prediction grid:

data(meuse.grid)
coordinates(meuse.grid) <- ~ x + y
gridded(meuse.grid) <- TRUE

2 The single best tree

Use all 155 to build and prune:

m.lzn <- rpart(logZn ~ ffreq + dist + soil + x + y,
                    data=meuse,
                    minsplit=2,
                    cp=0.005)
cp.table <- m.lzn[["cptable"]]
cp.min <- cp.table[which.min(cp.table[,"xerror"]),"CP"]
m.lzn <- prune(m.lzn, cp=cp.min)
rpart.plot(m.lzn)

grid.map <- predict(m.lzn, newdata=meuse.grid)
meuse.grid$pred <- grid.map
print(spplot(meuse.grid, zcol="pred"))

x <- m.lzn$variable.importance 
data.frame(variableImportance = 100 * x / sum(x))

##       variableImportance
## dist            45.53803
## x               16.88457
## soil            15.27468
## y               11.54119
## ffreq           10.76154

3 Repeated random sampling and tree construction

Randomly sample a specific number of the 155 points and use it to build and prune a regression tree; use the remaining points as out-of-bag evaluation. The model uses all the variable in the prediction grid: coordnates, relative distance from river, flooding frequency clqss, and soil type.

We write this as a function, so we can call it any number of times. It has one argument: the number of points to include in the calibration sample.

It returns information about the tree: the minimum complexity parameter, the out-of-bag error, and the variable importance, and the prediction over the grid. It also displays the tree and the prediction over the grid.

build.rpart <- function(n.cal) {
  meuse.cal.ix <- sample(1:dim(meuse)[1], size=n.cal, replace=FALSE)
  meuse.cal <- meuse[meuse.cal.ix,]
  meuse.val <- meuse[-meuse.cal.ix,]
  m.lzn.cal <- rpart(logZn ~ ffreq + dist + soil + x + y,
                    data=meuse.cal,
                    minsplit=2,
                    cp=0.005)
  cp.table <- m.lzn.cal[["cptable"]]
  cp.min <- cp.table[which.min(cp.table[,"xerror"]),"CP"]
  m.lzn.cal <- prune(m.lzn.cal, cp=cp.min)
  rpart.plot(m.lzn.cal)
  tmp <- m.lzn.cal$variable.importance
  v.imp <- data.frame(variableImportance = 100 * tmp / sum(tmp))
  grid.map <- predict(m.lzn.cal, newdata=meuse.grid)
  meuse.grid$pred <- grid.map
  val <- predict(m.lzn.cal, newdata=meuse.val)
  rmse <- sqrt((sum(meuse.val$logZn - val)^2)/length(meuse.val$logZn))
  print(spplot(meuse.grid, zcol="pred"))
  # return a list
  return(list(cp = cp.min, rmse = rmse, v.imp = v.imp, pred = grid.map))
}

Set up an object to collect the fitted paramter, and a field in the prediction grid to show the average prediction:

runs <- 64
fits.rf <- list(cp=vector(mode="numeric",length=runs),
             rmse=vector(mode="numeric",length=runs),
             v.imp=as.data.frame(matrix(data=0, nrow=runs, ncol=5))
)
names(fits.rf[["v.imp"]]) <- c("ffreq", "dist", "soil","x","y")
meuse.grid$logZn.avg <- 0  # will accumulate predictions and at the end divide by # of runs

Call the function and collect the fitted parameters. Save the graphics in a PDF for display.

for (run in 1:runs) {
  tmp <- build.rpart(140)
  fits.rf[["cp"]][run] <- tmp[[1]]
  fits.rf[["rmse"]][run] <- tmp[[2]]
  fits.rf[["v.imp"]][run,] <- t(tmp[[3]])  # returned a row, need a column
  # add to average prediction
  meuse.grid$logZn.avg <- meuse.grid$logZn.avg + tmp[[4]]
}