解释具有不同标记多类型点模式的点之间的空间关系 - 使用 spatstat - 并更好地可视化结果

编程问答 2022-05-02

问题描述

首先让我提到我没有接受过统计学方面的培训，而且我是这个领域的新手，无论是空间统计还是 R 编码。 但是，我正在尝试学习如何识别和表征点绘制模式的空间分布以及不同类型的点（标记）之间， 因为我的论文需要这样做。但是我现在已经到了需要帮助来了解我是否在正确的轨道上，或者更确切地说我是否使用了正确的命令的地步。 此外，我还有一些问题列在下面。

潜在的情况是我想破译考古发掘现场的内部模式。挖掘场地被划分为1平方米的网格。 每个点代表一个带有分类（因子值）标记的人工制品。我的主要目标是确定具有不同属性的工件之间的空间关系。 为此，我首先绘制了一个点模式图并进行了二次检验以识别空间特征（试图拒绝零假设，因此拒绝 CSR）， 然后执行 Ripley's K 函数来识别点过程是聚合还是分离。 样方检验的 p 值 Download the dataset

# ----------------------------------------------------------------------
# Prepare data
# ----------------------------------------------------------------------
    
# Read in data file
roma <- read.csv(file = "roma_points.csv",header = T) 

# Define the window (class owin) for analyses with three holes
W <- owin(poly=list(list(x=c(100,104,102,100),y=c(100,100,106,110,110)),list(x=c(100,101,101),y=c(102,103,102)),y=c(106,107,106)),y=c(107,108,107))))


# Define ppp object based on this project's point locations
roma.ppp <- ppp(x = roma$EAST,y = roma$NORTH,window = W,marks = factor(roma$MARK))

roma.ppp <- as.ppp(roma.ppp)
unitname(roma.ppp) <- c("metre","metres")

# Plotting point pattern map
plot(roma.ppp,cex = 0.5,main = "Point pattern map",axes = F,ann = F,cols = c("blue","black","red","green","sienna"))

# plot the point distribution of each mark
plot(split(roma.ppp),main = "Point pattern map of individual marks")

# ----------------------------------------------------------------------
# Chi-squared test with inhomogeneity taken into account
# ----------------------------------------------------------------------

roma.squared.inhom <- quadrat.test(roma.ppp,nx = 4,ny = 10,method = "MonteCarlo",lambda = density(roma.ppp))

    # plot density map and display grid info
plot(density(roma.ppp),main = "Density plot + chi-squared results for Inhomogeneous Density")
plot(roma.squared.inhom,col = "white",add = T)

Point pattern maps of individual marks

我的点模式总结：

Marked planar point pattern:  2351 points
Average intensity 81.06897 points per square metre

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 metres

Multitype:
  frequency proportion  intensity
A       671 0.28541050 23.1379300
B       454 0.19310930 15.6551700
C        27 0.01148447  0.9310345
D       191 0.08124202  6.5862070
E      1008 0.42875370 34.7586200

Window: polygonal boundary
single connected closed polygon with 14 vertices
enclosing rectangle: [100,104] x [100,110] metres
                     (4 x 10 metres)
Window area = 29 square metres
Unit of length: 1 metre
Fraction of frame area: 0.725

现在，我最终想要实现的是能够对不同的人工制品类型如何相互分配做出陈述/假设，更具体地说是例如 人工制品类型 A 聚集到人工制品类型 D？ 我（模糊！）的理解是，我需要对多类型模式的交叉类型 L 函数的非齐次版本进行估计：lcross.inhom。

现在，我的问题是：

我。正确的做法？
假设我最适合使用 lcross.inhom 是否正确，特别是考虑到强度，因为我的点模式不均匀？ 我读过有人建议为此执行一个可爱的二元 K 函数（localKcross 命令）。但我不明白结果，坦率地说，我并没有真正理解其中的区别。

二。如何正确执行Lcross.inhom？
关于计算 Lcross.inhom 函数。我认为我必须调整不同标记点的强度或拟合泊松模型以在存在空间不均匀性的情况下测试聚类（即零假设应该是不均匀的泊松过程）。 我尝试了两种不同的方法：

# ----------------------------------------------------------------------
# Ripley's K Function with inhomogeneity taken into account
# ----------------------------------------------------------------------

# Approach No1 ---------------------------------------------------------

# Calculating bandwidth sigma for the kernel estimator of point process
# intensity with cross validated bandwidth selection (bw.scott)
sigmaA <- with(split(roma.ppp),bw.scott(A))
sigmaB <- with(split(roma.ppp),bw.scott(B))
sigmaC <- with(split(roma.ppp),bw.scott(C))
sigmaD <- with(split(roma.ppp),bw.scott(D))

# Calculating with the Kernel smoothed intensity of point pattern
# (density.ppp) the values of the estimated intensity of the sub-
# process of points of type i and type j
lambdaA <- with(split(roma.ppp),density(A,sigmaA))
lambdaB <- with(split(roma.ppp),density(B,sigmaB))
lambdaC <- with(split(roma.ppp),density(C,sigmaC))
lambdaD <- with(split(roma.ppp),density(D,sigmaD))

# Calculating and plotting Lross.inhom with estimated intensities
# by non-parametric smoothing
par(mfrow = c(1,3))
plot(Lcross.inhom(roma.ppp,"A","B",lambdaA,lambdaB),. - r ~ r,main = "Approach 1 (A against B)",legend = FALSE)
plot(Lcross.inhom(roma.ppp,"C",lambdaC),main = "Approach 1 (A against C)","D",lambdaD),main = "Approach 1 (A against D)",legend = FALSE)



# Approach No2---------------------------------------------------------

# Fit parametric intensity model,by first fitting a 
# Poisson point process model to the observed data
fit <- ppm(roma.ppp ~marks * polynom(x,y,3))

# evaluate fitted intensities at data points
# (these are the intensities of the sub-processes of each type)
inten <- fitted(fit,dataonly=TRUE)

# split according to types of points
lambda <- split(inten,marks(roma.ppp))

# Calculate and plot Lcross.inhom for different marks
par(mfrow = c(1,lambda$A,lambda$B),main = "Approach 2 (A against B)",lambda$C),main = "Approach 2 (A against C)",lambda$D),main = "Approach 2 (A against D)",legend = FALSE)

Results of both approaches calculating Kcross.inhom

但是，这两种方法中的每一种的结果都非常不同。第一种方法表示 A 型和 B 型之间的聚类，但仅达到大约 0.3 m，而第二种方法表示这两种类型之间的聚类超过 1 米。在类型 A 和类型 D 之间非常相似。另一方面，在类型 A 和 C 之间，方法 1 表示再次聚类，直到大约 r=0.3 m，而方法 2 仅表示聚类到 r=0.8 m。 现在，我假设我可能没有足够的分数 C，因此得到有问题的结果。 但是对于其他类型，我不确定选择哪种方法。如果我将这些图与每个标记的点模式图进行比较，第一种方法似乎更合适。 不过，我并不完全确定。我是否在为我的点模式计算 Lcross.inhom 函数时犯了一个错误，或者我是否误解了如何正确解释这些图？因为结果只绘制到 r=1 米，所以通常很难在这里得出结论吗？我不知道如何改变它，我读过我可以用 nncross 为每个标记确定 r？

三。如何在点模式图上可视化结果？
有没有其他方法可以用点绘制的地图更好地显示 lcross.inhom 的结果？ 我想可视化可以观察到聚类和分散模式的位置，重要性级别和规模。 假设我想用我正在测试的两种类型绘制点模式图，然后突出显示它们彼此聚集的区域。 我已经读到我可以通过计算每个观察点的 K 函数来做到这一点，然后绘制它自己在给定距离 d 处拒绝零假设的显着性水平。 不幸的是，我完全不知道该怎么做。

我希望这不会太令人困惑，我要求的也不会太多。我已经搜索并阅读了很多，分别是 spatstat 书中的相应章节，但我必须承认，作为非统计学家，我发现大多数参考文献都相当繁重。 在这一点上，我只是只见树木不见森林。我希望这一切都有意义，并且您能够理解我到目前为止所做的/尝试的内容。我很抱歉阅读太长，但我尽量做到准确。我也很高兴收到有关如何改进我的第一篇帖子的提示或建议。

我的感谢和赞赏。祝一切顺利，保持安全
斯科蒂

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

distributionp-valuerrspatialspatstat