2.1 Building a regression tree

We will use the rpart function to build the regression tree. This uses a helper function rpart.control to control the behaviour of rpart.

Examine the help for rpart.control.

help(rpart.control)

Q: What is the meaning of the minsplit and cp parameters?

Build a regression tree for this, print its results, and display as a tree plot. Allow splits all the way to single observations, if the \(R^2\) increases by at least 0.5% at the split.

m.lzn.rp <- rpart(logZn ~ x + y + ffreq + dist.m + soil + lime + elev,
                  data=meuse,
                  minsplit=2,
                  cp=0.005)
print(m.lzn.rp)

## n= 155 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##   1) root 155 15.13633000 2.556160  
##     2) dist.m>=145 101  4.65211000 2.388584  
##       4) elev>=6.943 93  2.83545700 2.349952  
##         8) dist.m>=230 78  1.73584100 2.311717  
##          16) x>=179102.5 71  1.23063300 2.288808  
##            32) ffreq=2,3 40  0.36471060 2.223342  
##              64) y>=330942 21  0.14958940 2.167449 *
##              65) y< 330942 19  0.07701018 2.285117 *
##            33) ffreq=1 31  0.47328450 2.373281  
##              66) x< 179913.5 1  0.00000000 2.075547 *
##              67) x>=179913.5 30  0.38168440 2.383205  
##               134) x>=180346.5 25  0.20839560 2.355868  
##                 268) y< 333191 19  0.09154361 2.320174 *
##                 269) y>=333191 6  0.01598697 2.468900 *
##               135) x< 180346.5 5  0.06119400 2.519889 *
##          17) x< 179102.5 7  0.08998420 2.544084 *
##         9) dist.m< 230 15  0.39263770 2.548774  
##          18) x< 180563 10  0.23688250 2.494602  
##            36) x>=179405.5 7  0.10694750 2.420923 *
##            37) x< 179405.5 3  0.00326675 2.666521 *
##          19) x>=180563 5  0.06771894 2.657116 *
##       5) elev< 6.943 8  0.06436102 2.837680 *
##     3) dist.m< 145 54  2.34313500 2.869589  
##       6) elev>=8.1455 15  0.49618250 2.646331  
##        12) elev< 8.48 4  0.07502826 2.456464 *
##        13) elev>=8.48 11  0.22452210 2.715373 *
##       7) elev< 8.1455 39  0.81172600 2.955457  
##        14) dist.m>=75 11  0.08001190 2.848717 *
##        15) dist.m< 75 28  0.55715050 2.997391  
##          30) x< 179565 9  0.12124970 2.914541 *
##          31) x>=179565 19  0.34486030 3.036636  
##            62) x>=180194 12  0.12169000 2.966064 *
##            63) x< 180194 7  0.06095067 3.157617 *

rpart.plot(m.lzn.rp, digits=3, type=4, extra=1)

Q: Which variables were actually used in the tree?

Q: Which variable was used at the first split?

Q: Does that agree with your theory of the origin of the metals in these soils?

Q: What was the value of the splitting variable at that split?

Q: What was the average log10(Zn) before the split, and in the two leaves after the split?

Q: How many leaves (i.e., different predictions)? What proportion of the observations is this?

Q: What is the minimum number of observations among the leaves? The maximum?

2.2 Variable importance

Now we examine the variable importance:

print(tmp <- m.lzn.rp$variable.importance)

##   dist.m     elev     lime        x        y     soil    ffreq 
## 9.410837 6.134458 4.633959 2.097574 1.911179 1.169126 1.143282

data.frame(variableImportance = 100 * tmp / sum(tmp))

##        variableImportance
## dist.m          35.512035
## elev            23.148534
## lime            17.486364
## x                7.915251
## y                7.211884
## soil             4.411728
## ffreq            4.314205

Q: Which variable explained most of the metal concentration? What proportion? Were all predictor variables used somewhere in the tree?

2.3 Pruning the regression tree

Have we built a tree that is too complex, i.e., fitting noise not signal? In other words, fitting our dataset well but not the process?

Examine the cross-validation relative error vs. the complexity parameter.

Q: What complexity parameter minimizes the cross-validation relative error?

We can also compute this, rather than just looking at the graph.

plotcp(m.lzn.rp)

cp.table <- m.lzn.rp[["cptable"]]
cp.ix <- which.min(cp.table[,"xerror"])
print(cp.table[cp.ix,])

##          CP      nsplit   rel error      xerror        xstd 
## 0.009124468 8.000000000 0.153385306 0.292993740 0.032477480

cp.min <- cp.table[cp.ix,"CP"]

(Note that some other cross-validation errors are quite close to this minimum, so we might choose to use a higher value of the control parameter to have a more parsimonious model, that might extrapolate better.)

Prune the tree with this complexity parameter, and examine in the same way as above.

(Note we would get the same result if we re-fit with this parameter.)

(m.lzn.rpp <- prune(m.lzn.rp, cp=cp.min))

## n= 155 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 155 15.13633000 2.556160  
##    2) dist.m>=145 101  4.65211000 2.388584  
##      4) elev>=6.943 93  2.83545700 2.349952  
##        8) dist.m>=230 78  1.73584100 2.311717  
##         16) x>=179102.5 71  1.23063300 2.288808  
##           32) ffreq=2,3 40  0.36471060 2.223342 *
##           33) ffreq=1 31  0.47328450 2.373281 *
##         17) x< 179102.5 7  0.08998420 2.544084 *
##        9) dist.m< 230 15  0.39263770 2.548774 *
##      5) elev< 6.943 8  0.06436102 2.837680 *
##    3) dist.m< 145 54  2.34313500 2.869589  
##      6) elev>=8.1455 15  0.49618250 2.646331  
##       12) elev< 8.48 4  0.07502826 2.456464 *
##       13) elev>=8.48 11  0.22452210 2.715373 *
##      7) elev< 8.1455 39  0.81172600 2.955457  
##       14) dist.m>=75 11  0.08001190 2.848717 *
##       15) dist.m< 75 28  0.55715050 2.997391 *

rpart.plot(m.lzn.rpp, digits=3, type=4, extra=1)

Q: How does this tree differ from the full (un-pruned tree)? Look at the predictors actually used, number of splits, number of leaves, number of observations in each leaf.

Q: Which leaf predicts the highest concentration of log10(Zn)? Explain that by reading the predictors and their values down the tree from the root to this leaf. Does this fit your theory of the origin of these metals?

2.4 Goodness-of-fit

Now, compare the model predictions to the known values:

p.rpp <- predict(m.lzn.rpp, newdata=meuse)
length(unique(p.rpp))

## [1] 9

summary(r.rpp <- meuse$logZn - p.rpp)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.320887 -0.090478 -0.005171  0.000000  0.081155  0.293237

sqrt(sum(r.rpp^2)/length(r.rpp))

## [1] 0.1223873

plot(meuse$logZn ~ p.rpp, asp=1, pch=20, xlab="fitted", ylab="actual", xlim=c(2,3.3), ylim=c(2,3.3), main="log10(Zn), Meuse topsoils, Regression Tree")
grid()
abline(0,1)

Q: What is the RMSE?

Q: Why are there vertical lines at all the fitted values?

3.1 Building a classification tree

Flooding is influenced by the distance from the river and the elevation above the river. Build a model with these, using a complexity parameter of 2% increase in \(R^2\).

m.ff.rp <- rpart(ffreq ~ x + y + dist.m + elev,
                  data=meuse,
                  minsplit=2,
                  cp=0.02)
print(m.ff.rp)

## n= 155 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##   1) root 155 71 1 (0.54193548 0.30967742 0.14838710)  
##     2) y>=331984.5 63  6 1 (0.90476190 0.07936508 0.01587302)  
##       4) elev< 9.399 60  3 1 (0.95000000 0.03333333 0.01666667) *
##       5) elev>=9.399 3  0 2 (0.00000000 1.00000000 0.00000000) *
##     3) y< 331984.5 92 49 2 (0.29347826 0.46739130 0.23913043)  
##       6) elev< 7.552 24  2 1 (0.91666667 0.00000000 0.08333333)  
##        12) y>=330043 22  0 1 (1.00000000 0.00000000 0.00000000) *
##        13) y< 330043 2  0 3 (0.00000000 0.00000000 1.00000000) *
##       7) elev>=7.552 68 25 2 (0.07352941 0.63235294 0.29411765)  
##        14) y>=330330 54 13 2 (0.09259259 0.75925926 0.14814815)  
##          28) elev< 9.947 49  9 2 (0.10204082 0.81632653 0.08163265)  
##            56) y>=330644 35  5 2 (0.14285714 0.85714286 0.00000000)  
##             112) dist.m>=470 17  5 2 (0.29411765 0.70588235 0.00000000)  
##               224) elev< 8.526 3  0 1 (1.00000000 0.00000000 0.00000000) *
##               225) elev>=8.526 14  2 2 (0.14285714 0.85714286 0.00000000) *
##             113) dist.m< 470 18  0 2 (0.00000000 1.00000000 0.00000000) *
##            57) y< 330644 14  4 2 (0.00000000 0.71428571 0.28571429)  
##             114) elev< 8.7225 11  1 2 (0.00000000 0.90909091 0.09090909) *
##             115) elev>=8.7225 3  0 3 (0.00000000 0.00000000 1.00000000) *
##          29) elev>=9.947 5  1 3 (0.00000000 0.20000000 0.80000000) *
##        15) y< 330330 14  2 3 (0.00000000 0.14285714 0.85714286) *

rpart.plot(m.ff.rp)

Q: Were both predictors used?

Q: What was the first split?

3.2 Pruning the classification tree

Find the value of the complexity parameter that minimizes the cross-validation error, and prune the tree with that value.

printcp(m.ff.rp)

## 
## Classification tree:
## rpart(formula = ffreq ~ x + y + dist.m + elev, data = meuse, 
##     minsplit = 2, cp = 0.02)
## 
## Variables actually used in tree construction:
## [1] dist.m elev   y     
## 
## Root node error: 71/155 = 0.45806
## 
## n= 155 
## 
##         CP nsplit rel error  xerror     xstd
## 1 0.267606      0   1.00000 1.00000 0.087366
## 2 0.140845      2   0.46479 0.49296 0.073316
## 3 0.042254      3   0.32394 0.36620 0.065518
## 4 0.028169      5   0.23944 0.39437 0.067462
## 5 0.021127      6   0.21127 0.33803 0.063433
## 6 0.020000     10   0.12676 0.38028 0.066506

plotcp(m.ff.rp)

cp.table.class <- m.ff.rp[["cptable"]]
cp.ix <- which.min(cp.table.class[,"xerror"])
print(cp.table.class[cp.ix,])

##         CP     nsplit  rel error     xerror       xstd 
## 0.02112676 6.00000000 0.21126761 0.33802817 0.06343327

(cp.min <- cp.table.class[cp.ix,"CP"])

## [1] 0.02112676

m.ff.rpp <- prune(m.ff.rp, cp=cp.min)
print(m.ff.rpp)

## n= 155 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 155 71 1 (0.54193548 0.30967742 0.14838710)  
##    2) y>=331984.5 63  6 1 (0.90476190 0.07936508 0.01587302)  
##      4) elev< 9.399 60  3 1 (0.95000000 0.03333333 0.01666667) *
##      5) elev>=9.399 3  0 2 (0.00000000 1.00000000 0.00000000) *
##    3) y< 331984.5 92 49 2 (0.29347826 0.46739130 0.23913043)  
##      6) elev< 7.552 24  2 1 (0.91666667 0.00000000 0.08333333)  
##       12) y>=330043 22  0 1 (1.00000000 0.00000000 0.00000000) *
##       13) y< 330043 2  0 3 (0.00000000 0.00000000 1.00000000) *
##      7) elev>=7.552 68 25 2 (0.07352941 0.63235294 0.29411765)  
##       14) y>=330330 54 13 2 (0.09259259 0.75925926 0.14814815)  
##         28) elev< 9.947 49  9 2 (0.10204082 0.81632653 0.08163265) *
##         29) elev>=9.947 5  1 3 (0.00000000 0.20000000 0.80000000) *
##       15) y< 330330 14  2 3 (0.00000000 0.14285714 0.85714286) *

rpart.plot(m.ff.rpp, type=2, extra=1)

Q: Now how many splits? Which variables were used? Are all classes predicted? How pure are the leaves?

Q: How useful is this model? In other words, how well could flooding frequency class be mapped from knowing the distance to river and elevation?

Variable importance:

m.ff.rp$variable.importance

##     elev        y        x   dist.m 
## 38.26320 38.24051 16.58514 11.40992

m.ff.rpp$variable.importance

##        y     elev        x   dist.m 
## 36.74909 29.85313 14.99264  9.89731

In both the unpruned and pruned trees the y-coordinate and elevation are most important.

6.1 Sensitivity

Each RF is random, so several attempts will not be the same. How consistent are they?

Compute the RF several times and collect statistics:

n <- 24
rf.stats <- data.frame(rep=1:10, rsq=as.numeric(NA), mse=as.numeric(NA))
for (i in 1:n) {
  model.rf <- randomForest(logZn ~ ffreq + x + y + dist.m + elev + soil + lime,
                           data=meuse, importance=T, na.action=na.omit, mtry=5)
  summary(model.rf$rsq)
  summary(model.rf$mse)
  rf.stats[i, "mse"] <- median(summary(model.rf$mse))
  rf.stats[i, "rsq"] <- median(summary(model.rf$rsq))
}
summary(rf.stats[,2:3])

##       rsq              mse         
##  Min.   :0.7917   Min.   :0.01872  
##  1st Qu.:0.8002   1st Qu.:0.01921  
##  Median :0.8013   Median :0.01940  
##  Mean   :0.8013   Mean   :0.01940  
##  3rd Qu.:0.8033   3rd Qu.:0.01951  
##  Max.   :0.8083   Max.   :0.02034

hist(rf.stats[,"rsq"], xlab="RandomForest R^2")
rug(rf.stats[,"rsq"])

hist(rf.stats[,"mse"], xlab="RandomForest RMSE")
rug(rf.stats[,"mse"])

Q: How consistent are the 24 runs of the random forest?

6.2 Variable importance

Examine the variable importance. Type 1 is the increase in MSE if that predictor is assigned to its actuql point, rather than randomly assigned to a point. It is computed as follows (from the Help): “For each tree, the prediction error on the OOB portion of the data is recorded (…MSE for regression). Then the same is done after permuting each predictor variable. The difference between the two are then averaged over all trees, and normalized by the standard deviation of the differences”.

Type 2 is the increase in node purity. It is computed as follows: “The total decrease in node impurities from splitting on the [given] variable, averaged over all trees… measured by residual sum of squares.”

importance(m.lzn.rf, type=1)

##          %IncMSE
## ffreq  13.511952
## x      26.031058
## y      18.648671
## dist.m 57.516695
## elev   37.554231
## soil    6.934790
## lime    9.604052

importance(m.lzn.rf, type=2)

##        IncNodePurity
## ffreq      0.4501877
## x          1.1210360
## y          0.7244526
## dist.m     8.0261675
## elev       3.5421854
## soil       0.2220359
## lime       0.7047061

varImpPlot(m.lzn.rf, type=1)

Q: Does this agree with the variable importance from the single trees we built above?

Q: Why are some variables given importance, although they were not used in the regression tree models?

6.3 Examine variable importance and interaction with randomForestExplainer

Another way to look at this is by the depth in the tree at which each predictor is used. This is found with the min_depth_frame function of the randomForestExplainer package.

min_depth_frame <- min_depth_distribution(m.lzn.rf)
str(min_depth_frame)  # has results for all the trees

## 'data.frame':    3222 obs. of  3 variables:
##  $ tree         : int  1 1 1 1 1 1 2 2 2 2 ...
##  $ variable     : Factor w/ 7 levels "dist.m","elev",..: 1 2 3 4 5 6 1 2 4 5 ...
##  $ minimal_depth: num  0 1 1 4 2 3 1 0 1 2 ...
##  - attr(*, ".internal.selfref")=<externalptr>

plot_min_depth_distribution(min_depth_frame)

Q: Which predictor is used most as the root?

Q: Which predictor is used most at the second level?

Q: Which predictors are not used in any trees?

Q: Do these results agree with the variable importance plot?

This package gives a wider choice of importance measures:

importance_frame <- measure_importance(m.lzn.rf)
print(importance_frame)

##   variable mean_min_depth no_of_nodes mse_increase node_purity_increase
## 1   dist.m       0.342000        5437  0.081638354            8.0261675
## 2     elev       1.042000        6520  0.033787030            3.5421854
## 3    ffreq       3.625152         994  0.007941017            0.4501877
## 4     lime       4.404208         484  0.003762594            0.7047061
## 5     soil       5.189408         686  0.001803041            0.2220359
## 6        x       2.142000        5722  0.011567810            1.1210360
## 7        y       2.640000        5199  0.007384188            0.7244526
##   no_of_trees times_a_root       p_value
## 1         500          351 2.311138e-219
## 2         500          100  0.000000e+00
## 3         458            1  1.000000e+00
## 4         357           40  1.000000e+00
## 5         407            8  1.000000e+00
## 6         500            0 1.575661e-286
## 7         500            0 1.324765e-169

Q: Which predictor was used in the most nodes, over all trees?

Q: Which predictors were never used as the root of the tree?

Q: Which predictors were used in every tree?

The p-value measureswhether the observed number of successes (number of nodes in which \(X_j\) was used for splitting) exceeds the theoretical number of successes if they were random, following a binomial distribution. See https://cran.rstudio.com/web/packages/randomForestExplainer/vignettes/randomForestExplainer.html for details.

Another interesting plot is the “multiway importance plot”, which by default shows the number of times as root vs. the mean minimum depth for each predictor. There is typically a negative relation between times_a_root and mean_min_depth.

Q: Which predictor is most important by both measures?

Q: Is this a consistent negative relation? If not, where are the exceptions?

Another interesting plot is the pairwise comparison of different importance measures:

plot_importance_ggpairs(importance_frame)

Q: Which are the most and least related measures?

plot_multi_way_importance(m.lzn.rf, size_measure = "no_of_nodes")

Examine the interactions between predictors, in terms of their use within the forest:

interactions_frame <- min_depth_interactions(m.lzn.rf)
interactions_frame[order(interactions_frame$occurrences, decreasing = TRUE)[1:12], ]

##    variable root_variable mean_min_depth occurrences   interaction
## 43        y        dist.m      1.5595960         495      dist.m:y
## 8      elev        dist.m      0.3920364         494   dist.m:elev
## 1    dist.m        dist.m      1.0286343         492 dist.m:dist.m
## 36        x        dist.m      1.0994586         491      dist.m:x
## 9      elev          elev      1.1333818         483     elev:elev
## 44        y          elev      1.7066505         474        elev:y
## 37        x          elev      1.5122101         473        elev:x
## 2    dist.m          elev      1.2774303         456   elev:dist.m
## 15    ffreq        dist.m      3.9632323         395  dist.m:ffreq
## 48        y             x      1.9688081         389           x:y
## 41        x             x      2.2824646         388           x:x
## 13     elev             x      2.3496162         364        x:elev
##    uncond_mean_min_depth
## 43              2.640000
## 8               1.042000
## 1               0.342000
## 36              2.142000
## 9               1.042000
## 44              2.640000
## 37              2.142000
## 2               0.342000
## 15              3.625152
## 48              2.640000
## 41              2.142000
## 13              1.042000

plot_min_depth_interactions(interactions_frame)

Q: Which pair of predictors are most used near the root of the tree, and with the most co-occurrences?

“To further investigate the most frequent interaction we use the function plot_predict_interaction to plot the prediction of our forest on a grid of values for the components of each interaction. The function requires the forest, training data, variable to use on x and y-axis, respectively. In addition, one can also decrease the number of points in both dimensions of the grid from the default of 100 in case of insufficient memory using the parameter grid.”"

Here we see the predicted log(Zn) content for all combinations of distance from river and elevation:

plot_predict_interaction(m.lzn.rf, meuse, "dist.m", "elev")

6.4 Partial dependence plots

These plots show the effect of one predictor on the predictand, holding all the others constant at their median value. We show these in order of variable importance, as determined just above.

partialPlot(m.lzn.rf, pred.data=meuse, x.var="x")

partialPlot(m.lzn.rf, pred.data=meuse, x.var="y")

partialPlot(m.lzn.rf, pred.data=meuse, x.var="dist.m")

partialPlot(m.lzn.rf, pred.data=meuse, x.var="elev")

partialPlot(m.lzn.rf, pred.data=meuse, x.var="ffreq", which.class=1)

Q: What is the effect of just the distance to river on the RF prediction of log(Zn), when the other variables are held at their mean values (continuous variables, i.e., elev) or most common class (classified variables, i.e. ffreq, soil and lime). Is the effect monotonic (i.e., consistently increasing or decreasing)? Why or why not?

6.5 Goodness-of-fit and out-of-bag errors

Predict back onto the known points, and evaluate the goodness-of-fit:

p.rf <- predict(m.lzn.rf, newdata=meuse)
length(unique(p.rf))

## [1] 155

summary(r.rpp <- meuse$logZn - p.rf)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.1825439 -0.0364426 -0.0025826 -0.0004145  0.0329439  0.1577282

(rmse.rf <- sqrt(sum(r.rpp^2)/length(r.rpp)))

## [1] 0.06099567

plot(meuse$logZn ~ p.rf, asp=1, pch=20, xlab="fitted", ylab="actual", xlim=c(2,3.3),          ylim=c(2,3.3), main="log10(Zn), Meuse topsoils, Random Forest")
grid()
abline(0,1)

Q: What is the RMSE?

Q: How many different predictions were made? Why not just a few, as in regression trees?

A more realistic view of the predictive power of the model is the out-of-bag validation.

p.rf.oob <- predict(m.lzn.rf)
summary(r.rpp.oob <- meuse$logZn - p.rf.oob)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.426973 -0.082194 -0.010906 -0.001954  0.077214  0.377211

(rmse.oob <- sqrt(sum(r.rpp.oob^2)/length(r.rpp.oob)))

## [1] 0.1373769

plot(meuse$logZn ~ p.rf.oob, asp=1, pch=20,
     xlab="Out-of-bag cross-validation estimates",
     ylab="actual", xlim=c(2,3.3), ylim=c(2,3.3),
     main="log10(Zn), Meuse topsoils, Random Forest")
grid()
abline(0,1)

Q: What is the OOB RMSE? How does this compare with the fitted RMSE? Which is more realistic for predictions?

Q: How does this compare to the model fits?

Q: How does this compare with the RMSE from the Regression Tree?

6.6 Map over the grid with an RF model

We can try to build an RF model to map the grid, but again we can only use coordinates, flood frequency, soil type and relative distance to the river, because these are the only variables we have over the whole area, i.e., in the spatial object meuse.grid.

names(meuse.grid)

## [1] "x"      "y"      "part.a" "part.b" "dist"   "soil"   "ffreq"

m.lzn.rf.g <- randomForest(logZn ~ x + y + ffreq + dist + soil, data=meuse, importance=T, na.action=na.omit, mtry=3)
print(m.lzn.rf)

## 
## Call:
##  randomForest(formula = logZn ~ ffreq + x + y + dist.m + elev +      soil + lime, data = meuse, importance = T, mtry = 5, na.action = na.omit) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##           Mean of squared residuals: 0.01887242
##                     % Var explained: 80.67

(tmp <- importance(m.lzn.rf.g, type=1))

##        %IncMSE
## x     30.38927
## y     20.81755
## ffreq 37.19663
## dist  52.99835
## soil  11.68151

(100*tmp/sum(tmp))

##         %IncMSE
## x     19.851457
## y     13.598840
## ffreq 24.298292
## dist  34.620595
## soil   7.630817

varImpPlot(m.lzn.rf.g, type=1)

meuse.grid.sp$predictedLog10Zn.rf <- predict(m.lzn.rf.g, newdata=meuse.grid)
spplot(meuse.grid.sp, zcol="predictedLog10Zn.rf")

Q: Does this map look realistic? How does it compare with the map made with the single regressiom tree?

7.1 Unoptimized Cubist model

Now fit the model to the training set:

(c.model <- cubist(x = train.pred, y = train.resp))

## 
## Call:
## cubist.default(x = train.pred, y = train.resp)
## 
## Number of samples: 124 
## Number of predictors: 7 
## 
## Number of committees: 1 
## Number of rules: 2

summary(c.model)

## 
## Call:
## cubist.default(x = train.pred, y = train.resp)
## 
## 
## Cubist [Release 2.07 GPL Edition]  Sun Aug  4 18:41:04 2019
## ---------------------------------
## 
##     Target attribute `outcome'
## 
## Read 124 cases (8 attributes) from undefined.data
## 
## Model:
## 
##   Rule 1: [80 cases, mean 2.366983, range 2.053078 to 2.89098, est err 0.118638]
## 
##     if
##  dist.m > 140
##     then
##  outcome = -10.61724 + 0.000166 y - 0.000228 x - 0.111 elev
##            - 8e-05 dist.m
## 
##   Rule 2: [44 cases, mean 2.833035, range 2.281033 to 3.196176, est err 0.143328]
## 
##     if
##  dist.m <= 140
##     then
##  outcome = 3.586028 - 0.00129 dist.m - 0.088 elev
## 
## 
## Evaluation on training data (124 cases):
## 
##     Average  |error|           0.119965
##     Relative |error|               0.47
##     Correlation coefficient        0.87
## 
## 
##  Attribute usage:
##    Conds  Model
## 
##    100%   100%    dist.m
##           100%    elev
##            65%    x
##            65%    y
## 
## 
## Time: 0.0 secs

In this case there are only two leaves! each with a regression model. The split is on distance to river.

We predict onto the test set to get an estimate of predictive accuracy:

c.pred <- predict(c.model, test.pred)
## Test set RMSE
sqrt(mean((c.pred - test.resp)^2))

## [1] 0.1770856

## Test set R^2
cor(c.pred, test.resp)^2

## [1] 0.8244534

Now, fit a model on the entire set and then predict over the grid. Again we can only use predictors that are in the grid.

preds <- c("x","y","dist","soil","ffreq")
all.preds <- meuse[, preds]
all.resp <- meuse$logZn
c.model.grid <- cubist(x = all.preds, y = all.resp)
summary(c.model.grid)

## 
## Call:
## cubist.default(x = all.preds, y = all.resp)
## 
## 
## Cubist [Release 2.07 GPL Edition]  Sun Aug  4 18:41:04 2019
## ---------------------------------
## 
##     Target attribute `outcome'
## 
## Read 155 cases (6 attributes) from undefined.data
## 
## Model:
## 
##   Rule 1: [66 cases, mean 2.288309, range 2.053078 to 2.89098, est err 0.103603]
## 
##     if
##  x > 179095
##  dist > 0.211846
##     then
##  outcome = 2.406759 - 0.32 dist
## 
##   Rule 2: [9 cases, mean 2.596965, range 2.330414 to 2.832509, est err 0.116378]
## 
##     if
##  x <= 179095
##  dist > 0.211846
##     then
##  outcome = -277.415278 + 0.000847 y + 0.56 dist
## 
##   Rule 3: [80 cases, mean 2.772547, range 2.187521 to 3.264582, est err 0.157513]
## 
##     if
##  dist <= 0.211846
##     then
##  outcome = 2.632508 - 2.1 dist - 2.4e-05 x + 1.4e-05 y
## 
## 
## Evaluation on training data (155 cases):
## 
##     Average  |error|           0.111292
##     Relative |error|               0.41
##     Correlation coefficient        0.85
## 
## 
##  Attribute usage:
##    Conds  Model
## 
##    100%   100%    dist
##     48%    52%    x
##            57%    y
## 
## 
## Time: 0.0 secs

meuse.grid.sp$predictedLog10Zn.cubist <- predict(c.model.grid, newdata=meuse.grid[,preds])
spplot(meuse.grid.sp, zcol="predictedLog10Zn.cubist", col.regions=bpy.colors(64))

Here we can see the effect of the local regression models, as well as the hard split on x, y, and dist.

7.2 Optimizing cubist

Cubist can improve performance in two ways: (1) with “committees” and (2) adjusting by nearest neighbours:

The Cubist model can also use a boosting–like scheme called committees where iterative model trees are created in sequence. The first tree follows the procedure described in the last section. Subsequent trees are created using adjusted versions to the training set outcome: if the model over–predicted a value, the response is adjusted downward for the next model (and so on). Unlike traditional boosting, stage weights for each committee are not used to average the predictions from each model tree; the final prediction is a simple average of the predictions from each model tree.

Another innovation in Cubist using nearest–neighbors to adjust the predictions from the rule–based model. First, a model tree (with or without committees) is created. Once a sample is predicted by this model, Cubist can find it’s nearest neighbors and determine the average of these training set points.

To find the optimum values for these, we turn to the caret package.

The caret package, short for classification and regression training, contains numerous tools for developing predictive models using the rich set of models available in R. The package focuses on simplifying model training and tuning across a wide variety of modeling techniques.

See also: http://topepo.github.io/caret/model-training-and-tuning.html

The train function in caret can use many models, supplied by other packages, by specifying the method argument. Then the tuneGrid parameter refers to a grid of possible parameters to compare by cross-validation. Note we use all the points, because train does the splitting as it evaluates the alternative models.

require(caret)
test <- train(x = all.preds, y = all.resp, method="cubist", 
              tuneGrid = expand.grid(.committees = 1:12, 
                                     .neighbors = 0:5), 
              trControl = trainControl(method = 'cv'))
print(test)

## Cubist 
## 
## 155 samples
##   5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 140, 139, 140, 140, 139, 139, ... 
## Resampling results across tuning parameters:
## 
##   committees  neighbors  RMSE       Rsquared   MAE      
##    1          0          0.1883956  0.6644354  0.1458751
##    1          1          0.1588804  0.7350231  0.1166150
##    1          2          0.1460033  0.7801726  0.1091955
##    1          3          0.1465169  0.7847191  0.1106126
##    1          4          0.1493369  0.7795308  0.1132737
##    1          5          0.1520507  0.7743130  0.1169301
##    2          0          0.1683147  0.7411218  0.1315260
##    2          1          0.1608123  0.7397582  0.1172865
##    2          2          0.1462212  0.7862414  0.1088365
##    2          3          0.1448057  0.7933283  0.1096120
##    2          4          0.1457635  0.7923567  0.1110771
##    2          5          0.1458465  0.7932452  0.1128878
##    3          0          0.1676694  0.7274362  0.1300419
##    3          1          0.1594861  0.7355455  0.1159521
##    3          2          0.1463699  0.7792583  0.1092873
##    3          3          0.1450852  0.7868112  0.1097221
##    3          4          0.1463667  0.7846627  0.1116297
##    3          5          0.1471833  0.7842686  0.1135879
##    4          0          0.1722867  0.7230334  0.1350674
##    4          1          0.1595969  0.7374772  0.1157987
##    4          2          0.1465086  0.7796208  0.1087658
##    4          3          0.1457687  0.7860170  0.1102305
##    4          4          0.1471952  0.7840891  0.1124495
##    4          5          0.1479295  0.7841296  0.1146606
##    5          0          0.1692093  0.7264544  0.1308310
##    5          1          0.1573260  0.7432324  0.1134389
##    5          2          0.1451288  0.7830943  0.1085275
##    5          3          0.1445451  0.7892661  0.1097826
##    5          4          0.1458363  0.7875755  0.1119257
##    5          5          0.1465786  0.7876978  0.1137389
##    6          0          0.1681138  0.7348997  0.1306287
##    6          1          0.1577049  0.7448716  0.1139049
##    6          2          0.1453464  0.7845837  0.1083253
##    6          3          0.1450470  0.7898445  0.1096985
##    6          4          0.1465139  0.7878958  0.1119542
##    6          5          0.1472043  0.7882202  0.1139688
##    7          0          0.1689879  0.7319173  0.1298003
##    7          1          0.1590891  0.7424403  0.1141648
##    7          2          0.1478159  0.7790716  0.1101016
##    7          3          0.1472611  0.7850995  0.1114656
##    7          4          0.1486386  0.7831451  0.1129634
##    7          5          0.1494694  0.7834350  0.1150561
##    8          0          0.1692387  0.7303566  0.1308911
##    8          1          0.1583597  0.7451302  0.1140818
##    8          2          0.1469568  0.7816830  0.1096589
##    8          3          0.1467558  0.7864582  0.1110293
##    8          4          0.1482476  0.7842532  0.1130076
##    8          5          0.1490496  0.7844580  0.1148313
##    9          0          0.1680918  0.7334097  0.1293326
##    9          1          0.1580984  0.7445659  0.1135899
##    9          2          0.1460286  0.7835949  0.1082919
##    9          3          0.1456577  0.7891344  0.1096461
##    9          4          0.1470649  0.7872591  0.1117024
##    9          5          0.1473304  0.7886694  0.1135409
##   10          0          0.1667062  0.7361165  0.1293054
##   10          1          0.1583984  0.7437456  0.1141123
##   10          2          0.1463530  0.7827190  0.1089020
##   10          3          0.1461691  0.7873274  0.1105037
##   10          4          0.1475471  0.7854459  0.1122338
##   10          5          0.1484767  0.7850997  0.1142535
##   11          0          0.1667215  0.7397498  0.1287537
##   11          1          0.1583368  0.7437772  0.1138371
##   11          2          0.1460518  0.7836939  0.1080391
##   11          3          0.1456496  0.7891811  0.1092283
##   11          4          0.1468634  0.7879053  0.1108272
##   11          5          0.1471398  0.7893048  0.1128886
##   12          0          0.1645851  0.7456427  0.1277227
##   12          1          0.1583615  0.7432244  0.1137673
##   12          2          0.1456015  0.7842648  0.1080848
##   12          3          0.1451584  0.7898198  0.1092122
##   12          4          0.1462251  0.7892730  0.1110120
##   12          5          0.1469560  0.7895430  0.1129932
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were committees = 5 and neighbors = 3.

plot(test, metric="RMSE")

plot(test, metric="Rsquared")

Of course this is random in the splits, so each run will be different. In this case the graph shows clearly that more committees boost the performance considerably, from around 0.19 to 0.15 cross-validation RMSE. The nearest-neighbours also help considerably, up to about four, and then the improvement is less.

So we fit a model with these recommended parameters. Note the neighbour count is only used in prediction, not model building.

c.model.grid.2 <- cubist(x = all.preds, y = all.resp, committees=4)
summary(c.model.grid.2)

## 
## Call:
## cubist.default(x = all.preds, y = all.resp, committees = 4)
## 
## 
## Cubist [Release 2.07 GPL Edition]  Sun Aug  4 18:41:13 2019
## ---------------------------------
## 
##     Target attribute `outcome'
## 
## Read 155 cases (6 attributes) from undefined.data
## 
## Model 1:
## 
##   Rule 1/1: [66 cases, mean 2.288309, range 2.053078 to 2.89098, est err 0.103603]
## 
##     if
##  x > 179095
##  dist > 0.211846
##     then
##  outcome = 2.406759 - 0.32 dist
## 
##   Rule 1/2: [9 cases, mean 2.596965, range 2.330414 to 2.832509, est err 0.116378]
## 
##     if
##  x <= 179095
##  dist > 0.211846
##     then
##  outcome = -277.415278 + 0.000847 y + 0.56 dist
## 
##   Rule 1/3: [80 cases, mean 2.772547, range 2.187521 to 3.264582, est err 0.157513]
## 
##     if
##  dist <= 0.211846
##     then
##  outcome = 2.632508 - 2.1 dist - 2.4e-05 x + 1.4e-05 y
## 
## Model 2:
## 
##   Rule 2/1: [45 cases, mean 2.418724, range 2.10721 to 2.893762, est err 0.182228]
## 
##     if
##  x <= 179826
##  ffreq in {2, 3}
##     then
##  outcome = 128.701732 - 0.000705 x
## 
##   Rule 2/2: [121 cases, mean 2.443053, range 2.053078 to 3.055378, est err 0.181513]
## 
##     if
##  dist > 0.0703468
##     then
##  outcome = 30.512065 - 0.87 dist - 0.000154 x
## 
##   Rule 2/3: [55 cases, mean 2.543648, range 2.075547 to 3.055378, est err 0.125950]
## 
##     if
##  dist > 0.0703468
##  ffreq = 1
##     then
##  outcome = 37.730889 - 0.000314 x - 0.35 dist + 6.5e-05 y
## 
##   Rule 2/4: [34 cases, mean 2.958686, range 2.574031 to 3.264582, est err 0.139639]
## 
##     if
##  dist <= 0.0703468
##     then
##  outcome = 2.982852 - 0.36 dist
## 
## Model 3:
## 
##   Rule 3/1: [75 cases, mean 2.325347, range 2.053078 to 2.89098, est err 0.150262]
## 
##     if
##  dist > 0.211846
##     then
##  outcome = 2.263175 - 0.11 dist
## 
##   Rule 3/2: [80 cases, mean 2.772547, range 2.187521 to 3.264582, est err 0.156918]
## 
##     if
##  dist <= 0.211846
##     then
##  outcome = 2.978074 - 2.42 dist
## 
## Model 4:
## 
##   Rule 4/1: [121 cases, mean 2.443053, range 2.053078 to 3.055378, est err 0.161020]
## 
##     if
##  dist > 0.0703468
##     then
##  outcome = 4.601694 - 0.61 dist - 0.000101 x + 4.9e-05 y
## 
##   Rule 4/2: [51 cases, mean 2.505994, range 2.075547 to 3.055378, est err 0.196863]
## 
##     if
##  x <= 179731
##  dist > 0.0703468
##     then
##  outcome = 175.779227 - 0.000967 x
## 
##   Rule 4/3: [55 cases, mean 2.543648, range 2.075547 to 3.055378, est err 0.147504]
## 
##     if
##  dist > 0.0703468
##  ffreq = 1
##     then
##  outcome = 53.629078 - 0.000343 x - 0.44 dist + 3.3e-05 y
## 
##   Rule 4/4: [34 cases, mean 2.958686, range 2.574031 to 3.264582, est err 0.140965]
## 
##     if
##  dist <= 0.0703468
##     then
##  outcome = 2.975269 - 0.08 dist
## 
## 
## Evaluation on training data (155 cases):
## 
##     Average  |error|           0.110135
##     Relative |error|               0.40
##     Correlation coefficient        0.84
## 
## 
##  Attribute usage:
##    Conds  Model
## 
##     95%    88%    dist
##     21%    64%    x
##     19%           ffreq
##            39%    y
## 
## 
## Time: 0.0 secs

meuse.grid.sp$predictedLog10Zn.cubist.2 <- 
  predict(c.model.grid.2, newdata=meuse.grid[,preds],
          neighbors=5)
summary(meuse.grid.sp$predictedLog10Zn.cubist.2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.954   2.228   2.346   2.445   2.640   3.154

spplot(meuse.grid.sp, zcol="predictedLog10Zn.cubist.2", col.regions=bpy.colors(64),
       main="Optimized Cubist prediction")

A much nicer looking map.

Variable importance. The usage field of the fitted model shows the usage of each variable in either the rule conditions or the (terminal) linear model.

c.model.grid$usage

##   Conditions Model Variable
## 1        100   100     dist
## 2         48    52        x
## 3          0    57        y
## 4          0     0     soil
## 5          0     0    ffreq

c.model.grid.2$usage

##   Conditions Model Variable
## 1         95    88     dist
## 2         21    64        x
## 3         19     0    ffreq
## 4          0    39        y
## 5          0     0     soil

caret::varImp(c.model.grid)

##       Overall
## dist    100.0
## x        50.0
## y        28.5
## soil      0.0
## ffreq     0.0

caret::varImp(c.model.grid.2)

##       Overall
## dist     91.5
## x        42.5
## ffreq     9.5
## y        19.5
## soil      0.0

Distance is always used in the single model, but sometimes not in the committee model.

Notice how flood frequency came into the model when a committee was formed, i.e., subsidiary rules to adjust mis-matches from higher-level rules.

It is interesting to see how the different models fit a single predictor, in this case distance to river:

distances <- all.preds[, "dist", drop = FALSE]
rules_only <- cubist(x = distances, y = all.resp)
rules_and_com <- cubist(x = distances, y = all.resp, committees = 4)

predictions <- distances
predictions$true <- all.resp
predictions$rules <- predict(rules_only, distances, neighbors = 0)
predictions$rules_neigh <- predict(rules_only, distances, neighbors = 5)
predictions$committees <- predict(rules_and_com, distances, neighbors = 0)
predictions$committees_neigh <- predict(rules_and_com, distances, neighbors = 5)
summary(predictions)

##       dist              true           rules        rules_neigh   
##  Min.   :0.00000   Min.   :2.053   Min.   :2.069   Min.   :2.046  
##  1st Qu.:0.07569   1st Qu.:2.297   1st Qu.:2.311   1st Qu.:2.331  
##  Median :0.21184   Median :2.513   Median :2.482   Median :2.541  
##  Mean   :0.24002   Mean   :2.556   Mean   :2.545   Mean   :2.556  
##  3rd Qu.:0.36407   3rd Qu.:2.829   3rd Qu.:2.804   3rd Qu.:2.791  
##  Max.   :0.88039   Max.   :3.265   Max.   :2.982   Max.   :2.992  
##    committees    committees_neigh
##  Min.   :2.043   Min.   :2.035   
##  1st Qu.:2.347   1st Qu.:2.330   
##  Median :2.495   Median :2.524   
##  Mean   :2.545   Mean   :2.555   
##  3rd Qu.:2.771   3rd Qu.:2.791   
##  Max.   :2.924   Max.   :2.990

plot(true ~ dist, data=predictions, pch=20, col="darkgray",
     xlab="covariate: normalized distance to river",
     ylab="log10(Zn)", main="Cubist models")
grid()
m <- matrix(cbind(predictions$dist, predictions$rules), ncol=2)
ix <- order(m[,1]); m <- m[ix,]
lines(m, col=1, lwd=2)
m <- matrix(cbind(predictions$dist, predictions$rules_neigh), ncol=2)
ix <- order(m[,1]); m <- m[ix,]
lines(m, col=2, lwd=2)
m <- matrix(cbind(predictions$dist, predictions$committees), ncol=2)
ix <- order(m[,1]); m <- m[ix,]
lines(m, col=3, lwd=2)
m <- matrix(cbind(predictions$dist, predictions$committees_neigh), ncol=2)
ix <- order(m[,1]); m <- m[ix,]
lines(m, col=4, lwd=2)
legend("topright", c("rules", "rules + neighbour",
                     "committee", "committee + neighbour"),
       pch="----", col=1:4)

Notice how the committee (green) smoothes out the centre of distance range from the single-rules model (black).

Notice how the neighbours come much closer to some of the points – too close?

9.1 Method 1 – buffers to quantiles

library(ggplot2)
g <- ggplot(meuse, mapping=aes(x=logZn))
g + geom_histogram(bins=16, fill="lightgreen", color="blue") + geom_rug()

(q.lzn <- quantile(meuse$logZn, seq(0,1,by=0.0625)))
# must include.lowest to have lowest value in some class
classes.lzn.q <- cut(meuse$logZn, breaks=q.lzn, ordered_result=TRUE, include.lowest=TRUE)
levels(classes.lzn.q)

##       0%    6.25%    12.5%   18.75%      25%   31.25%    37.5%   43.75% 
## 2.053078 2.139461 2.201372 2.261554 2.296665 2.322219 2.362643 2.414125 
##      50%   56.25%    62.5%   68.75%      75%   81.25%    87.5%   93.75% 
## 2.513218 2.603414 2.680704 2.754247 2.828981 2.873820 2.920515 3.056094 
##     100% 
## 3.264582 
##  [1] "[2.05,2.14]" "(2.14,2.2]"  "(2.2,2.26]"  "(2.26,2.3]"  "(2.3,2.32]" 
##  [6] "(2.32,2.36]" "(2.36,2.41]" "(2.41,2.51]" "(2.51,2.6]"  "(2.6,2.68]" 
## [11] "(2.68,2.75]" "(2.75,2.83]" "(2.83,2.87]" "(2.87,2.92]" "(2.92,3.06]"
## [16] "(3.06,3.26]"

coordinates(meuse) <- ~x + y
grid.dist.lzn <- GSIF::buffer.dist(meuse["logZn"], meuse.grid.sp, classes.lzn.q)
plot(stack(grid.dist.lzn))

buffer.dists <- over(meuse, grid.dist.lzn)
str(buffer.dists)

## 'data.frame':    155 obs. of  16 variables:
##  $ layer.1 : num  1897 1809 1865 1878 1737 ...
##  $ layer.2 : num  1784 1695 1755 1772 1633 ...
##  $ layer.3 : num  735 680 621 561 402 ...
##  $ layer.4 : num  468 412 362 330 179 ...
##  $ layer.5 : num  804 730 721 699 544 ...
##  $ layer.6 : num  685 612 601 582 431 ...
##  $ layer.7 : num  268 283 126 0 160 ...
##  $ layer.8 : num  369 330 233 160 0 ...
##  $ layer.9 : num  268 226 160 170 126 ...
##  $ layer.10: num  243 160 226 305 283 ...
##  $ layer.11: num  369 283 344 400 330 ...
##  $ layer.12: num  144 160 0 126 233 ...
##  $ layer.13: num  734 645 738 794 700 ...
##  $ layer.14: num  1531 1442 1537 1585 1471 ...
##  $ layer.15: num  0 89.4 144.2 268.3 368.8 ...
##  $ layer.16: num  89.4 0 160 282.8 344.1 ...

meuse@data <- cbind(meuse@data, buffer.dists)

buffer.dists[1,]

##    layer.1  layer.2  layer.3 layer.4 layer.5  layer.6  layer.7  layer.8
## 1 1897.367 1783.928 735.3911 468.188  803.99 684.6897 268.3282 368.7818
##    layer.9 layer.10 layer.11 layer.12 layer.13 layer.14 layer.15 layer.16
## 1 268.3282 243.3105 368.7818 144.2221 734.3024 1531.013        0 89.44272

dn <- paste(names(grid.dist.lzn), collapse="+")
(fm <- as.formula(paste("logZn ~", dn)))

## logZn ~ layer.1 + layer.2 + layer.3 + layer.4 + layer.5 + layer.6 + 
##     layer.7 + layer.8 + layer.9 + layer.10 + layer.11 + layer.12 + 
##     layer.13 + layer.14 + layer.15 + layer.16

(rf <- randomForest(fm , meuse@data, importance=TRUE, min.split=6, mtry=6, ntree=640))

## 
## Call:
##  randomForest(formula = fm, data = meuse@data, importance = TRUE,      min.split = 6, mtry = 6, ntree = 640) 
##                Type of random forest: regression
##                      Number of trees: 640
## No. of variables tried at each split: 6
## 
##           Mean of squared residuals: 0.02514371
##                     % Var explained: 74.25

pred.rf <- predict(rf, newdata=meuse@data)

plot(rf)

importance(rf, type=1)

##            %IncMSE
## layer.1  15.125917
## layer.2  14.959521
## layer.3  13.808773
## layer.4  11.403248
## layer.5  11.637966
## layer.6  17.236908
## layer.7   7.091440
## layer.8  10.915838
## layer.9   9.779422
## layer.10 10.792435
## layer.11  6.995092
## layer.12 17.543799
## layer.13 11.360035
## layer.14 13.343180
## layer.15 18.142337
## layer.16 23.833595

varImpPlot(rf, type=1)

plot(meuse$logZn ~ pred.rf, asp=1, pch=20, xlab="Random forest fit", ylab="Actual value", main="Zn, log(mg kg-1)")
abline(0,1); grid()

(rmse.rf <- sqrt(sum((pred.rf-meuse$logZn)^2)/length(pred.rf)))

## [1] 0.06472746

#
oob.rf <- predict(rf)
plot(meuse$logZn ~ oob.rf, asp=1, pch=20, xlab="Out-of-bag prediction", ylab="Actual value", main="Zn, log(mg kg-1)")
abline(0,1); grid()

(rmse.oob <- sqrt(sum((oob.rf-meuse$logZn)^2)/length(oob.rf)))

## [1] 0.1585677

pred.grid <- predict(rf, newdata=grid.dist.lzn@data)
summary(pred.grid)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.111   2.316   2.437   2.489   2.650   3.169

meuse.grid.sp$logZn.rf <- pred.grid
breaks <- seq(2, 3.3, by=.05)
p1 <- spplot(meuse.grid.sp, zcol="logZn.rf", main="Zn, log(mg kg-1)", sub="RRF 16 distance buffers", at=breaks)
print(p1)

require(gstat)
v <- variogram(logZn ~ 1, meuse)
(vmf <- fit.variogram(v, model=vgm(0.1, "Sph", 1000, 0.02)))

##   model       psill    range
## 1   Nug 0.009554751   0.0000
## 2   Sph 0.111394769 896.9909

plot(v, model=vmf, plot.numbers=T)

k.cv <- krige.cv(logZn ~ 1, meuse, model=vmf, verbose=FALSE)
plot(meuse$logZn ~ k.cv$var1.pred, pch=20, xlab="Cross-validation fit", ylab="Actual value", main="Zn, log(mg kg-1)", asp=1)
abline(0,1); grid()

(rmse.kcv <- sqrt(sum((k.cv$var1.pred-meuse$logZn)^2)/length(k.cv$var1.pred)))

## [1] 0.170157

k.ok <- krige(logZn ~ 1, meuse, newdata=meuse.grid.sp, model=vmf)

## [using ordinary kriging]

summary(k.ok$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.074   2.275   2.420   2.479   2.680   3.231

p2 <- spplot(k.ok, zcol="var1.pred", main="Zn, log(mg kg-1)", sub="Ordinary Kriging", at=breaks)
print(p2, split=c(1,1,2,1), more=T)
print(p1, split=c(2,1,2,1), more=F)

k.ok$diff.srf.ok <- meuse.grid.sp$logZn.rf - k.ok$var1.pred
summary(k.ok$diff.srf.ok)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.237562 -0.040157  0.004313  0.010131  0.048444  0.403515

p3 <- spplot(k.ok, zcol="diff.srf.ok", main="Spatial RF - OK predictions, log(Zn)", col.regions=topo.colors(36))
print(p3)

9.2 Method 2 – buffers around each point

The concept of quantile buffers is similar to the empirical variogram, i.e., assume stationarity of spatial correlation structure.

Another approach is to use the distance from each known point as a predictor. Obviously this increases the size of the problem. This is the approach taken in the Hengl paper cited above.

Set up the buffer distances and visualize a few of them:

grid.dist.0 <- GSIF::buffer.dist(meuse["logZn"], meuse.grid.sp, as.factor(1:nrow(meuse)))
dim(grid.dist.0)

## [1] 3103  155

summary(grid.dist.0$layer.1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    1677    2685    2556    3506    4569

plot(stack(grid.dist.0[1:9]))

Extract the distances to each buffer at each point, and add these as fields to the points. Note that we first have to remove the previous buffer distances (the 16 quantiles).

buffer.dists <- over(meuse, grid.dist.0)
dim(buffer.dists)

## [1] 155 155

# distances from point 1 to all the other points
(buffer.dists[1,1:10])

##   layer.1  layer.2 layer.3  layer.4  layer.5  layer.6  layer.7  layer.8
## 1       0 89.44272 144.222 268.3282 368.7818 481.6638 268.3282 243.3105
##   layer.9 layer.10
## 1     400  468.188

names(meuse@data)

##  [1] "cadmium"  "copper"   "lead"     "zinc"     "elev"     "dist"    
##  [7] "om"       "ffreq"    "soil"     "lime"     "landuse"  "dist.m"  
## [13] "logZn"    "logCu"    "logPb"    "logCd"    "layer.1"  "layer.2" 
## [19] "layer.3"  "layer.4"  "layer.5"  "layer.6"  "layer.7"  "layer.8" 
## [25] "layer.9"  "layer.10" "layer.11" "layer.12" "layer.13" "layer.14"
## [31] "layer.15" "layer.16"

meuse@data <- meuse@data[,1:16]
meuse@data <- cbind(meuse@data, buffer.dists)

So now we have 155 buffer distances. Set up a formula for the RF and build the RF. Note that the default mtry will be 52, i.e., distance to 1/3 of the points. We increase this to 2/3.

dn0 <- paste(names(grid.dist.lzn), collapse="+") 
fm0 <- as.formula(paste("logZn ~ ", dn0))    
(rf <- randomForest(fm0 , meuse@data, importance=TRUE, 
                    mtry=102, min.split=6, ntree=640))

## 
## Call:
##  randomForest(formula = fm0, data = meuse@data, importance = TRUE,      mtry = 102, min.split = 6, ntree = 640) 
##                Type of random forest: regression
##                      Number of trees: 640
## No. of variables tried at each split: 16
## 
##           Mean of squared residuals: 0.09413569
##                     % Var explained: 3.6

This does not perform as well as using the quantiles, the \(R^2\) is about 10% lower.

Mapping:

pred.grid <- predict(rf, newdata=grid.dist.lzn@data)
summary(pred.grid)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.299   2.477   2.596   2.590   2.697   2.897

meuse.grid.sp$logZn.rf0 <- pred.grid
breaks <- seq(2, 3.3, by=.05)
p2 <- spplot(meuse.grid.sp, zcol="logZn.rf0", main="Zn, log(mg kg-1)",
             sub="RF distance all points", at=breaks)
# compare with 16-quantile buffers
print(p2, split=c(1,1,2,1), more=T)
print(p1, split=c(2,1,2,1), more=F)

A clearly inferior map. So we don’t care about individual point values, rather about the distance to quantiles.

10.1 Mapping classes

Again we need to restrict the RF model to variables also available in the prediction grid.

(m.ffreq.rf.2 <- randomForest(ffreq ~ x + y + dist + soil, data=meuse, importance=T, na.action=na.omit, mtry=3))

## 
## Call:
##  randomForest(formula = ffreq ~ x + y + dist + soil, data = meuse,      importance = T, mtry = 3, na.action = na.omit) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 28.39%
## Confusion matrix:
##    1  2  3 class.error
## 1 71 11  2   0.1547619
## 2 14 29  5   0.3958333
## 3  3  9 11   0.5217391

This is less accurate than the model with elevation. Anyway, we show the map:

meuse.grid.sp$predictedffreq.rf <- predict(m.ffreq.rf.2, newdata=meuse.grid)
spplot(meuse.grid.sp, zcol="predictedffreq.rf", main="Most probable flood frequency class")

Not a nice-looking map. It disagrees strongly with the actual flooding-frequency map.

summary(meuse.grid.sp$ffreq.tf <-
          as.factor((meuse.grid.sp$predictedffreq.rf==
                       meuse.grid.sp$ffreq)))

## FALSE  TRUE 
##  1452  1651

spplot(meuse.grid.sp, zcol="ffreq.tf",
       main="Accuracy of RF model")

Probability of class 1 at each grid cell:

meuse.grid.sp$probffreq <-
  predict(m.ffreq.rf.2, newdata=meuse.grid, type="prob")
meuse.grid.sp$probffreq.rf.1 <- meuse.grid.sp$probffreq[,1]
meuse.grid.sp$probffreq.rf.2 <- meuse.grid.sp$probffreq[,2]
meuse.grid.sp$probffreq.rf.3 <- meuse.grid.sp$probffreq[,3]
p1 <- spplot(meuse.grid.sp, zcol="probffreq.rf.1",
             main="Prob. class 1")
p2 <- spplot(meuse.grid.sp, zcol="probffreq.rf.2",
             main="Prob. class 2")
p3 <- spplot(meuse.grid.sp, zcol="probffreq.rf.3",
             main="Prob. class 3")
print(p1, split=c(1,1,3,1), more=T)
print(p2, split=c(2,1,3,1), more=T)
print(p3, split=c(3,1,3,1), more=F)

11.1 Mapping

We first have to add the covariates to the prediction grid.

grid.dist.lzn.covar <- grid.dist.lzn
grid.dist.lzn.covar$dist <- meuse.grid$dist
grid.dist.lzn.covar$ffreq <- meuse.grid$ffreq
grid.dist.lzn.covar$soil <- meuse.grid$soil

Now predict over the grid and display:

pred.grid.c <- predict(rf.c, newdata=grid.dist.lzn.covar@data)
summary(pred.grid.c)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.280   2.345   2.422   2.491   2.615   2.919

meuse.grid.sp$logZn.rf.c <- pred.grid.c
breaks <- seq(2, 3.3, by=.05)
p1.c <- spplot(meuse.grid.sp, zcol="logZn.rf.c", main="Zn, log(mg kg-1)", sub="Random forest on distance buffers & covariates", at=breaks)
print(p1.c, split=c(1,1,2,1), more=T)
print(p1, split=c(2,1,2,1), more=F)

Compute and display the differences:

meuse.grid.sp$logZn.rf.diff <- meuse.grid.sp$logZn.rf.c -  meuse.grid.sp$logZn.rf
summary(meuse.grid.sp$logZn.rf.diff)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.4434319 -0.0854825  0.0005576  0.0022438  0.0907812  0.4409751

p3.c <- spplot(meuse.grid.sp, zcol="logZn.rf.diff", main="Spatial RF with - without covariates, log(Zn)", col.regions=topo.colors(36))
print(p3.c)

Q: Compare this to the maps made by RF/Spatial Distances without covariates. Where does using the covariates result in higher and lower predicted values? Try to explain why.

Challenge: Map logZn over the grid by Kriging with External Drift (KED), using the same three covariates (dist, ffreq, soil), and compare with the map made by spatial random forest with covariates. Reminder: the variogram must be computed for the linear model residuals using these covariates.

Q: Describe the differences between the residual and ordinary variogram.

Q: Describe the differences between the OK and KED maps.

Q: Compare the KED map with the spatial random forest + covariates map.