inst/doc/cast04-AOA-tutorial.R

In CAST: 'caret' Applications for Spatial-Temporal Models

## ----setup, echo=FALSE--------------------------------------------------------
knitr::opts_chunk$set(fig.width = 8.83)

## ----message = FALSE, warning=FALSE-------------------------------------------
library(CAST)
library(caret)
library(terra)
library(sf)
library(viridis)
library(gridExtra)

## ----message = FALSE,include=FALSE, warning=FALSE-----------------------------
RMSE = function(a, b){
    sqrt(mean((a - b)^2,na.rm=T))
}

## ----message = FALSE, warning=FALSE-------------------------------------------
predictors <- rast(system.file("extdata","bioclim.tif",package="CAST"))
plot(predictors,col=viridis(100))

## ----message = FALSE, warning=FALSE-------------------------------------------

generate_random_response <- function(raster, predictornames =
names(raster), seed = sample(seq(1000), 1)){
  operands_1 = c("+", "-", "*", "/")
  operands_2 = c("^1","^2")
  
  expression <- paste(as.character(predictornames, sep=""))
  # assign random power to predictors
  set.seed(seed)
  expression <- paste(expression,
                      sample(operands_2, length(predictornames),
replace = TRUE),
                      sep = "")
  
  # assign random math function between predictors (expect after the last one)
  set.seed(seed)
  expression[-length(expression)] <- paste(expression[-
length(expression)],
                                           sample(operands_1,
length(predictornames)-1, replace = TRUE),
                                           sep = " ")
  print(paste0(expression, collapse = " "))
  # collapse
  e = paste0("raster$", expression, collapse = " ")
  
  response = eval(parse(text = e))
  names(response) <- "response"
  return(response)
  
}

## ----message = FALSE, warning=FALSE-------------------------------------------
response <- generate_random_response (predictors, seed = 10)
plot(response,col=viridis(100),main="virtual response")

## ----message = FALSE, warning=FALSE-------------------------------------------
mask <- predictors[[1]]
values(mask)[!is.na(values(mask))] <- 1
mask <- st_as_sf(as.polygons(mask))
mask <- st_make_valid(mask)

## ----message = FALSE, warning=FALSE-------------------------------------------
set.seed(15)
samplepoints <- st_as_sf(st_sample(mask,20,"random"))

plot(response,col=viridis(100))
plot(samplepoints,col="red",add=T,pch=3)

## ----message = FALSE, warning=FALSE-------------------------------------------
trainDat <- extract(predictors,samplepoints,na.rm=FALSE)
trainDat$response <- extract(response,samplepoints,na.rm=FALSE, ID=FALSE)$response
trainDat <- na.omit(trainDat)

## ----message = FALSE, warning=FALSE-------------------------------------------
set.seed(10)
model <- train(trainDat[,names(predictors)],
               trainDat$response,
               method="rf",
               importance=TRUE,
               trControl = trainControl(method="cv"))
print(model)


## ----message = FALSE, warning=FALSE-------------------------------------------
plot(varImp(model,scale = F),col="black")

## ----message = FALSE, warning=FALSE-------------------------------------------
prediction <- predict(predictors,model,na.rm=T)
truediff <- abs(prediction-response)
plot(rast(list(prediction,response)),main=c("prediction","reference"))

## ----message = FALSE, warning=FALSE-------------------------------------------
AOA <- aoa(predictors, model, LPD = TRUE, verbose = FALSE)
class(AOA)
names(AOA)
print(AOA)

## ----message = FALSE, warning=FALSE-------------------------------------------
plot(AOA)

## ----message = FALSE, warning=FALSE,  fig.show="hold", out.width="50%"--------
plot(truediff,col=viridis(100),main="true prediction error")
plot(AOA$DI,col=viridis(100),main="DI")
plot(AOA$LPD,col=viridis(100),main="LPD")
plot(prediction, col=viridis(100),main="prediction for AOA")
plot(AOA$AOA,col=c("grey","transparent"),add=T,plg=list(x="topleft",box.col="black",bty="o",title="AOA"))

## ----message = FALSE, warning=FALSE-------------------------------------------
set.seed(25)
samplepoints <- clustered_sample(mask,75,15,radius=25000)

plot(response,col=viridis(100))
plot(samplepoints,col="red",add=T,pch=3)


## ----message = FALSE, warning=FALSE-------------------------------------------

trainDat <- extract(predictors,samplepoints,na.rm=FALSE)
trainDat$response <- extract(response,samplepoints,na.rm=FALSE)$response
trainDat <- data.frame(trainDat,samplepoints)
trainDat <- na.omit(trainDat)

## ----message = FALSE, warning=FALSE-------------------------------------------
set.seed(10)
model_random <- train(trainDat[,names(predictors)],
               trainDat$response,
               method="rf",
               importance=TRUE,
               trControl = trainControl(method="cv"))
prediction_random <- predict(predictors,model_random,na.rm=TRUE)
print(model_random)

## ----message = FALSE, warning=FALSE-------------------------------------------
folds <- CreateSpacetimeFolds(trainDat, spacevar="parent",k=10)
set.seed(15)
model <- train(trainDat[,names(predictors)],
                 trainDat$response,
                     method="rf",
                 importance=TRUE,
                 tuneGrid = expand.grid(mtry = c(2:length(names(predictors)))),
                 trControl = trainControl(method="cv",index=folds$index))
  print(model)
  
prediction <- predict(predictors,model,na.rm=TRUE)

## ----message = FALSE, warning=FALSE-------------------------------------------
AOA_spatial <- aoa(predictors, model, LPD = TRUE, verbose = FALSE)

AOA_random <- aoa(predictors, model_random, LPD = FALSE, verbose = FALSE)

## ----message = FALSE, warning=FALSE,  fig.show="hold", out.width="50%"--------
plot(AOA_spatial$DI,col=viridis(100),main="DI")
plot(AOA_spatial$LPD,col=viridis(100),main="LPD")
plot(prediction, col=viridis(100),main="prediction for AOA \n(spatial CV error applies)")
plot(AOA_spatial$AOA,col=c("grey","transparent"),add=TRUE,plg=list(x="topleft",box.col="black",bty="o",bg = 'white',title="AOA"))
plot(prediction_random, col=viridis(100),main="prediction for AOA \n(random CV error applies)")
plot(AOA_random$AOA,col=c("grey","transparent"),add=TRUE,plg=list(x="topleft",box.col="black",bty="o",bg = 'white',title="AOA"))

## ----message = FALSE, warning=FALSE-------------------------------------------
grid.arrange(plot(AOA_spatial, variable = "DI") + ggplot2::ggtitle("Spatial CV"),
             plot(AOA_random, variable = "DI") + ggplot2::ggtitle("Random CV"), ncol = 2)

## ----message = FALSE, warning=FALSE-------------------------------------------
###for the spatial CV:
RMSE(values(prediction)[values(AOA_spatial$AOA)==1],
     values(response)[values(AOA_spatial$AOA)==1])
RMSE(values(prediction)[values(AOA_spatial$AOA)==0],
     values(response)[values(AOA_spatial$AOA)==0])
model$results

###and for the random CV:
RMSE(values(prediction_random)[values(AOA_random$AOA)==1],
     values(response)[values(AOA_random$AOA)==1])
RMSE(values(prediction_random)[values(AOA_random$AOA)==0],
     values(response)[values(AOA_random$AOA)==0])
model_random$results

## ----message = FALSE, warning=FALSE-------------------------------------------
DI_RMSE_relation <- errorProfiles(model, AOA_spatial$parameters, multiCV=TRUE,
                                    window.size = 5, length.out = 5, variable = "DI")
plot(DI_RMSE_relation)

LPD_RMSE_relation <- errorProfiles(model, AOA_spatial$parameters, multiCV=TRUE,
                                    window.size = 5, length.out = 5, variable = "LPD")
plot(LPD_RMSE_relation)

DI_expected_RMSE = terra::predict(AOA_spatial$DI, DI_RMSE_relation)
LPD_expected_RMSE = terra::predict(AOA_spatial$LPD, LPD_RMSE_relation)

# account for multiCV changing the DI threshold
DI_updated_AOA = AOA_spatial$DI > attr(DI_RMSE_relation, "AOA_threshold")

# account for multiCV changing the DI threshold
LPD_updated_AOA = AOA_spatial$DI > attr(LPD_RMSE_relation, "AOA_threshold")



plot(DI_expected_RMSE,col=viridis(100),main="DI expected RMSE")
plot(DI_updated_AOA, col=c("grey","transparent"),add=TRUE,plg=list(x="topleft",box.col="black",bty="o",bg = 'white',title="AOA"))

plot(LPD_expected_RMSE,col=viridis(100),main="LPD expected RMSE")
plot(LPD_updated_AOA, col=c("grey","transparent"),add=TRUE,plg=list(x="topleft",box.col="black",bty="o",bg = 'white',title="AOA"))

## ----message = FALSE, warning=FALSE-------------------------------------------
data(cookfarm)
# calculate average of VW for each sampling site:
dat <- aggregate(cookfarm[,c("VW","Easting","Northing")],by=list(as.character(cookfarm$SOURCEID)),mean)
# create sf object from the data:
pts <- st_as_sf(dat,coords=c("Easting","Northing"))

##### Extract Predictors for the locations of the sampling points
studyArea <- rast(system.file("extdata","predictors_2012-03-25.tif",package="CAST"))
st_crs(pts) <- crs(studyArea)
trainDat <- extract(studyArea,pts,na.rm=FALSE)
pts$ID <- 1:nrow(pts)
trainDat <- merge(trainDat,pts,by.x="ID",by.y="ID")
# The final training dataset with potential predictors and VW:
head(trainDat)

## ----message = FALSE, warning=FALSE-------------------------------------------
predictors <- c("DEM","NDRE.Sd","TWI","Bt")
response <- "VW"

model <- train(trainDat[,predictors],trainDat[,response],
               method="rf",tuneLength=3,importance=TRUE,
               trControl=trainControl(method="LOOCV"))
model

## ----message = FALSE, warning=FALSE-------------------------------------------
#Predictors:
plot(stretch(studyArea[[predictors]]))

#prediction:
prediction <- predict(studyArea,model,na.rm=TRUE)

## ----message = FALSE, warning=FALSE,  fig.show="hold", out.width="30%"--------
AOA <- aoa(studyArea, model, LPD = TRUE, verbose = FALSE)

#### Plot results:
plot(AOA$DI,col=viridis(100),main="DI with sampling locations (red)")
plot(pts,zcol="ID",col="red",add=TRUE)

plot(AOA$LPD,col=viridis(100),main="LPD with sampling locations (red)")
plot(pts,zcol="ID",col="red",add=TRUE)

plot(prediction, col=viridis(100),main="prediction for AOA \n(LOOCV error applies)")
plot(AOA$AOA,col=c("grey","transparent"),add=TRUE,plg=list(x="topleft",box.col="black",bty="o", bg = 'white', title="AOA"))