R/gk.R

In vincnardelli/robspdep: robspdep: Robust Spatial Dependence

#' Permutation test for GK statistic
#'
#'
#'  A permutation test for GK statistic calculated by using nsim random
#'  permutations of x for the given spatial weighting scheme generated with spdep,
#'  to establish the rank of the observed statistic in relation to the nsim simulated values.
#'
#' @param x a numeric vector the same length as the neighbours list in spdep listw
#' @param listw a \code{listw} object created with spdep for example by \code{spdep::nb2listw}
#' @param alternative alternative a character string specifying the alternative hypothesis, must be "greater" (default), or "less".
#' @param nsim  number of permutations
#'
#' @return  A list with \code{statistic}, \code{p.value} and \code{res}
#' @export
#' @import boot
#' @importFrom stats mad median punif
#'
#' @examples
#' gk(sampledata$x,sampledata$w)


gk <- function(x, listw, alternative = "greater", nsim=999){

  gk_boot <- function(x, i, x_lag, listw, nsim=nsim, ...) {
    x <- x[i]
    a<-1/mad(x)
    x_lag<- spdep::lag.listw(x=listw,var=x)
    b<-1/mad(x_lag)
    madsum <- mad(a*x+b*x_lag)^2
    maddif <- mad(a*x-b*x_lag)^2
    num<-madsum-maddif
    den<-madsum+maddif
    return(num/den)
  }


  x_lag<-spdep::lag.listw(x=listw,var=x)

  res3 <- boot::boot(x, statistic = gk_boot,
                     R = nsim,
                     sim = "permutation",
                     listw = listw,
                     x_lag=x_lag)


  a<-1/mad(x)
  b<-1/mad(x_lag)
  num<-(mad(a*x+b*x_lag)^2)-(mad(a*x-b*x_lag)^2)
  den<-(mad(a*x+b*x_lag)^2)+(mad(a*x-b*x_lag)^2)
  true <- num/den
  true

  res <- res3$t
  res[nsim + 1] <- true
  rankres <- rank(res)
  xrank <- rankres[length(res)]
  diff <- nsim - xrank
  diff <- ifelse(diff > 0, diff, 0)
  if (alternative == "less")
    pval <- punif((diff + 1)/(nsim + 1), lower.tail = FALSE)
  else if (alternative == "greater")
    pval <- punif((diff + 1)/(nsim + 1))

  statistic <- res[nsim + 1]
  statistic
  pval

  ret <- list(statistic = statistic,
              p.value = pval,
              res = res)
  class(ret) <- "robspdep.gk"

  return(ret)
}