R/occuRN.R

In unmarked: Models for Data from Unmarked Animals

# Fit the Occupancy model of Royle and Nichols

occuRN <- function(formula, data, K = 25, starts, method = "BFGS",
    se = TRUE, engine=c("C","R"), threads = 1, ...)
{
  formulas <- split_formula(formula)
  names(formulas) <- c("det", "state")
  check_no_support(formulas)

  if(!is(data, "unmarkedFrameOccu"))
    stop("Data is not an unmarkedFrameOccu object.")

  engine <- match.arg(engine, c("C", "R"))
  dm <- getDesign(data, formulas)
  y <- dm$y
  y <- truncateToBinary(y)

  J <- ncol(y)
  M <- nrow(y)

  occParms <- colnames(dm$X_state)
  detParms <- colnames(dm$X_det)
  nDP <- ncol(dm$X_det)
  nOP <- ncol(dm$X_state)

  nP <- nDP + nOP
  if(!missing(starts) && length(starts) != nP)
	   stop(paste("The number of starting values should be", nP))

  y.ji <- as.vector(y)
  navec <- is.na(y.ji)
  n <- 0:K

  nll_R <- function(parms, f = "Poisson")
  {

    ## compute individual level detection probabilities
    r.ij <- matrix(plogis(dm$X_det %*% parms[(nOP + 1) : nP] + dm$offset_det), M, J,
      byrow = TRUE)

    ## compute list of detection probabilities along N
    p.ij.list <- lapply(n, function(k) 1 - (1 - r.ij)^k)

    ## compute P(y_{ij} | N) (cell probabilities) along N
    cp.ij.list <- lapply(p.ij.list, function(pmat) pmat^y * (1-pmat)^(1-y))

    ## replace NA cell probabilities with 1.
    cp.ij.list <- lapply(cp.ij.list, function(cpmat) {
      cpmat[navec] <- 1
      cpmat
    })

    ## multiply across J to get P(y_i | N) along N
    cp.in <- sapply(cp.ij.list, rowProds)

    ## compute P(N = n | lambda_i) along i
    lambda.i <- exp(dm$X_state %*% parms[1 : nOP] + dm$offset_state)
    lambda.in <- sapply(n, function(x) dpois(x, lambda.i))

    ## integrate over P(y_i | N = n) * P(N = n | lambda_i) wrt n
    like.i <- rowSums(cp.in * lambda.in)

    -sum(log(like.i))
  }

  if(engine=="R"){
    nll <- nll_R
  } else{
    n_param <- c(nOP, nDP)
    Kmin <- apply(y, 1, function(x) max(x, na.rm=TRUE))
    nll <- function(params){
      nll_occuRN(params, n_param, y, dm$X_state, dm$X_det, dm$offset_state, dm$offset_det,
                 K, Kmin, threads)
    }
  }

	if(missing(starts)) starts <- rep(0, nP)
  fm <- optim(starts, nll, method = method, hessian = se, ...)
  covMat <- invertHessian(fm, nP, se)
  ests <- fm$par
  fmAIC <- 2 * fm$value + 2 * nP # + 2 * nP * (nP + 1) / (M - nP - 1)
  names(ests) <- c(occParms, detParms)

  stateEstimates <- unmarkedEstimate(name = "Abundance",
      short.name = "lam",
      estimates = ests[1:nOP],
      covMat = as.matrix(covMat[1:nOP,1:nOP]), invlink = "exp",
      invlinkGrad = "exp")

  detEstimates <- unmarkedEstimate(name = "Detection", short.name = "p",
      estimates = ests[(nOP + 1) : nP],
      covMat = as.matrix(covMat[(nOP + 1) : nP, (nOP + 1) : nP]),
      invlink = "logistic", invlinkGrad = "logistic.grad")

  estimateList <- unmarkedEstimateList(list(state=stateEstimates,
          det=detEstimates))

  umfit <- new("unmarkedFitOccuRN", fitType = "occuRN",
      call = match.call(), formula = formula, formlist = formulas, data = data,
      sitesRemoved = dm$removed.sites, estimates = estimateList,
      AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, K = K)

  return(umfit)
}