R/Ktest.r

In dbmss: Distance-Based Measures of Spatial Structures

Ktest <- function(X, r) {

  # Verify X
  if (!inherits(X, "ppp")) {
    stop("X is not of class ppp")
  }
  # Verify r
  if (!is.numeric(r) && !is.vector(r)) {
    stop("r must be a numeric vector")
  }
  # Eliminate 0 because there is no neighbor at distance 0
  r <- r[r != 0]
  if (length(r) < 2) {
    stop(paste("r has length", length(r), "- must be at least 2"))
  }
  if (any(diff(r) <= 0)) {
    stop("successive values of r must be increasing")
  }


  # Utilities
  # r : radius
  # w, l : width and length of the window
  #-------------------------------------------------------------------------------
  # ern : Probability for a point to be in another's neighborhood.
  # (= pi * r^2 / (w*l) corrected from edge effects)
  ern <- function(r, w, l) {
    return(r * r / w / l * (pi - 4 * r / 3 * (1 / w + 1 / l) + r * r / 2 / w / l))
  }
  #-------------------------------------------------------------------------------
  # espKc : Expectation of K, when density is known
  espKc <- function(r,w,l) {
    return(w * l * ern(r,w,l))
  }
  #-------------------------------------------------------------------------------
  # espKi : Expectation of K, when density is unknown.
  # rho*w*l is estimated by the number of points (N).
  # Difference between espKi and espKc goes to 0 when points are more than 20.
  espKi <- function(r, N, w, l) {
    return(espKc(r, w, l) * (1 - (1 + N) * exp(-N)))
  }
  #-------------------------------------------------------------------------------
  # eh : Expectation of h_1^2(U, r)
  eh <- function(r, w, l) {
    return(
      r^5 * (w + l) / 2 / w^3 / l^3 * (8 * pi / 3 - 256 / 45) +
      r^6 / w^3 / l^3 * (11 * pi / 48 + 8 / 9 - 16 / 9 * (w + l) * ( w + l) / w / l) +
      4 / 3 * r^7 * (w + l) / w^4 / l^4 - r^8 / 4 / w^4 / l^4
    )
  }
  # ------------------------------------------------------------------------------
  # sigmaKi : Variance matrix for the estimator of K (unknown rho),
  #           normalized by multiplication by par w*l*rho
  # vec     : Vector of distances to compute K
  #           (example : c(1, 2) to calculate K(1) and K(2)
  # N       : number of points
  # w       : width of the rectangle window
  # l       : length of the rectangle window
  sigmaKi <- function(vec, N, w, l) {
    # Intermediate computations
    d <- length(vec)
    ern_ <- ern(vec, w, l)

    c1 <- 2 * w * w * l * l / (N * (N - 1))
    c2 <- w * w * l * l * exp(-N) * (1 + N) * (1 - exp(-N) - N * exp(-N))
    c3 <- 2 * (N - 2) * c1

    # Preparation of a square matrix
    sigmaKi_ <- matrix(nrow = d, ncol = d)

    # Top half of the matrix
    for (i in seq_len(d - 1)) {
      for (j in seq_len(d)) {
        covh1_ <- covh1(r1 = vec[i], r2 = vec[j], w = w, l = l)
        sigmaKi_[i, j] <- c1 * ern_[min(i, j)] +
          (c2 - c1) * ern_[i] * ern_[j] + c3 * covh1_
      }
    }
    # Diagonal
    for (i in seq_len(d)) {
      sigmaKi_[i, i] <- c1 * ern_[i] + (c2 - c1) * ern_[i] * ern_[i] +
        c3 * eh(vec[i], w, l)
    }
    # Bottom half
    for (j in seq_len(d - 1)) {
      for (i in j:d) {
        sigmaKi_[i, j] <- sigmaKi_[j, i]
      }
    }
    # Result
    return(sigmaKi_)
  }
  #-------------------------------------------------------------------------------
  # invsqrtmat : Transforms a matrix into the square root of its inverse
  #             such as invsqrtmat %*% mat %*% t(invsqrtmat) = Id
  invsqrtmat <- function(mat) {
    if (length(mat) > 1) {
      e <- eigen(mat)
      # Eigen vectors
      p <- e$vectors
      # Square roots of eigen values
      d <- sqrt(e$values)
      # put into a diagonal matrix
      rd <- diag(d)
      # Resolution
      return(solve(p %*% rd))
    } else {
      return(1 / sqrt(mat))
    }
  }
  #-------------------------------------------------------------------------------
  # covh1  : covariance of h1
  # r1, r2 : distances
  # w      : width of the rectangle window
  # l      : length of the rectangle window
  covh1 <- function(r1, r2, w, l) {
    # Values of the product in the corner.
    # Coordinates are x'=(n-xi)/r2 without normalization in r'^2r^2/w^2/l^2
    # This function must be inside covh1 because it needs to use local variables,
    # that can not be passed as parameters because adaptIntegrate forbides it.
    corner <- function(x){
      (foncA1(ra2 * x / ra1, ra1, w, l) +
        foncA2(ra2 * x / ra1, ra1, w, l) +
        foncA3(ra2 * x / ra1, ra1, w, l) +
        foncA4(ra2 * x / ra1, ra1, w, l)) *
        (foncA3(x, ra2, w, l) + foncA4(x, ra2, w, l)
      )
    }

    # r values must be ordered
    ra1 <- min(r1, r2)
    ra2 <- max(r1, r2)

    # Size ratios parameters
    rapr <- ra1 / ra2
    r12 <- ra1 * ra1 / w / l
    r22 <- ra2 * ra2 / w / l

    # Biases, normalized by n^2 / r^2
    b1 <- brn(ra1, w, l)
    b2 <- brn(ra2, w, l)

    # Numerical computing of the elliptic integral
    int2 <- stats::integrate(
      integrand3,
      lower = 0,
      upper = 1,
      r1 = ra1,
      r2 = ra2
    )
    intcorner <- cubature::adaptIntegrate(
      corner,
      lowerLimit = c(0, 0),
      upperLimit = c(1, 1)
    )

    # line 1
    covh1_ <- (w - 2 * r2) / w * (l - 2 * r2) / l * b1 * b2
    # line 2
    covh1_ <- covh1_ + 2 * (w + l - 4 * r2) * r2 / w / l * b1 * (b2 - foncG(1))
    # line 3
    covh1_ <- covh1_ +
      2 * (w + l - 4 * r2) * r1 / w / l * (int2$value - b2 * foncG(1)) +
      4 * r22 * intcorner$integral

    # multiplication by the common factor
    covh1_ <- covh1_ * r12 * r22

    # Result
    return(covh1_)
  }

  foncg   <- function(x){
    if (x <= 1) {
      return(acos(x) - x * sqrt(1 - x * x))
    } else {
      return(0)
    }
  }
  foncg01 <- function(x){
    acos(x) - x * sqrt(1 - x * x)
  }
  foncG   <- function(x){
    x * acos(x) - sqrt(1 - x * x) * (2 + x * x) / 3 + 2 / 3
  }

  # Values of h1 on different zones without normalization by r^2/w/l
  brn    <- function(r, w, l){
    4 * r * (w + l) / (3 * w * l) - r * r / (2 * w * l)
  }
  indic  <- function(a){
    as.numeric(a)
  }
  foncA1 <- function(x, r, w, l){
    brn(r, w, l) * indic(x[1] >= 1) * indic(x[2] >= 1)
  }
  foncA2 <- function(x, r, w, l){
    (brn(r, w, l) - foncg(x[2])) * indic(x[1] >= 1) * indic(x[2] < 1) +
      (brn(r, w, l) - foncg(x[1])) * indic(x[2] >= 1) * indic(x[1] < 1)
  }
  foncA3 <- function(x, r, w, l){
    (brn(r, w, l) - foncg(x[1]) -
       foncg(x[2])) * indic(x[1] < 1) * indic(x[2] < 1) * indic(x[1]^2 + x[2]^2 > 1)
  }
  foncA4 <- function(x, r, w, l){
    (brn(r, w, l) + x[1] * x[2] -
      (foncg(x[1]) + foncg(x[2])) / 2 - pi / 4) * indic(x[1] < 1) *
      indic(x[2] < 1) * indic(x[1]^2 + x[2]^2 <= 1)
  }

  integrand3 <- function(x, r1, r2){
    foncg01(r1 * x / r2) * foncg01(x)
  }
  #-------------------------------------------------------------------------------
  # End of utilities

  # fails if the window is not a rectangle
  if (!is.rectangle(X$window)) stop("The window must be a rectangle to apply Ktest.")

  #  value = NA if the number of points is less than 2
  if (X$n > 1) {
    # Width - Length
    l <- X$window$xrange[2] - X$window$xrange[1]
    w <- X$window$yrange[2] - X$window$yrange[1]
    # espKi : expectation of K, rho unknown,
    # calculated according to the number of points.
    espKi_ <- espKi(r, X$n, w, l)

    # Computed variance. Depends on the observed number of points and the window.
    sigmaKi_ <- sigmaKi(r, X$n, w, l)

    # Distance matrix.
    # Requires a lot of RAM. Limited to 8,000 points with 32-bit versions of R.
    pairdist_ <- pairdist.ppp(X)

    # Estimation of K
    # NbPairs : number of pairs of points less than r apart
    #           (it's a vector, one value for each r)
    # pairdist_ > 0 eliminates distance from a point to itself
    # * 1.0 is to convert values to real and avoid infinite values
    # in NbPairs * w * l later
    NbPairs <- vapply(
      r,
      FUN = function(d) sum(pairdist_ > 0 & pairdist_ < d) * 1.0,
      FUN.VALUE = 0.0
    )
    # Kest gets the estimator of K, centered on the expected value.
    Kest <- NbPairs * w * l / (X$n * (X$n - 1)) - espKi_

    TestVector <- invsqrtmat(sigmaKi_) %*% Kest
    return(1 - stats::pchisq(sum(TestVector * TestVector), length(r)))
  }
  else return(NA)
}