---
title: "Mapping classes"
author: "D G Rossiter"
date: "`r Sys.Date()`"
output:
   html_document:
     toc: TRUE
     toc_float: TRUE  
     theme: "lumen"
     code_folding: show
     number_section: TRUE
     fig_height: 4
     fig_width: 6
     fig_align: 'center'
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=5, fig.height=5, fig.align='center', fig.path='rmdfigs/', warning=FALSE, message=FALSE)
```

This tutorial shows how to map class probabilities over a grid and then make a map of the most probable, along with a map of prediction reliability, represented by the maximum probability of any class.

It also has a section on using classification trees and random forests to classify.

# Dataset

Load the example dataset; this is in the `gstat` package:

```{r}
require(sp)
require(gstat)
data(jura) # loads several objects, see ?jura
summary(jura.pred)
# convert to spatial objects by identifying which fields are coordinates
coordinates(jura.pred) <- ~Xloc + Yloc
coordinates(jura.grid) <- ~Xloc + Yloc
spplot(jura.pred, zcol="Landuse", key.space="right")
```

# Modeling spatial dependence

We will try to predict land use `Landuse` over the grid by OK, separately for each class.
We use indicator kriging (IK) despite its theoretical drawbacks and somewhat (!) _ad hoc_ nature, just to get some probabilities as examples.

References:

* Bierkens, M. F. P., & Burrough, P. A. (1993a). The Indicator Approach to Categorical Soil Data .1. Theory. Journal of Soil Science, 44(2), 361–368.
* Bierkens, M. F. P., & Burrough, P. A. (1993b). The indicator approach to categorical soil data. II. Application to mapping and land use suitability analysis. Journal of Soil Science, 44(2), 369–381.
* Goovaerts, P. (1997). Geostatistics for natural resources evaluation. New York; Oxford: Oxford University Press.

Make separate indicators for each class:

```{r}
table(jura.pred$Landuse)
jura.pred$Forest <- (jura.pred$Landuse=="Forest")
table(jura.pred$Forest)
jura.pred$Pasture <- (jura.pred$Landuse=="Pasture")
table(jura.pred$Pasture)
jura.pred$Meadow <- (jura.pred$Landuse=="Meadow")
table(jura.pred$Meadow)
jura.pred$Tillage <- (jura.pred$Landuse=="Tillage")
table(jura.pred$Tillage)
```

Per-class empirical variograms:

```{r fig.width=6, fig.height=4, fig.align="left", out.width='50%', fig.show="hold"}
v.forest <- variogram(Forest ~ 1, data=jura.pred)
plot(v.forest, pl=T) 
v.pasture <- variogram(Pasture ~ 1, data=jura.pred)
plot(v.pasture, pl=T)
v.meadow <- variogram(Meadow ~ 1, data=jura.pred)
plot(v.meadow, pl=T)
v.tillage <- variogram(Tillage ~ 1, data=jura.pred)
plot(v.tillage, pl=T)
```

None of these have really clear structures, but anyway we model them, to illustrate the methods.

```{r fig.width=6, fig.height=4, fig.align="left", out.width='50%',fig.show="hold"}
(vmf.forest <- fit.variogram(v.forest, model=vgm(0.15, "Exp", 1, 0)))
plot(v.forest, pl=T, model=vmf.forest) 
(vmf.pasture <- fit.variogram(v.pasture, model=vgm(0.15, "Exp", 1, 0)))
plot(v.pasture, pl=T, model=vmf.pasture) 
(vmf.meadow <- fit.variogram(v.meadow, model=vgm(0.15, "Exp", 1, 0)))
plot(v.meadow, pl=T, model=vmf.meadow) 
(vmf.tillage <- fit.variogram(v.tillage, model=vgm(0.15, "Exp", 1, 0)))
plot(v.tillage, pl=T, model=vmf.tillage) 
```

# Mapping individual classes

Now use the variograms to krige over the grid; this gives the probability of each class at each location.

```{r}
ik.forest <- krige(Forest ~ 1, loc=jura.pred, model=vmf.forest, newdata=jura.grid)
ik.pasture <- krige(Pasture ~ 1, loc=jura.pred, model=vmf.pasture, newdata=jura.grid)
ik.meadow <- krige(Meadow ~ 1, loc=jura.pred, model=vmf.meadow, newdata=jura.grid)
ik.tillage <- krige(Tillage ~ 1, loc=jura.pred, model=vmf.tillage, newdata=jura.grid)
```

Kriging is not a convex predictor, so it is possible to get values outside of the data range, in this case $[0,1]$. So limit predictions outside this range to the extremes:


```{r}
ik.forest$var1.pred <- pmax(ik.forest$var1.pred, 0)
ik.forest$var1.pred <- pmin(ik.forest$var1.pred, 1)
ik.pasture$var1.pred <- pmax(ik.pasture$var1.pred, 0)
ik.pasture$var1.pred <- pmin(ik.pasture$var1.pred, 1)
ik.meadow$var1.pred <- pmax(ik.meadow$var1.pred, 0)
ik.meadow$var1.pred <- pmin(ik.meadow$var1.pred, 1)
ik.tillage$var1.pred <- pmax(ik.tillage$var1.pred, 0)
ik.tillage$var1.pred <- pmin(ik.tillage$var1.pred, 1)
```

Now plot these:


```{r}
ramp <- seq(0,1,by=0.1)
spplot(ik.forest, zcol="var1.pred", key.space="right", main="Prob(Forest)",
       cuts=ramp, col.regions=rev(terrain.colors(64)))
spplot(ik.pasture, zcol="var1.pred", key.space="right", main="Prob(Pasture)",
       cuts=ramp, col.regions=rev(terrain.colors(64)))
spplot(ik.meadow, zcol="var1.pred", key.space="right", main="Prob(Meadow)",
       cuts=ramp, col.regions=rev(terrain.colors(64)))
spplot(ik.tillage, zcol="var1.pred", key.space="right", main="Prob(Tillage)",
       cuts=ramp, col.regions=rev(terrain.colors(64)))
```

Clearly meadow is highly probable in most of the area, and tillage hardly over any area.


# Most probable class

To find the most-probable value, combine the predictions into one dataframe:

```{r}
ik <- data.frame(p.forest=ik.forest$var1.pred,
                 p.pasture=ik.pasture$var1.pred,
                 p.meadow=ik.meadow$var1.pred,
                 p.tillage=ik.tillage$var1.pred)
summary(ik)
```

Convert to a raster stack, i.e., a stack of individual layers. First each data frame has to be marked as `gridded`, so it becomes a `SpatialPixelsDataFrame`, and then converted to raster layer. The first variable, i.e., the prediction `var1.pred`, will be the raster layer value.

```{r}
require(raster)
gridded(ik.forest) <- TRUE
r.forest <- raster(ik.forest)
# image(r.forest, asp=1)
gridded(ik.pasture) <- TRUE
r.pasture <- raster(ik.pasture)
# image(r.pasture, asp=1)
gridded(ik.meadow) <- TRUE
r.meadow <- raster(ik.meadow)
# image(r.meadow, asp=1)
gridded(ik.tillage) <- TRUE
r.tillage <- raster(ik.tillage)
# image(r.tillage, asp=1)
```

Combine predictions into one object (both `raster` and `sp`) and display:

```{r out.width='80%'}
(r <- raster::stack(r.tillage, r.pasture, r.meadow, r.tillage))
names(r) <- c("prob.forest", "prob.pasture", "prob.meadow", "prob.tillage")
r.spdf <- as(r, "SpatialPixelsDataFrame")
spplot(r.spdf)
```

Clearly, meadow dominates, and there are no high probabilities of forest or tillage.

Now we can use the `which.max` function of the `raster` library on this stack to find the most probable class at each grid cell:

```{r}
mpc <- raster::which.max(r)
table(mpc@data@values)
mpc.df <- as.data.frame(mpc)
table(mpc.df)
```

Note that tillage is never the most probable class.

Now make a map of the most probable class:

```{r}
# image(mpc, asp=1, col=1:4)
# better graphics using sp, not raster
mpc.spdf <- as(mpc, "SpatialPixelsDataFrame")
str(mpc.spdf)
mpc.spdf@data$layer <- as.factor(mpc.spdf@data$layer)
levels(mpc.spdf$layer)  <- c("Forest", "Pasture", "Meadow")
spplot(mpc.spdf, col.regions=c("darkgreen", "lightblue", "green"))
```

# Maximum classification probability

We can also make a map of the maximum probability of _any_ class, to show where we are more or less confident in the result.

```{r fig.width=7, fig.height=5}
# apply the `max` function along the rows of the matrix
max.prob <- apply(r.spdf@data, 1, max)
summary(max.prob)
hist(max.prob) # quite some locations with low probability of any class
```

```{r}
(m <- max(r))
m.spdf <- as(m, "SpatialPixelsDataFrame")
spplot(m.spdf, key.space="right", main="Maximum class probability")
```

The areas with low probability (e.g., right-centre) are not at all reliably predicted.

# Uncertainty

The individual probabilities can be combined into one measure of uncertainty, for example the Shannon entropy $-\sum p_i \cdot \log_2 p_i$.

However, this assumes that $\sum_i p_i = 1$, which is not the case here since each class was predicted separately. The solution is to normalize all sums to 1, proportionally.

We first see the total probability, of all four classes, for the prediction grid:

```{r}
sum.prob <- apply(r.spdf@data, 1, sum)
summary(sum.prob)
r.spdf$sum.prob <- sum.prob
spplot(r.spdf, zcol="sum.prob", key.space="right", main="Sum of class probabilities")
```

Some areas are much more likely to be classified into any class than others. Some are over-determined, i.e., several classes have high probability.

Now normalize per grid cell:

```{r}
# see how this works on the first six cells
str(r.spdf@data)
# normalize by the sum
(test <- r.spdf@data[1:6,1:4]/r.spdf@data[1:6,5])
# the probabilities should sum to 1
apply(test, 1, sum)
# yes, so apply to all cells
summary(normalized <- r.spdf@data[,1:4]/r.spdf@data[,5])
# function to compute Shannon entropy from a vector
compute.shannon <- function(v) {
  # zero probabilities do not enter into the calculation
  ix <- which(v==0)
  if (length(ix > 0)) { v <- v[-ix] }
  return(-sum(v * log2(v)))
}
## avoid 0's
summary(shannon.entropy <- apply(normalized, 1, compute.shannon))
```

Display on a map:

```{r}
r.spdf$shannon <- shannon.entropy
spplot(r.spdf, zcol="shannon", key.space="right", main="Shannon entropy", col.regions=topo.colors(64))
```

We see some areas quite certain (low entropy), mostly in the meadows -- these had high individual probability of meadow and low for the others. The largest confusion (high entropy) is in the right-centre where no class had high probability.

# Classification tree

Somewhat an artificial example, but how well can we predict the land use from the metals, position, and rock type?

```{r}
library(rpart)
library(rpart.plot)
```

Bring the metric coordinates back into the prediction frame:

```{r}
jura.df <- as.data.frame(jura.pred)
```

## Prior probabilities = observations

Classification tree:

Note "The default priors are proportional to the data counts". This can be changed with `parms=`, optional parameters for the splitting function.

In this case Meadow is over half the observations, Tillage very few.

```{r}
table(jura.df$Landuse)
```


"The splitting index can be `gini` [default] or `information`".

Here we specify `cp=0` to build as large a tree as possible, but limit the classification bins to 3 observations each.

```{r fig.width=6, fig.height=3}
ct.lu <- rpart(Landuse ~ Xloc + Yloc + Rock + Cd + Co + Cr + Cu + Ni + Pb + Zn,
               data=jura.df,
               method="class", model=TRUE,
               minbucket=3, cp=0)
printcp(ct.lu)
plotcp(ct.lu)
```

Clearly over-fit, prune back to 12 splits, at `cp=0.01` approx.

```{r}
ct.lu <- prune(ct.lu, cp=0.01)
```

```{r fig.width=12, fig.height=9}
cp.table.class <- ct.lu[["cptable"]]
# total variance explained is sum of the CP for all split
paste("Proportion of variance explained:", round(sum(cp.table.class[,"CP"]),3))
par(mfrow=c(1,2))
rpart.plot(ct.lu, type=4, extra=1)
rpart.plot(ct.lu, type=4, extra=4)
par(mfrow=c(1,1))
```

See what is predicted by this model:

```{r}
ct.lu.pred <- predict(ct.lu, newdata=jura.df, type="class")
table(ct.lu.pred)
table(jura.df$Landuse)
```

Notice that tillage is never predicted; the predicted proportions are close to the observed.

Misclassifications:

```{r}
table(jura.df$Landuse, ct.lu.pred)
```



Details of the tree:

```{r}
summary(ct.lu)
```

## Equal prior probabilities

Suppose this is a biased sample (it isn't but...) and we do not want to use the observed proportions as priors. So set them equal. Use the same complexity parameter as found with the pruning, above.

```{r fig.width=6, fig.height=3}
(p <- nlevels(jura.df$Landuse))
ct.lu.eq <- rpart(Landuse ~ Xloc + Yloc + Rock + Cd + Co + Cr + Cu + Ni + Pb + Zn,
               data=jura.df,
               method="class", model=TRUE,
               minbucket=3, cp=0.01,
               parms=list(prior=rep(1/p, p)))
printcp(ct.lu.eq)
plotcp(ct.lu.eq)
```

```{r fig.width=12, fig.height=9}
cp.table.class <- ct.lu.eq[["cptable"]]
# total variance explained is sum of the CP for all split
paste("Proportion of variance explained:", round(sum(cp.table.class[,"CP"]),3))
par(mfrow=c(1,2))
rpart.plot(ct.lu.eq, type=4, extra=1)
rpart.plot(ct.lu.eq, type=4, extra=4)
par(mfrow=c(1,1))
```

```{r}
summary(ct.lu.eq)
```

See what is predicted by this model:

```{r}
ct.lu.eq.pred <- predict(ct.lu.eq, newdata=jura.df, type="class")
table(ct.lu.eq.pred)
table(jura.df$Landuse)
```

Now tillage is over-predicted as are forest and pasture, at the expense of meadow. These do not agree with the observations.

Misclassifications:

```{r}
table(jura.df$Landuse, ct.lu.eq.pred)
```

All actual tillage are correctly predicted, but in addition five meadow and one pasture also are predicted as tillage. A large number of meadows are predicted as forest and pasture.

This shows the effect of specifying prior probabilities.

# Random forest

```{r}
library(randomForest)
library(randomForestExplainer)
```

## Default weighting

Accept default `mtry` = $\lfloor{\sqrt{p}\rfloor}$, in this case `3`; do not specify class weights `classwt`, so all are equally-weighted. These are equal priors but implicitly weighted by the number of observations, because the random selection for each tree will more or less reproduce that proportion.



```{r}
rf.lu <- randomForest(Landuse ~ Xloc + Yloc + Rock + Cd + Co + Cr + Cu + Ni + Pb + Zn,
               data=jura.df, ntree=1024,
               importance=TRUE)
print(rf.lu)
```

Again, meadow is over-predicted and the others under-predicted; tillage is never predicted.

Variable importance:

```{r fig.width=10, fig.height=6}
varImpPlot(rf.lu)
```

Cu is important, it picks out several of the meadows.

The Gini impurity Gini impurity is how often a randomly chosen observation *would* be incorrectly classified by the forest *if* it were randomly classified according to the proportion of observations. So here the original unbalance of the observations is influential.

Cross-validation errors as a function of the number of trees:

```{r fig.width=8, fig.height=4}
plot(rf.lu)
palette()[1:5]
```

Overall OOB error is the black line; then the classes are in order forest (red), pasture (green), meadow (blue) -- the lowest by far; tillage (cyan) -- complete error.


If more variables are tried at each split, the result is closer to a single classification tree. Show this with 6 out of the 10, rather than 3:

```{r fig.width=10, fig.height=6}
rf.lu.6 <- randomForest(Landuse ~ 
                          Xloc + Yloc + Rock + Cd + Co + Cr + Cu + Ni + Pb + Zn,
               data=jura.df, ntree=1024, mtry=6,
               importance=TRUE)
print(rf.lu.6)
varImpPlot(rf.lu.6)
```

Error rate is higher; Cu becomes even more important.

Cross-validation errors as a function of the number of trees:

```{r fig.width=8, fig.height=4}
plot(rf.lu.6)
```

## Equal priors

For random forest we have to specify class weights; to make these essentially equal priors we inverse-weight by the number of observations.

```{r}
(weights <- 1/as.numeric(table(jura.df$Landuse)))
rf.lu.eq <- randomForest(Landuse ~
                           Xloc + Yloc + Rock + Cd + Co + Cr + Cu + Ni + Pb + Zn,
               data=jura.df, ntree=1024,
               classwt=weights,
               importance=TRUE)
print(rf.lu.eq)
```

This has little effect, in fact it results in more observations being classified as tillage, which seems counter-intuitive.

```{r}
table(predict(rf.lu))
table(predict(rf.lu.eq))
```