Main.Rmd

---
title: "Organic carbon in surface sediments - Main"
output: html_notebook
---

# Install packages

```{r packages, message=FALSE}
rm(list=ls())

library(terra)
library(dplyr)
library(caret)
library(ggplot2)
library(sf)
library(CAST)
library(lwgeom)
library(geosphere)
library(quantregForest)
library(doParallel)
library(ModelMetrics)
library(forcats)
```


# Preparation

## Which sediment depth interval?

0: 0 - 10 cm

```{r depth_interval}
# upper limit
du <- 0

#lower limit
dl <- 10
```


## Load required data

If more than one sediment depth interval is predicted, it might be useful to use the prediction of the uppermost layer as a predictor.
*Note: The name of the OC prediction of the sediment interval above needs to be added manually.*

```{r load_data}
predictors <- rast("input/predictors.tif")
AoI <- read_sf("input/AoI.shp")

if(du != 0){
  OCabove <- rast(paste0("output/OC0-10cm_median_", date, ".tif"))
  predictors <- c(predictors, OCabove)
  names(predictors)[[length(names(predictors))]] <- "OC"
  rm(OCabove)
} 

OC <- read_sf(paste0("input/OC", du, ".shp"))
OC <- vect(OC)

names(predictors)
```


## Type of response?

Define which property is used as response data.

```{r response_type}
resp_type <- "TOC"
resp_unit <- "weight-%"
```


## Extract predictors

```{r extract_predictors}
OC <- terra::extract(predictors, OC, bind = TRUE)
OC <- terra::na.omit(OC, field = "", geom = TRUE)

plot(predictors$BATH)
plot(OC, pch = 20, col = "black", cex = 0.5, add = TRUE)
```


## Create a regression matrix

```{r regression matrix}
rm_oc <- as.data.frame(OC)
rm_oc <- rm_oc[-1]
names(rm_oc)[1] <- "TOC_perc"

summary(rm_oc)
```


# Data exploration

## Histogram

```{r hist_oc_content}
hist(rm_oc$TOC_perc, breaks = 40, main = "", xlab = paste0(resp_type, " (", resp_unit, ")"))
```


## Distances in environmental space

Distances in environmental (feature) space are computed.

```{r env_space_dist, message=FALSE}
dist_env <- plot_geodist(st_as_sf(OC), predictors,
                     type = "feature",
                     showPlot = FALSE)

dist_env$plot
dist_env$plot + scale_x_log10()
```


## Distances in geographic space

Distances in geographic space are computed.

```{r geogr_space_dist, message=FALSE}
dist_geogr <- plot_geodist(st_as_sf(OC), predictors,
                     type = "geo",
                     unit="km",
                     showPlot = FALSE)

dist_geogr$plot
dist_geogr$plot + scale_x_log10()
```


# Quantile Regression Forest model

## Creating spatial blocks

Spatial blocks and folds are created. The folds will be used in a spatial k-fold cross validation. The k-fold nearest neighbour distance matching algorithm is used here.

```{r nndm}
k <- 10 # Number of folds
knndmfolds <- knndm(tpoints = st_as_sf(OC),
                    modeldomain = AoI,
                    k = k,
                    samplesize = 2000)
```


## Distances in geographic space including CV distances

```{r geogr_space_dist2, message=FALSE}
dist_geogr2 <- plot_geodist(st_as_sf(OC), predictors,
                     cvfolds= knndmfolds$indx_test,
                     type = "geo",
                     unit="km",
                     showPlot = FALSE)

dist_geogr2$plot
dist_geogr2$plot + scale_x_log10()
```


## Model tuning

A Quantile Regression Forest model is tuned. Predictor variables are selected in a forward feature selection approach and various values of the mtry parameter are tested in a spatial k-fold cross validation.

The maximum number of iterations to be performed can be calculated upfront, based on the number of pre-selected predictors:

```{r max_iter}
factorial(length(names(predictors)))/(factorial(2)*factorial(length(names(predictors))-2)) + sum(c((length(names(predictors))-2):1))
```


### Forward feature selection

The best combination of predictor variables (features) is found in a forward feature selection process.

```{r ffs, message=FALSE, warning=FALSE}
nCores <- detectCores()
cl <- makePSOCKcluster(nCores - 1)
registerDoParallel(cl)

set.seed(42)

model <- ffs(rm_oc[names(predictors)],
               rm_oc$TOC_perc,
               metric = "Rsquared",
               method="qrf",
               what = 0.5,
               replace = FALSE,
               importance = TRUE,
               trControl = trainControl(method="CV",
                                        number = k,
                                        savePredictions = "final",
                                        index = knndmfolds$indx_train, 
                                        allowParallel = TRUE),
               verbose = TRUE)

stopCluster(cl)

model

sel_preds <- model$selectedvars
```


### FFS plot

Plot of R2 over the model runs.

```{r ffs_plot}
plot_ffs(model)
```


## Validation statistics

The validation results of the optimal RF model.

Note that these are the statistics based on the predicted values of the selected model. These differ from the values from the tuning (above), which are the means of the k predictions based on the folds.

```{r validation_stats}
t <- data.frame(model$pred$pred, model$pred$obs)

validation <- data.frame(mse=numeric(), rmse=numeric(), r2=numeric())
validation[1,1] <- round(sum(t$model.pred.obs - t$model.pred.pred)/nrow(t), 3)
validation[1,2] <- round(rmse(t$model.pred.obs, t$model.pred.pred), 3)
validation[1,3] <- round(cor(t$model.pred.obs, t$model.pred.pred)^2, 3)

colnames(validation) <- c("ME", "RMSE", "r2")
rownames(validation) <- NULL
validation
```


## Validation plot

```{r validation_plot, message=FALSE}
ggplot(t, aes(x = model.pred.pred, y = model.pred.obs)) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_abline(intercept = 0, slope = 1, colour = "grey", linewidth = 1.2) +
  theme_bw() +
  scale_x_continuous(name = "Predicted value") +
  scale_y_continuous(name = "Observed value") +
  ggtitle(paste0(resp_type, " (", resp_unit, ")"))
```


## Variable importance

```{r variable_importance_plot, warning=FALSE}
imp <- varImp(model$finalModel)
imp$Predictor <- rownames(imp)
rownames(imp) <- NULL
imp <- imp[order(imp[1], decreasing = TRUE), c(2, 1)]
colnames(imp)[2] <- "IncMSE"
imp

impfig <- imp %>%
  mutate(Predictor = fct_reorder(Predictor, IncMSE)) %>%
  ggplot( aes(x=Predictor, y=IncMSE)) +
  geom_bar(stat="identity", fill="#f68060", alpha=.6, width=.4) +
  coord_flip() +
  xlab("") +
  ylab("% increase in MSE") +
  ggtitle("Organic carbon content") +
  theme_bw()
    
impfig

jpeg(filename = "output/TOC_VarImp.jpg", width = 12, height = 12, units = "cm", res = 500)
impfig
dev.off()
```


## Distances in environmental space including CV distances

```{r env_space_dist2, message=FALSE}
dist_env2 <- plot_geodist(st_as_sf(OC), predictors,
                     type = "feature",
                     cvfolds= knndmfolds$indx_test,
                     variables = sel_preds,
                     showPlot = FALSE)

dist_env2$plot
dist_env2$plot + scale_x_log10()
```


## Partial dependence

Partial dependence plots give a graphical depiction of the marginal effect of a variable on the response.

```{r partial_plots}
m2 <- model$finalModel
class(m2) <- "randomForest"

for (i in 1:length(sel_preds)) {
  partialPlot(x = m2, pred.data = rm_oc, x.var = sel_preds[i], main = "", xlab = sel_preds[i], ylab = paste0(resp_type, " (", resp_unit, ")"))
}

for (i in 1:length(sel_preds)) {
  jpeg(filename = paste0("output/", sel_preds[i], ".jpg"), width = 12, height = 12, units = "cm", res = 500)
  par(mar = c(4.5,5,1,1))
  partialPlot(x = m2, pred.data = rm_oc, x.var = sel_preds[i], ylim=c(0.8,1.75), main = "", xlab = sel_preds[i], ylab = "Organic carbon content (weight-%)")
  dev.off()
}
```


# Predict QRF model

## Predict OC content

Organic carbon content is predicted. Median values of the QRF distribution are calculated as central values. The 90% prediction interval and the prediction interval ratio are calculated as measures of uncertainty.

```{r predict_oc}
preds <- raster::stack(predictors[[sel_preds]])
OC_med <- predict(preds, model$finalModel, what = 0.5)
OC_p95 <- predict(preds, model$finalModel, what = 0.95)
OC_p5 <- predict(preds, model$finalModel, what = 0.05)
OC_pi90 <- OC_p95 - OC_p5
OC_pir <- OC_pi90 / OC_med

hist(OC_med, breaks = 20, main = "", xlab = paste0(resp_type, " (", resp_unit, ")"))
```


## Area of applicability

```{r aoa}
OC_trainDI <- trainDI(model = model,
                        variables = sel_preds)
print(OC_trainDI)

OC_aoa <- aoa(newdata = predictors, 
                model = model,
                trainDI = OC_trainDI,
                variables = sel_preds,
)

plot(OC_aoa)
```


## Plot results

```{r plot_results}
plot(OC_med, main = "OC median")
plot(OC_pi90, main = "90% prediction interval")
plot(OC_pir, main = "Prediction interval ratio")
plot(OC_aoa$DI, main = "Dissimilarity index")
plot(OC_aoa$AOA, main = "Area of applicability")

fr <- freq(OC_aoa$AOA)
print(paste0("AOA = ", round(100*fr$count[2]/ sum(fr$count),2), "% of pixels"))
```


## Convert AOA from raster to polygon

```{r aoa_poly}
aoa_poly <- as.polygons(OC_aoa$AOA, dissolve = TRUE)
plot(aoa_poly)

write_sf(st_as_sf(aoa_poly), dsn = "output", layer = paste0("OC", du, "-", dl, "cm_AOA_", Sys.Date()), driver = "ESRI Shapefile")
```


## Export results

```{r export_results}
writeRaster(OC_med, paste0("output/OC", du, "-", dl, "cm_median_", Sys.Date(), ".tif"))
writeRaster(OC_p5, paste0("output/OC", du, "-", dl, "cm_P5_", Sys.Date(), ".tif"))
writeRaster(OC_p95, paste0("output/OC", du, "-", dl, "cm_P95_", Sys.Date(), ".tif"))
writeRaster(OC_pi90, paste0("output/OC", du, "-", dl, "cm_PI90_", Sys.Date(), ".tif"))
writeRaster(OC_pir, paste0("output/OC", du, "-", dl, "cm_PIR_", Sys.Date(), ".tif"))
#writeRaster(OC_aoa$DI, paste0("output/OC", du, "-", dl, "cm__DI_", Sys.Date(), ".tif"))
writeRaster(OC_aoa$AOA, paste0("output/OC", du, "-", dl, "cm_AOA_", Sys.Date(), ".tif"))
```


## Output a log file

```{r log}
sink(file = paste0("output/ModelLog_", du, "-", dl, "_", Sys.Date(), ".txt"))
print("Selected Predictors")
sel_preds
model
print("Final Model")
paste0("ME = ", validation[1,1])
paste0("RMSE = ", validation[1,2])
paste0("R2 = ", validation[1,3])
paste0("AOA = ", round(100*fr$count[2]/ sum(fr$count),2), "% of pixels")
sink()
```


# Finishing off

# Save QRF model

```{r save_model}
saveRDS(model, "qrfmodel.rds")
```


## Save session info

```{r save_session_info}
sessionInfo <- sessionInfo()
save(sessionInfo, file = "sessionInfo_main.Rdata")
rm("sessionInfo")
```


## Save global environment

```{r save_global_env}
save.image(file = "globEnv_main.RData")
```