Search code examples
rspatial

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

I need to integrate the area under the curve for the O-ring statistic in Rstudio. However, the package spatialEco does not report the actual values of the O-ring statistic as you would see in the Ripley's K function from the package spatstat. Here is the code to get to the point where I am at.

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)

This produces the O-ring plot: enter image description here

Is it possible to integrate this plot to estimate the area under this curve?

Solution

  • If you peek inside the o.ring function, you see this:

    
## > 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")
}

    Now execute these steps yourself:

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

print.data.frame(head(O))

    Outputs:

    
> 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

    I guess, though I 'd need to look at plot.fv to be sure, that trans and iso is either of the solid or dashed line. Do you know which one you want to integrate?

    In any case, these are your respective areas:

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