[GIS] R package spatstat: How to use point process model covariate as factor when pixel image values are numeric

rspatial statisticsspatstat

I am trying to model a point process with an image covariate using the ppm() function in the spatstat package in R. I convert my raster to an im object
for use with spatstat, and I run into a problem using the im as a covariate in the model. The pixel values are numeric, but these are actually just codes
for different landscape zones so the crux of the problem is getting the model to read the pixel values as factor rather than numeric. I have tried the following
two approaches (R code and data are presented below). The first consists of converting the raster values from numeric to factor prior to converting the
raster object to the im object. Using the as.factor() function this seems to have the desired effect of converting the values to factor. However, when I
run the ppm model with this covariate, the ppm() function does not include a parameter for each factor level in the model (compared to a reference level). Rather
it treats the covariate as numeric with just the one parameter for the one covariate. The second approach was to runt the ppm model with factor(covariate)
used to specify the covariate in the formula argument, rather than just the covariate itself. This actually works in fitting the model, giving me a parameter
for each factor level compared to the reference. However, when I run predict.ppm() to get my predictions it fails because I used factor() in the formula
argument. The qustion is, how can I run the ppm model such that it recognizes the values of the covariate image as factor and, thus, fitting a model with a
parameter estimate for each factor level (minus the reference) and allowing prediction with predict.ppm.

The point process data is in csv format here:
https://www.dropbox.com/s/tp1opzsmc14e2hb/EbolaData_AnalyticSet_8.8.14.csv?dl=0

The tiff file for the covariate is here:
https://www.dropbox.com/s/0fyt0jflokrpp5z/anthrome2000_global_5min.tif?dl=0

And the R code is as follows:

library(raster)
library(spatstat)
library(geostatsp)

# First set the geographic extent we'll be using
e <- extent(-20, 60, -40, 35)

# Then read in the point process data in .csv file:
outbreaks <- read.csv("EbolaData_AnalyticSet_8.8.14.csv")
coordinates(outbreaks) <- ~Longitude+Latitude
proj4string(outbreaks) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") 

# Then read in the anthropogenic biome tiff and crop it
anthro2000 <- raster("anthrome2000_global_5min.tif")
anthro2000.africa <- crop(anthro2000, e)

# Then we define oour point process and window for spatstat:
SP <- as(outbreaks, "SpatialPoints")
outbreaks.ppp <- as(SP, "ppp")

# Now let's create a window
SP.win <- as(e, "SpatialPolygons")
W <- as(SP.win, "owin")

# Before creating the im object, let's convert the pixel values in raster to factor:
is.factor(anthro2000.africa)
f <- as.factor(anthro2000.africa)
is.factor(f)
rat <- levels(f)[[1]]
rat$anthro <- c("Urban", "Mixed Settle", "Rice Villages", "Irrigated villages", "Rainfed villages", "Pastoral vilalges",
    "Resid. irrig. cropland", "Resid. rainfed cropland", "Pop. cropland", "Remote cropland", 
    "Resid. rangeland", "Pop. rangeland", "Remote rangeland", "Resid. forests", "Pop. forests",
    "Remote forests", "Inhabited treeless and barren", "Wild forests", "Wild treeless and Barren")
rat$code <- c(11,12,21,22,23,24,31,32,33,34,41,42,43,51,52,53,54,61,62)
levels(f) <- rat

# Now let's convert that raster to an im object for use in spatstat:
anthro2000.africa.img <- asImRaster(f)

# Everything is now set up for the ppm models

# Aprroach 1
spatial.m.1 <- ppm(outbreaks.1.ppp, ~ Cov1, covariates = list(Cov1 = anthro2000.africa.img))
spatial.m.1 # Notice the model is fitted, however the pixel values of the covariate are not interepreted as factor


# Approach 2:
spatial.m.2 <- ppm(outbreaks.ppp, ~ factor(Cov1), covariates = list(Cov1 = anthro2000.africa.img)) # Noticce the use of factor() here to force the covariate as factor
spatial.m.2 # Here the model is fitted correctly with covariate as factor and thus each factor level has a parameter estimate in the model (relative to the ref)

# However, the approach does not allow me to run the predictions:
p <- predict.ppm(spatial.m.2, covariates = list(Cov1 = anthro2000.africa.img))

Best Answer

The problem here is a bug in asImRaster, a function from the package geostatsp designed to convert between image formats. After the command anthro2000.africa.img <- asImRaster(f), if you had printed the object anthro2000.africa.img, the printout would have revealed that the pixel values are numeric. The factor levels have been discarded by asImRaster.

You can work around this bug by reinstating the levels:

levels(anthro2000.africa.img) <- levels(f).

[Always print the resulting object to check it is what you expected.]

Then you can proceed to fit the model using ppm as you did before.

[Note: contrary to one of the comments above, this is a perfectly valid point process model. It is explained in detail in chapter 9 of the spatstat book.]

Related Solutions

[GIS] Converting raster to im object for point process model covariate in R

If you would check help for as.im() you would notice that raster data type is not supported. We should thank @mdsumner for bringing up my old question (I do not remember the occasion that forced me to ask it - suppose it was just some testing and I wouldn't recall that question by myself). I used there a geostatsp package that provides needed functionality.

In your case it is just:

library(geostatsp)
mgvf.img <- as.im(mgvf.2001.africa)

Notice that if you don't want to use geostatsp (for no good reason), as.im() will work with matrices. So you will need to as.matrix(mgvf.2001.africa) and then rotate it (several times probably - I'm too sleepy to test it) with

rotate = function(mat) t(mat[nrow(mat):1,,drop=FALSE]) # mat - matrix

and then pass rotated matrix to as.im().

[GIS] using raster.predict for Cubist model, factor levels list

I stick with what I said in my previous comment on this post. I believe that the issue is that you are forcing factors rather than letting the raster predict function handle them. The raster predict wrapper code has explicit factor handling build in.

You should coerce the appropriate variables in your training data to factors and not your rasters. To factorize a raster it has to be pulled into memory, which is not necessary and defeats a memory safe approach using raster. The model object stores class types so, the raster predict function knows if a given variable, corresponding to a raster in the stack, is a factor or numeric. Factorial levels that are present in the new data (raster stack) but not in the training data are assigned NA (NoData).

Following your coding logic it would be more appropriate to use an sp class SpatialPixelsDataFrame or SpatiaGridDataFrame object where, each column represents a covariate and factor coercion is appropriate. In this case, you would just predict a new column to the object.