0

我正在尝试使用 R 中 spatstat 包中的 ppm() 函数对带有图像协变量的点过程进行建模。我将栅格转换为 im 对象以与 spatstat 一起使用，但使用 im 作为协变量时遇到了问题在模型中。像素值是数字的，但这些实际上只是不同景观区域的代码，所以问题的关键是让模型将像素值作为因子而不是数字来读取。我尝试了以下两种方法（R代码和数据如下所示）。第一个包括在将栅格对象转换为 im 对象之前将栅格值从数值转换为因子。使用 as.factor() 函数，这似乎具有将值转换为因子的预期效果。但是，当我使用这个协变量运行 ppm 模型时，ppm() 函数不包括模型中每个因子水平的参数（与参考水平相比）。相反，它将协变量视为数字，只有一个协变量的一个参数。第二种方法是运行 ppm 模型，其中因子（协变量）用于在公式参数中指定协变量，而不仅仅是协变量本身。这实际上适用于拟合模型，与参考相比，为每个因子水平提供了一个参数。但是，当我运行 predict.ppm() 来获得我的预测时，它会失败，因为我在公式参数中使用了 factor()。问题是，我如何运行 ppm 模型，使其将协变量图像的值识别为因子，从而为每个因子水平（减去参考值）拟合具有参数估计的模型，并允许使用 predict.ppm 进行预测()。

此处的点过程数据为 csv 格式：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

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

1 回答 1

1

问题是 R 没有因子值矩阵，因此将因子放入 中总是有点古怪im，但是一旦完成，一切都会正常工作。我的解决方案是将整数值光栅转换为im格式并从那里处理所有内容（我不是光栅包的常规用户）。

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

下面的语法ppm要求您运行 spatstat 1.37-0 或更高版本。此外，我正在使用Window1.38-0 中的新通用函数，您需要对旧版本稍作不同的处理（我强烈推荐最新版本 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)
于 2014-10-02T19:44:47.760 回答