R包spatstat:当像素图像值是数字时，如何使用点过程模型协变量作为因子 [英] R package spatstat: How to use point process model covariate as factor when pixel image values are numeric

问题描述

我正在尝试使用R中spatstat包中的ppm()函数使用图像协变量对点过程进行建模.我将栅格转换为im对象 与spatstat一起使用时，在使用im作为模型的协变量时遇到了问题.像素值是数字，但实际上只是代码 对于不同的景观区域，因此问题的症结在于使模型读取像素值作为因子而不是数值.我尝试了以下 两种方法(R代码和数据在下面介绍).第一种方法是先将栅格值从数字转换为因数，然后再转换 将栅格对象更改为im对象.使用as.factor()函数似乎具有将值转换为factor的预期效果.但是，当我 使用此协变量运行ppm模型，ppm()函数不包含模型中每个因子水平的参数(与参考水平相比).相当 它将协变量视为数值，仅将一个协变量的一个参数作为参数.第二种方法是使用系数(协变量)运行ppm模型 用于在公式参数中指定协变量，而不仅仅是协变量本身.这实际上可以拟合模型，给我一个参数 与参考相比，每个因子水平.但是，当我运行predict.ppm()进行预测时，它失败了，因为我在公式中使用了factor() 争论.问题是，我该如何运行ppm模型，使其将协变量图像的值识别为因子，从而使模型适合 每个因子水平的参数估计值(减去参考值)，并允许使用predict.ppm()进行预测.

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 run 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().

此处的点流程数据为csv格式: https://www.dropbox.com/s/tp1opzsmc14e2hb/EbolaData_AnalyticSet_8.8.14.csv?dl=0

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

协变量的tiff文件在此处: https://www .dropbox.com/s/0fyt0jflokrpp5z/anthrome2000_global_5min.tif?dl = 0

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

R代码如下:

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

推荐答案

问题是R没有因子值矩阵，因此将因子放入im总是有点古怪，但是一旦完成一切正常.我的解决方案是将整数值栅格转换为im格式并在那里处理所有内容(我不是栅格数据包的常规用户).

The problem is that R doesn't have factor valued matrices, so it always a bit quirky to get factors into an im, but once it is done everything works as it should. My solution was just to convert the integer valued raster into the im format and handle everything from there (I'm not a regular user of the raster package).

我必须加载maptools库才能使命令SP <- as(outbreaks, "SpatialPoints")正常工作.另外，由于第一列中有一些奇怪的字符(我们也不会使用)，R无法读取给定的csv文件，因此我必须删除这些字符才能正常工作.

I had to load the maptools library for the command SP <- as(outbreaks, "SpatialPoints") to work. Also, R couldn't read the given csv file due to some strange characters in the first column (which we don't use anyway), so I had to remove these for everything to work.

ppm的以下语法要求您正在运行spatstat 1.37-0或更高版本.此外，我正在使用1.38-0中的新通用函数Window，对于较旧的版本，您需要做些不同的操作(我强烈建议您使用最新版本1.38-1):

The syntax below for ppm requires that you are running spatstat 1.37-0 or newer. Furthermore, I'm using a new generic function Window from 1.38-0, and you need to do this slightly differently for older versions (I highly recommend the newest version 1.38-1):

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

# 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 we define our point process (with the bounding box as temporary window)
# for spatstat:
SP <- as(outbreaks, "SpatialPoints")
outbreaks.ppp <- as(SP, "ppp")

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

# Now let's convert that raster to an im object for use in spatstat
# (and make it into a factor):
anthro <- asImRaster(anthro)
anthro <- eval.im(as.factor(anthro))
levels(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")

# Make Africa into the observation window (of type mask):
Window(outbreaks.ppp) <- Window(anthro)

# See the data we have read in:
plot(anthro)
plot(outbreaks.ppp, add = TRUE)

# Fit model and predict:
spatial.m.1 <- ppm(outbreaks.ppp ~ anthro)
p <- predict.ppm(spatial.m.1)

这篇关于R包spatstat:当像素图像值是数字时，如何使用点过程模型协变量作为因子的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！