# How to validate or make predictions from a spatial Cox PH survival model built with R-INLA?

I'm trying to model deforestation as a survival analysis. I have a raster map where unaffected areas are zero and deforested pixels have values 1-20 depending on the year of deforestation (2001-20). I take 10km square grids as 'individuals' and 'death' is where the number of deforested pixels reaches a certain level. The other grid squares are 'survivors' (i.e. there is no such thing as 'lost to follow-up'). I have 4 covariates: difficulty of access, mean elevation, mean population, and land cover type.

Adding a spatial term to the regression seems to make a big difference - not surprising, as forest tends to be cleared near to already-cleared areas. What I'd like to do is to build a model from the data to 2010 and use that to see how well it predicts deforestation from 2010-20. There have been various posts on this subject, but they tend to involve the survival and rms packages. Can anyone suggest where to start, given an INLA object? I know Cox models only output hazard rations ... My code was:

    library(INLA)
library(rgdal)
library(rgeos)
library(sp)
library(raster)
library(spdep)

# read in grids with coordinates, covariates, def_count = count of 'deforested' pixels
# def_mean = mean of non-zero 'deforested' pixel values (= year of deforestation) for grid
grid2 <- readOGR(paste(getwd(), "/grid_5km_with_cost_dem_popn_aspect_globcover_PNG.shp", sep = ""))
# convert from m to km
grid2 <- spTransform(grid2, "+proj=utm +zone=54 +south +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=km +no_defs")

library(spatialEco)
# adjust values of grids with no deforestation
grid2@data$$def_mean[is.na(grid2@data$$def_mean)] <- 0
grid2@data$$def_count[is.na(grid2@data$$def_count)] <- 0
# eliminate grids with NA values for covariates
grid3 <- sp.na.omit(grid2, col.name = "cost_mean")
grid4 <- sp.na.omit(grid3, col.name = "elev_range")
grid5 <- sp.na.omit(grid4, col.name = "globcov_ma")

# 'deforested' grids defined as those with >10% of the maximum number of 'deforestation' pixels found in any grid
maxval <- max(grid5@data$$def_count)/10 # status: 1 = event (deforestation), 0 = censored ('still intact') grid5@data$$status <- ifelse(grid5@data$$def_count >= maxval, 1, 0) # set response variable (mean year of deforestation activity, censoring status) surv <- inla.surv(grid5@data$$def_mean, grid5@data$$status) # set up unique ID for each grid grid5@data$$id2 <- 1:nrow(grid5@data)
# use grid IDs and coordinates to set up INLA lattice graph
lattice_temp <- poly2nb(grid5, row.names = grid5@data$$id2) nb2INLA(paste(getwd(), "/lattice.graph", sep = ""), lattice_temp) Lattice.adj <- paste(getwd(), "/lattice.graph", sep = "") # set up intercept intercept1 = rep(1, nrow(grid5@data)) # INLA requires a dataframe as input grid10 <- grid5@data # z-scale the covariates grid10[c("elev_mean", "popn_mean", "cost_mean")] <- scale(grid10[c("elev_mean", "popn_mean", "cost_mean")]) # globcov_ma is a land cover category variable (e.g. forest, mangrove, cropland) - must be a factor grid10$$globcov_ma <- as.factor(grid10$globcov_ma) # formula and model coxinlaZ <- inla.coxph(surv ~ -1 + intercept1 + cost_mean + elev_mean + popn_mean + globcov_ma + f(id2, model = "bym", graph = Lattice.adj, scale.model = TRUE), list(surv = surv, cost_mean=grid10$$cost_mean, elev_mean=grid10$$elev_mean, popn_mean=grid10$$popn_mean, popn_mean=grid10$$popn_mean, globcov_ma=grid10$$globcov_ma, id2=grid10$$id2, intercept1 = intercept1)) rzPNGcov3 <- inla(coxinlaZ$$formula, family = coxinlaZ$$family, data=c(as.list(coxinlaZ$$data), coxinlaZ$$data.list), E = coxinlaZ$E,
control.compute = list(mlik = T))