我正在尝试在使用 "spatialEco" 包时整合 oring 统计曲线下的面积

I'm trying to integrate the area under the curve of the oring statistic in using "spatialEco" package

我需要在 Rstudio 中为 O 形圈统计数据整合曲线下面积。但是,spatialEco 包不会报告 O 形环统计的实际值,正如您在 spatstat 包的 Ripley's K 函数中看到的那样。这是到达我所在位置的代码。

library('spatstat')
library('ggplot2')
library('spatialEco')

set.seed(seed=24)

radiusCluster<-100
lambdaParent<-.02
lambdaDaughter<-30
hosts<-1000
randmod<-0
dim<-2000


numbparents<-rpois(1,lambdaParent*dim)

xxParent<-runif(numbparents,0+radiusCluster,dim-radiusCluster)
yyParent<-runif(numbparents,0+radiusCluster,dim-radiusCluster)

numbdaughter<-rpois(numbparents,(lambdaDaughter))
sumdaughter<-sum(numbdaughter)



thetaLandscape<-2*pi*runif(sumdaughter)

rho<-radiusCluster*sqrt(runif(sumdaughter))


xx0=rho*cos(thetaLandscape)
yy0=rho*sin(thetaLandscape)


xx<-rep(xxParent,numbdaughter)
yy<-rep(yyParent,numbdaughter)

xx<-xx+xx0

yy<-yy+yy0
cds<-data.frame(xx,yy)
is_outlier<-function(x){
  x > dim| x < 0
}
cds<-cds[!(is_outlier(cds$xx)|is_outlier(cds$yy)),]
while (nrow(cds)<hosts){
 dif<-hosts-nrow(cds)
  extraparentxx<-sample(xxParent,dif,replace = TRUE)
  extraparentyy<-sample(yyParent,dif,replace = TRUE)
  extrathetaLandscape<-2*pi*runif(dif)
  extrarho<-radiusCluster*sqrt(runif(dif))
  newextracoodsxx<-extrarho*cos(extrathetaLandscape)
  newextracoodsyy<-extrarho*sin(extrathetaLandscape)
  extraxx<-extraparentxx+newextracoodsxx
  extrayy<-extraparentyy+newextracoodsyy
  cdsextra<-data.frame(xx=extraxx,yy=extrayy)
  cds<-rbind(cds,cdsextra)
}


sampleselect<-sample(1:nrow(cds),hosts,replace=F)
cds<-cds%>%slice(sampleselect)

randfunction<-function(x){
  x<-runif(length(x),0,dim)
}
randselect<-sample(1:nrow(cds),floor(hosts*randmod),replace=F)
cds[randselect,]<-apply(cds[randselect,],1,randfunction)

landscape<-ppp(x=cds$xx,y=cds$yy,window=owin(xrange=c(0,dim),yrange=c(0,dim)))
ggplot(data.frame(landscape))+geom_point(aes(x=x,y=y))+coord_equal()+theme_minimal()

ostat<-o.ring(landscape,inhomogenous=FALSE)

这会生成 O 形环图:

是否可以整合此图来估计此曲线下的面积?

如果您查看 o.ring 函数内部,您会看到:


## > o.ring

function (x, inhomogeneous = FALSE, ...) 
{
    if (inhomogeneous) {
        g <- spatstat.core::pcfinhom(x, ...)
    }
    else {
        g <- spatstat.core::pcf(x, ...)
    }
    lambda <- summary(x)$intensity
    O <- spatstat.core::eval.fv(lambda * g)
    graphics::plot(O, ylab = "O-ring(r)", main = "O-ring")
}

现在自己执行这些步骤:


x <- landscape
g <- spatstat.core::pcf(x)
lambda <- summary(x)$intensity
O <- spatstat.core::eval.fv(lambda * g)

print.data.frame(head(O))

输出:


> print.data.frame( head(O) )
          r    theo       trans         iso
1 0.0000000 0.00025         Inf         Inf
2 0.9765625 0.00025 0.002681005 0.002671475
3 1.9531250 0.00025 0.001666361 0.001659885
4 2.9296875 0.00025 0.001351382 0.001345655
5 3.9062500 0.00025 0.001212683 0.001207095
6 4.8828125 0.00025 0.001141121 0.001135419

我想,尽管我需要查看 plot.fv 才能确定,trans 和 iso 是实线还是虚线。您知道要集成哪一个吗?

无论如何,这些都是你们各自的领域:


## filter out the Inf values
my.O <- O %>% filter( is.finite(trans) & is.finite(iso) )

library(zoo)

sum( diff(my.O$r) * rollmean( my.O$trans, 2 ) )

sum( diff(my.O$r) * rollmean( my.O$iso, 2 ) )