Dissever Tutorial

Malone, B. and Roudier, P.

2018-01-15

library(dissever)

library(raster)
library(viridis)
library(dplyr)
library(magrittr)

Here is the coarse model to dissever:

plot(edgeroi$carbon, col = viridis(100))

and here are the covariates that will be used for the disseveration:

plot(edgeroi$predictors, col = viridis(100))

Example

In this section, we will dissever the sample dataset with the available environmental covariates, using 4 different regression methods:

It is simple to switch between regression methods using the method option of dissever.

But first let’s setup some parameters:

min_iter <- 5 # Minimum number of iterations
max_iter <- 10 # Maximum number of iterations
p_train <- 0.025 # Subsampling of the initial data

We can then launch the 4 different disseveration procedures:

# Random Forest
res_rf <- dissever(
    coarse = edgeroi$carbon, # stack of fine resolution covariates
    fine = edgeroi$predictors, # coarse resolution raster
    method = "rf", # regression method used for disseveration
    p = p_train, # proportion of pixels sampled for training regression model
    min_iter = min_iter, # minimum iterations
    max_iter = max_iter # maximum iterations
  )

# Cubist
res_cubist <- dissever(
  coarse = edgeroi$carbon, 
  fine = edgeroi$predictors, 
  method = "cubist", 
  p = p_train, 
  min_iter = min_iter,
  max_iter = max_iter
)

# GAM
res_gam <- dissever(
  coarse = edgeroi$carbon, 
  fine = edgeroi$predictors, 
  method = "gamSpline", 
  p = p_train, 
  min_iter = min_iter,
  max_iter = max_iter
)

# Linear model
res_lm <- dissever(
  coarse = edgeroi$carbon, 
  fine = edgeroi$predictors, 
  method = "lm", 
  p = p_train, 
  min_iter = min_iter,
  max_iter = max_iter
)

Analysing the results

The object returned by dissever is an object of class dissever. It’s basically a list with 3 elements: * fit: a train object storing the regression model used in the final disseveration * map: a RasterLayer object storing the dissevered map * perf: a data.frame object storing the evolution of the RMSE (and confidence intervals) with the disseveration iterations.

The dissever result has a plot function that can plot either the map or the performance results, depending on the type option.

Maps

Below are the dissevered maps for all 4 regression methods:

# Plotting maps
par(mfrow = c(2, 2))
plot(res_rf, type = 'map', main = "Random Forest")
plot(res_cubist, type = 'map', main = "Cubist")
plot(res_gam, type = 'map', main = "GAM")
plot(res_lm, type = 'map', main = "Linear Model")

Convergence

We can also analyse and plot the optimisation of the dissevering model:

# Plot performance
par(mfrow = c(2, 2))
plot(res_rf, type = 'perf', main = "Random Forest")
plot(res_cubist, type = 'perf', main = "Cubist")
plot(res_gam, type = 'perf', main = "GAM")
plot(res_lm, type = 'perf', main = "Linear Model")

Performance of the regression

The tools provided by the caret package can also be leveraged, e.g. to plot the observed vs. predicted values:

# Plot preds vs obs
preds <- extractPrediction(list(res_rf$fit, res_cubist$fit, res_gam$fit, res_lm$fit))
plotObsVsPred(preds)

The models can be compared using the preds object that we just derived using the extractPrediction function. In this case I’m computing the R2, the RMSE and the CCC for each model:

# Compare models
perf <- preds %>% 
  group_by(model, dataType) %>% 
  summarise(
    rsq = cor(obs, pred)^2, 
    rmse = sqrt(mean((pred - obs)^2))
  )
perf
## # A tibble: 4 x 4
## # Groups:   model [?]
##   model     dataType   rsq  rmse
##   <fctr>    <fctr>   <dbl> <dbl>
## 1 cubist    Training 0.259  5.44
## 2 gamSpline Training 0.398  4.33
## 3 lm        Training 0.398  4.36
## 4 rf        Training 0.939  1.58

Ensemble approach

We can either take the best option (in this case Random Forest – but the results probably show that we should add some kind of validation process), or use a simple ensemble-type approach. For the latter, I am computing weights for each of the maps based on the CCC of the 4 different regression models:

# We can weight results with Rsquared
w <- perf$rsq / sum(perf$rsq)

# Make stack of weighted predictions and compute sum
l_maps <- list(res_cubist$map, res_gam$map, res_lm$map, res_rf$map)
ens <- lapply(1:4, function(x) l_maps[[x]] * w[x]) %>%
  stack %>%
  sum

Ensemble modelling seem to work better when the models are not too correlated – i.e. when they capture different facets of the data:

s_res <- stack(l_maps)
names(s_res) <- c('Cubist', 'GAM', 'Linear Model', 'Random Forest')
s_res %>%
  as.data.frame %>%
  na.exclude %>%
  cor
##                  Cubist       GAM Linear.Model Random.Forest
## Cubist        1.0000000 0.6854560    0.6435265     0.6997470
## GAM           0.6854560 1.0000000    0.9392300     0.8466576
## Linear.Model  0.6435265 0.9392300    1.0000000     0.8228443
## Random.Forest 0.6997470 0.8466576    0.8228443     1.0000000

Here is the resulting map:

# Plot result
plot(ens, col = viridis(100), main = "CCC Weighted Average")