R/gencrowd.R

In vincnardelli/crowdpostsampler: Spatial Crowdsourcing Sampling simulatior

#' Model Function
#'
#' This function allows you to express your love of cats.
#' @param NN number of units in the areas
#' @param nn number of sample units in the areas
#' @param lambda spatial autocorrelation variable
#'
#' @importFrom spdep dnearneigh nb2mat mat2listw moran.test
#' @importFrom graphics plot points
#' @importFrom stats rnorm runif var
#'
#' @export


gencrowd <- function(nn, NN, lambda){
  # Units observed in each stratum c1, c2, c3, c4 through crowdsourcing
  # they are given in this exercise, but in practice this should be evaluated ex-post
  # counting the units sampled in each stratum after allocating units with simple random sample
  # disregarding the different size of the strata
  n <- sum(nn)
  cc <- rep(n/4, 4)
  N <- sum(NN)
  N1 <- NN[1]
  N2 <- NN[2]
  N3 <- NN[3]
  N4 <- NN[4]
  quadrants <- list(
    list(c(-0.5, 0), c(-.5, 0)),
    list(c(-0.5, 0), c(0, 0.5)),
    list(c(0, 0.5), c(-0.5, 0)),
    list(c(0, 0.5), c(0, 0.5))
  )
  coord<-c()
  for(i in 1:4){
    xx1 <- runif(NN[i], quadrants[[i]][[1]][1], quadrants[[i]][[1]][2])
    xx2 <- runif(NN[i], quadrants[[i]][[2]][1], quadrants[[i]][[2]][2])
    coord <- rbind(coord, cbind(xx1, xx2))
  }

  # generation of a spatially correlated target variable y
  # W matrix is generated as a threshold distance. Threshold depends on the number of points
  lim<-200/N
  nb<-spdep::dnearneigh(coord,0,lim, longlat=FALSE)
  w<-spdep::nb2mat(nb)
  II<-array(1,N)
  I<-diag(II)
  G<-(I-lambda*w)
  # provisional population mean is MUP
  MUP<-10
  yy<-rnorm(N, MUP,0.1*MUP)
  eps<-rnorm(N,0,1)
  GINV<-solve(G)
  y<-GINV%*%yy+GINV%*%eps
  w_list<-mat2listw(w)
  # check of spatial correlation
  Moran_y<- spdep::moran.test(y,w_list)
  # ex-post true population mean is MU. This is the value that we want to estimate through sample observations.
  MU<-mean(y)
  #creating 4 sets of data for each stratum
  y_list <- list(y[1:N1], y[(N1+1):(N1+N2)], y[(N1+N2+1):(N1+N2+N3)], y[(N1+N2+N3+1):(N1+N2+N3+N4)])

  out <- list(
    "NN"= NN,
    "nn"= nn,
    "lambda" = lambda,
    "moran" = Moran_y,
    "MU" = MU,
    "y" = y_list,
    "coord" = coord

  )
  class(out) <- "crowd"
  return(out)
}