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:

    # 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")
    # 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,
id2=grid10$id2, intercept1 = intercept1))
    rzPNGcov3 <- inla(coxinlaZ$formula,
             family = coxinlaZ$family,
                 E = coxinlaZ$E,
                 control.compute = list(mlik = T))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.