使用 spatstat 中的 "envelope" 函数对空间点模式超帧进行基于仿真的假设检验
Simulation-based hypothesis testing on spatial point pattern hyperframes using "envelope" function in spatstat
我想使用复制的空间点模式在 spatstat 中执行假设检验。 spatstat 有精彩的文档，您可以在此处找到有关复制点模式分析的详细信息：https://cran.r-project.org/web/packages/spatstat/vignettes/replicated.pdf
几个映射的森林图将代表在不同位置工作的点过程。每个点图案都有标记，其中每个标记代表一个树种。此外，每个点模式都与一个协变量光栅图像相关联。首先，我想创建一个空模型，假设每个标记与协变量都有不同的关系。然后我想使用这个空模型通过模拟随机标记假设并绘制生成的模拟包络线来测试某些物种是否相互关联（或避免）。 （随机标记假说指出标记是随机分配给点模式中的点的。）
首先，我将展示我通常如何使用单点模式进行此分析。然后我将解释我在使用超帧执行相同分析时遇到的两个问题。
假设您有一个标记点图案，其中每个标记都是树种。我将使用 spatstat 中可用的 "lansing" 数据集：
library(spatstat)
data(lansing)
par(mar=rep(0.5,4))
plot(split(lansing),main="")
现在假设您想查看一些空间协变量（例如，土壤养分或水分），因此您创建了测量的核平滑密度的栅格图像：
sim1 <- rpoispp(function(x,y) {500 * exp(-3*x)}, win=owin(c(0,1),c(0,1)))
sim1 <- density(sim1)
首先创建空模型：
single.mod <- ppm(lansing ~ marks*sim1)
其中 "ppm" 函数将 "marks" 识别为具有物种名称的列，并将 "sim1" 识别为协变量。
然后执行基于模拟的假设检验，您感兴趣的是黑橡树和枫树是否出现在同一位置。
single.E<-envelope.ppm(single.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(lansing)))
par(mar=rep(4,4))
plot(E,legend=FALSE,ylab="L-function",xlab="Spatial scale (m)",main="Testing random label hypothesis \nfor single point pattern")
这很好用。现在，如果我们出去采样更多的地块以使我们的分析更加稳健，我们可以将额外的地块合并到超框架中，其中每个地块都有自己的点模式，并被视为 "experimental" 重复。每个地块也将获得自己的空间协变量：
sim2 <- rpoispp(function(x,y) {500 * exp(-2*y)}, win=owin(c(0,1),c(0,1)))
sim2 <- density(sim2)
sim3 <- rpoispp(500, win=owin(c(0,1),c(0,1)))
sim3 <- density(sim3)
hyper <- hyperframe(pp=list(lansing,lansing,lansing),sims=list(sim1,sim2,sim3))
Sim2 和 sim3 是我们收集的另外两个图的空间协变量，"hyperframe" 函数将三点模式及其关联的空间协变量组合到一个超帧中。
我想使用 "mppm" 构建模型（用于为多点模式创建模型），其中每个点模式都由它们的空间协变量解释，"sims":
hyper.mod <- mppm(pp ~ sims, data = hyper)
当我试图让每个标记与协变量具有不同的关系时，出现了第一个问题：
int.mod <- mppm(pp ~ marks*sims, data=hyper)
出现以下错误消息：
Error in checkvars(formula, data.sumry$col.names, extra = c("x", "y", : Variable "marks" in formula is not one of the names in data"
我得到同样的错误使用：
int.mod <- mppm(pp ~ pp$marks*sims, data=hyper)
第二个问题是弄清楚如何运行 在超帧上进行基于模拟的假设检验。让我们使用有效的超帧模型 (hyper.mod) 来试试这个：
E.hyper <- envelope(hyper.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(pp)))
您收到一条错误消息：
Error in UseMethod("envelope") : no applicable method for 'envelope' applied to an object of class "c('mppm', 'list')"
暗示 "envelope" 不适用于 mppm 对象（仅 ppp 或 ppm）。我怀疑有一种方法可以巧妙地绕过这个限制，但我还没有找到它。任何建议或指导都会非常有帮助！
第一个错误是spatstat开发版几天前修复的bug。如果您有 devtools 软件包，您可以轻松获得最新（也是最好的）spatstat：
devtools::install_github('spatstat/spatstat')
如果这不是您的选择（以及您使用的平台），请告诉我。
第二个错误确实是因为 envelope 没有为 class mppm 实现，所以你现在必须像你说的那样设计一个解决方法。
我认为你到目前为止所做的有几个问题：模型是不均匀的，但你使用 Lcross 而不是 Lcross.inhom 我认为你的 single.E相当于（你根本不使用拟合模型的地方）：
single.E <- envelope(lansing, Lcross, i="blackoak", j="maple", nsim=39, rank=1, global=TRUE, correction="best", simulate=expression(rlabel(lansing)))
请告诉我您的进展情况。我可能会找时间提供更多有关缺失 envelope.mppm 解决方法的详细信息（每个模式的汇集汇总函数）。
spatstat中没有envelope.mppm，因为一些统计问题（与多重假设检验相关）尚未解决。
最快的解决方案可能是使用 cdf.test.mppm，它将执行测试并给出一些图形输出。
从不同的点模式中汇集包络（例如 K 函数）以获得单个包络是合理的提供 拟合模型意味着不同的模式在统计上应该是等效的。例如，如果模型包含对于不同模式不同的协变量，那将是无效的。
更好的策略可能是为每个模式绘制全局包络线，并使用产品规则进行多次测试。假设数据中有 M 个点模式，并且您想要检验（拟合模型的）显着性水平 alpha（通常 alpha = 0.05）。然后，您将构建 M 个信封，每个模式一个，每个信封的显着性水平 gamma = 1 - (1-alpha)^(1/M)。每个信封将由 envelope.ppp 生成，大约 global=TRUE 和 nsim = 1/gamma - 1。示例：如果 M = 10 且 alpha = 0.05，则 gamma=1 - 0.95^(1/10) = 0.0051 所以 nsim=1/gamma-1 = 194.45，称其为 nsim=195。由于您要制作全局包络，您可能需要两倍的模拟次数，如 envelope 的帮助中所解释的那样。所以如下，其中fit是拟合模型，H是数据的原始超帧：
sims <- simulate(fit, nsim=2*195)
SIMS <- list()
for(i in 1:nrow(sims)) SIMS[[i]] <- as.solist(sims[i,,drop=TRUE])
Hplus <- cbind(H, hyperframe(Sims=SIMS))
这通过添加一列模拟模式来扩充原始超帧（列中的每个条目都是一个包含 2*195 模式的 'solist'）。然后执行（其中 X 是包含原始点模式数据集的 H 的列）
EE <- with(Hplus, envelope(X, Lest, global=TRUE, nsim=195, simulate=Sims))
plot(EE)
这会生成一个包含许多包络面板的绘图。他们的解释是，如果任何一个观察到的 L 函数在相应的包络线之外徘徊，则结果很重要 - 检验拒绝零假设，即拟合模型为真。
我想使用复制的空间点模式在 spatstat 中执行假设检验。 spatstat 有精彩的文档，您可以在此处找到有关复制点模式分析的详细信息：https://cran.r-project.org/web/packages/spatstat/vignettes/replicated.pdf
几个映射的森林图将代表在不同位置工作的点过程。每个点图案都有标记，其中每个标记代表一个树种。此外，每个点模式都与一个协变量光栅图像相关联。首先，我想创建一个空模型，假设每个标记与协变量都有不同的关系。然后我想使用这个空模型通过模拟随机标记假设并绘制生成的模拟包络线来测试某些物种是否相互关联（或避免）。 （随机标记假说指出标记是随机分配给点模式中的点的。）
首先，我将展示我通常如何使用单点模式进行此分析。然后我将解释我在使用超帧执行相同分析时遇到的两个问题。
假设您有一个标记点图案，其中每个标记都是树种。我将使用 spatstat 中可用的 "lansing" 数据集：
library(spatstat)
data(lansing)
par(mar=rep(0.5,4))
plot(split(lansing),main="")
现在假设您想查看一些空间协变量（例如，土壤养分或水分），因此您创建了测量的核平滑密度的栅格图像：
sim1 <- rpoispp(function(x,y) {500 * exp(-3*x)}, win=owin(c(0,1),c(0,1)))
sim1 <- density(sim1)
首先创建空模型：
single.mod <- ppm(lansing ~ marks*sim1)
其中 "ppm" 函数将 "marks" 识别为具有物种名称的列，并将 "sim1" 识别为协变量。
然后执行基于模拟的假设检验，您感兴趣的是黑橡树和枫树是否出现在同一位置。
single.E<-envelope.ppm(single.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(lansing)))
par(mar=rep(4,4))
plot(E,legend=FALSE,ylab="L-function",xlab="Spatial scale (m)",main="Testing random label hypothesis \nfor single point pattern")
这很好用。现在，如果我们出去采样更多的地块以使我们的分析更加稳健，我们可以将额外的地块合并到超框架中，其中每个地块都有自己的点模式，并被视为 "experimental" 重复。每个地块也将获得自己的空间协变量：
sim2 <- rpoispp(function(x,y) {500 * exp(-2*y)}, win=owin(c(0,1),c(0,1)))
sim2 <- density(sim2)
sim3 <- rpoispp(500, win=owin(c(0,1),c(0,1)))
sim3 <- density(sim3)
hyper <- hyperframe(pp=list(lansing,lansing,lansing),sims=list(sim1,sim2,sim3))
Sim2 和 sim3 是我们收集的另外两个图的空间协变量，"hyperframe" 函数将三点模式及其关联的空间协变量组合到一个超帧中。
我想使用 "mppm" 构建模型（用于为多点模式创建模型），其中每个点模式都由它们的空间协变量解释，"sims":
hyper.mod <- mppm(pp ~ sims, data = hyper)
当我试图让每个标记与协变量具有不同的关系时，出现了第一个问题：
int.mod <- mppm(pp ~ marks*sims, data=hyper)
出现以下错误消息：
Error in checkvars(formula, data.sumry$col.names, extra = c("x", "y", : Variable "marks" in formula is not one of the names in data"
我得到同样的错误使用：
int.mod <- mppm(pp ~ pp$marks*sims, data=hyper)
第二个问题是弄清楚如何运行 在超帧上进行基于模拟的假设检验。让我们使用有效的超帧模型 (hyper.mod) 来试试这个：
E.hyper <- envelope(hyper.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(pp)))
您收到一条错误消息：
Error in UseMethod("envelope") : no applicable method for 'envelope' applied to an object of class "c('mppm', 'list')"
暗示 "envelope" 不适用于 mppm 对象（仅 ppp 或 ppm）。我怀疑有一种方法可以巧妙地绕过这个限制，但我还没有找到它。任何建议或指导都会非常有帮助！
第一个错误是spatstat开发版几天前修复的bug。如果您有 devtools 软件包，您可以轻松获得最新（也是最好的）spatstat：
devtools::install_github('spatstat/spatstat')
如果这不是您的选择（以及您使用的平台），请告诉我。
第二个错误确实是因为 envelope 没有为 class mppm 实现，所以你现在必须像你说的那样设计一个解决方法。
我认为你到目前为止所做的有几个问题：模型是不均匀的，但你使用 Lcross 而不是 Lcross.inhom 我认为你的 single.E相当于（你根本不使用拟合模型的地方）：
single.E <- envelope(lansing, Lcross, i="blackoak", j="maple", nsim=39, rank=1, global=TRUE, correction="best", simulate=expression(rlabel(lansing)))
请告诉我您的进展情况。我可能会找时间提供更多有关缺失 envelope.mppm 解决方法的详细信息（每个模式的汇集汇总函数）。
spatstat中没有envelope.mppm，因为一些统计问题（与多重假设检验相关）尚未解决。
最快的解决方案可能是使用 cdf.test.mppm，它将执行测试并给出一些图形输出。
从不同的点模式中汇集包络（例如 K 函数）以获得单个包络是合理的提供 拟合模型意味着不同的模式在统计上应该是等效的。例如，如果模型包含对于不同模式不同的协变量，那将是无效的。
更好的策略可能是为每个模式绘制全局包络线，并使用产品规则进行多次测试。假设数据中有 M 个点模式，并且您想要检验（拟合模型的）显着性水平 alpha（通常 alpha = 0.05）。然后，您将构建 M 个信封，每个模式一个，每个信封的显着性水平 gamma = 1 - (1-alpha)^(1/M)。每个信封将由 envelope.ppp 生成，大约 global=TRUE 和 nsim = 1/gamma - 1。示例：如果 M = 10 且 alpha = 0.05，则 gamma=1 - 0.95^(1/10) = 0.0051 所以 nsim=1/gamma-1 = 194.45，称其为 nsim=195。由于您要制作全局包络，您可能需要两倍的模拟次数，如 envelope 的帮助中所解释的那样。所以如下，其中fit是拟合模型，H是数据的原始超帧：
sims <- simulate(fit, nsim=2*195)
SIMS <- list()
for(i in 1:nrow(sims)) SIMS[[i]] <- as.solist(sims[i,,drop=TRUE])
Hplus <- cbind(H, hyperframe(Sims=SIMS))
这通过添加一列模拟模式来扩充原始超帧（列中的每个条目都是一个包含 2*195 模式的 'solist'）。然后执行（其中 X 是包含原始点模式数据集的 H 的列）
EE <- with(Hplus, envelope(X, Lest, global=TRUE, nsim=195, simulate=Sims))
plot(EE)
这会生成一个包含许多包络面板的绘图。他们的解释是，如果任何一个观察到的 L 函数在相应的包络线之外徘徊，则结果很重要 - 检验拒绝零假设，即拟合模型为真。