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))