R/IDS.R

In rbchan/unmarked: Models for Data from Unmarked Animals

setClassUnion("unmarkedFrameOrNULL", members=c("unmarkedFrame", "NULL"))

setClass("unmarkedFitIDS",
    representation(
        formlist = "list",
        keyfun = "character",
        K = "numeric",
        dataPC = "unmarkedFrameOrNULL",
        dataOC = "unmarkedFrameOrNULL",
        maxDist = "list",
        surveyDurations = "list",
        unitsOut = "character"),
        contains = "unmarkedFit")

get_ds_info <- function(db){
  J <- length(db) - 1
  a <- u <- rep(NA, J)
  a[1] <- pi*db[2]^2
  if(J > 1){
    for (j in 2:J){
      a[j] <- pi*db[j+1]^2 - sum(a[1:(j-1)])
    }
  }

  total_area <- sum(a)
  u <- a/total_area
  w <- diff(db)
  list(total_area=total_area, a=a, u=u, w=w)
}

IDS <- function(lambdaformula = ~1,
                detformulaDS = ~1,
                detformulaPC = NULL,
                detformulaOC = NULL,
                dataDS,
                dataPC = NULL,
                dataOC = NULL,
                availformula = NULL,
                durationDS = NULL,
                durationPC = NULL,
                durationOC = NULL,
                keyfun = "halfnorm",
                maxDistPC,
                maxDistOC,
                K = 100,
                unitsOut = "ha",
                starts = NULL,
                method = "BFGS",
                ...
                ){

  # Process inputs-------------------------------------------------------------

  form_hds <- as.formula(paste(c(as.character(detformulaDS),
                                 as.character(lambdaformula)), collapse=""))

  if(is.null(detformulaPC)){
    form_pc <- as.formula(paste(c(as.character(detformulaDS),
                                 as.character(lambdaformula)), collapse=""))
  } else{
    form_pc <- as.formula(paste(c(as.character(detformulaPC),
                                 as.character(lambdaformula)), collapse=""))
  }

  if(is.null(detformulaOC)){
    form_oc <- as.formula(paste(c(as.character(detformulaDS),
                                 as.character(lambdaformula)), collapse=""))
  } else{
    form_oc <- as.formula(paste(c(as.character(detformulaOC),
                                 as.character(lambdaformula)), collapse=""))
  }

  formlist <- list(lam=lambdaformula, ds=form_hds, pc=form_pc, oc=form_oc,
                   phi=availformula)

  stopifnot(inherits(dataDS, "unmarkedFrameDS"))
  stopifnot(inherits(dataPC, c("unmarkedFramePCount", "NULL")))
  stopifnot(inherits(dataOC, c("unmarkedFrameOccu", "NULL")))
  #stopifnot(!is.null(dataPC) | !is.null(dataOC))

  has_avail <- FALSE
  if(!is.null(durationDS) | !is.null(durationPC) | !is.null(durationOC)){
    has_avail <- TRUE
    if(is.null(availformula)) availformula <- ~1
    form_avail <- as.formula(paste("~1", paste(as.character(availformula), collapse="")))
    stopifnot(!is.null(durationDS))
    if(!is.null(dataPC)) stopifnot(!is.null(durationPC))
    if(!is.null(dataOC)) stopifnot(!is.null(durationOC))
  }

  if(has_avail & !is.null(dataOC)){
    stop("Availability estimation doesn't work with detection-nondetection data", call.=FALSE)
  }

  stopifnot(is.null(durationDS) || (length(durationDS) == numSites(dataDS)))
  stopifnot(is.null(durationPC) || (length(durationPC) == numSites(dataPC)))
  stopifnot(is.null(durationOC) || (length(durationOC) == numSites(dataOC)))
  surveyDurations <- list(ds=durationDS, pc=durationPC, oc=durationOC)

  stopifnot(keyfun %in% c("halfnorm", "exp"))
  keyidx <- switch(keyfun, "halfnorm"={1}, "exp"={2})

  if(missing(maxDistPC)) maxDistPC <- max(dataDS@dist.breaks)
  if(missing(maxDistOC)) maxDistOC <- max(dataDS@dist.breaks)


  # Design matrices------------------------------------------------------------

  # Need to add offset support here eventually
  gd_hds <- getDesign(dataDS, form_hds)
  if(any(is.na(gd_hds$y))){
    stop("Missing values in only some distance bins is not supported", call.=FALSE)
  } 
  ds_hds <- get_ds_info(dataDS@dist.breaks)
  Xavail_ds <- matrix(0,0,0)
  if(has_avail) Xavail_ds <- getDesign(dataDS, form_avail)$X
  if(is.null(durationDS)) durationDS <- rep(0,0)

  gd_pc <- list(y=matrix(0,0,0), X=matrix(0,0,0), V=matrix(0,0,0))
  ds_pc <- list(total_area=0, db=c(0,0), a=0, w=0, u=0)
  Xavail_pc <- matrix(0,0,0)
  if(is.null(durationPC)) durationPC <- rep(0,0)
  if(!is.null(dataPC)){
    gd_pc <- getDesign(dataPC, form_pc)
    ds_pc <- get_ds_info(c(0, maxDistPC))
    if(has_avail) Xavail_pc <- getDesign(dataPC, form_avail)$X
  }

  gd_oc <- list(y=matrix(0,0,0), X=matrix(0,0,0), V=matrix(0,0,0))
  ds_oc <- list(total_area=0, db=c(0,0), a=0, w=0, u=0)
  Kmin_oc <- rep(0,0)
  Xavail_oc <- matrix(0,0,0)
  if(is.null(durationOC)) durationOC <- rep(0,0)
  if(!is.null(dataOC)){
    gd_oc <- getDesign(dataOC, form_oc)
    ds_oc <- get_ds_info(c(0, maxDistOC))
    Kmin_oc <- apply(gd_oc$y, 1, max, na.rm=T)
    if(has_avail) Xavail_oc <- getDesign(dataOC, form_avail)$X
  }

  # Density conversion and unequal area correction
  lam_adjust <- c(ds_hds$total_area, ds_pc$total_area, ds_oc$total_area)
  names(lam_adjust) <- c("hds", "pc", "oc")

  switch(dataDS@unitsIn,
    m = lam_adjust <- lam_adjust / 1e6,
    km = lam_adjust <- lam_adjust)

  stopifnot(unitsOut %in% c("m","ha","kmsq"))
  switch(unitsOut,
    m = lam_adjust <- lam_adjust * 1e6,
    ha = lam_adjust <- lam_adjust * 100,
    kmsq = lam_adjust <- lam_adjust)

  # Parameter stuff------------------------------------------------------------
  # Doesn't support hazard
  pind_mat <- matrix(0, nrow=5, ncol=2)
  pind_mat[1,] <- c(1, ncol(gd_hds$X))
  pind_mat[2,] <- max(pind_mat) + c(1, ncol(gd_hds$V))
  if(!is.null(detformulaPC) & !is.null(dataPC)){
    pind_mat[3,] <- max(pind_mat) + c(1, ncol(gd_pc$V))
  }
  if(!is.null(detformulaOC) & !is.null(dataOC)){
    pind_mat[4,] <- max(pind_mat) + c(1, ncol(gd_oc$V))
  }
  if(has_avail){
    pind_mat[5,] <- max(pind_mat) + c(1, ncol(Xavail_ds))
  }

  if(is.null(starts)){
    lam_init <- log(mean(apply(dataDS@y, 1, sum, na.rm=TRUE)) / lam_adjust[1])
    params_tmb <- list(beta_lam = c(lam_init, rep(0, ncol(gd_hds$X)-1)),
                     beta_hds = c(log(median(dataDS@dist.breaks)),rep(0, ncol(gd_hds$V)-1)),
                     beta_pc = rep(0,0),
                     beta_oc = rep(0,0),
                     beta_avail = rep(0,0))

    if(!is.null(detformulaPC) & !is.null(dataPC)){
      params_tmb$beta_pc <- c(log(maxDistPC/2), rep(0, ncol(gd_pc$V)-1))
    }
    if(!is.null(detformulaOC) & !is.null(dataOC)){
      params_tmb$beta_oc <- c(log(maxDistOC/2), rep(0, ncol(gd_oc$V)-1))
    }
    if(has_avail){
      params_tmb$beta_avail <- rep(0, ncol(Xavail_ds))
    }
    starts <- unlist(params_tmb)
  } else {
    if(length(starts) != max(pind_mat)){
      stop("Length of starts should be ", max(pind_mat), call.=FALSE)
    }
    params_tmb <- list(beta_lam = starts[pind_mat[1,1]:pind_mat[1,2]],
                     beta_hds = starts[pind_mat[2,1]:pind_mat[2,2]],
                     beta_pc = rep(0,0),
                     beta_oc = rep(0,0),
                     beta_avail = rep(0,0))

    if(!is.null(detformulaPC) & !is.null(dataPC)){
      params_tmb$beta_pc <- starts[pind_mat[3,1]:pind_mat[3,2]]
    }
    if(!is.null(detformulaOC) & !is.null(dataOC)){
      params_tmb$beta_oc <- starts[pind_mat[4,1]:pind_mat[4,2]]
    }
    if(has_avail){
      params_tmb$beta_avail <- starts[pind_mat[5,1]:pind_mat[5,2]]
    }
  }

  # Construct TMB data list----------------------------------------------------
  tmb_dat <- list(
    # Common
    pind = pind_mat, lam_adjust = lam_adjust,

    # HDS data
    y_hds = gd_hds$y, X_hds = gd_hds$X, V_hds = gd_hds$V, key_hds = keyidx,
    db_hds = dataDS@dist.breaks, a_hds = ds_hds$a, w_hds = ds_hds$w,
    u_hds = ds_hds$u,

    # PC data
    y_pc = gd_pc$y, X_pc = gd_pc$X, V_pc = gd_pc$V, key_pc = keyidx,
    db_pc = c(0, maxDistPC), a_pc = ds_pc$a, w_pc = ds_pc$w, u_pc = ds_pc$u,

    # occ data
    y_oc = gd_oc$y, X_oc = gd_oc$X, V_oc = gd_oc$V, key_oc = keyidx,
    db_oc = c(0, maxDistOC), a_oc = ds_oc$a, w_oc = ds_oc$w, u_oc = ds_oc$u,
    K_oc = K, Kmin_oc = Kmin_oc,

    # avail data
    durationDS = durationDS, durationPC = durationPC, durationOC = durationOC,
    Xa_hds = Xavail_ds, Xa_pc = Xavail_pc, Xa_oc = Xavail_oc
  )

  tmb_obj <- TMB::MakeADFun(data = c(model = "tmb_IDS", tmb_dat), parameters = params_tmb,
                            DLL = "unmarked_TMBExports", silent=TRUE)

  opt <- optim(unlist(params_tmb), fn=tmb_obj$fn, gr=tmb_obj$gr, method=method, ...)

  fmAIC <- 2 * opt$value + 2 * length(unlist(params_tmb))

  sdr <- TMB::sdreport(tmb_obj)

  lam_coef <- get_coef_info(sdr, "lam", colnames(gd_hds$X),
                                       pind_mat[1,1]:pind_mat[1,2])

  lam_est <- unmarkedEstimate(name="Density", short.name="lam",
    estimates = lam_coef$ests, covMat = lam_coef$cov, fixed=1:ncol(gd_hds$X),
    invlink = "exp", invlinkGrad = "exp")

  dist_coef <- get_coef_info(sdr, "hds", colnames(gd_hds$V),
                                       pind_mat[2,1]:pind_mat[2,2])

  dist_est <- unmarkedEstimate(name="Distance sampling detection", short.name="ds",
    estimates = dist_coef$ests, covMat = dist_coef$cov, fixed=1:ncol(gd_hds$V),
    invlink = "exp", invlinkGrad = "exp")

  est_list <- list(lam=lam_est, ds=dist_est)

  if(!is.null(detformulaPC) & !is.null(dataPC)){
    pc_coef <- get_coef_info(sdr, "pc", colnames(gd_pc$V),
                                        pind_mat[3,1]:pind_mat[3,2])

    pc_est <- unmarkedEstimate(name="Point count detection", short.name="pc",
      estimates = pc_coef$ests, covMat = pc_coef$cov, fixed=1:ncol(gd_pc$V),
      invlink = "exp", invlinkGrad = "exp")
    est_list <- c(est_list, list(pc=pc_est))
  }

  if(!is.null(detformulaOC) & !is.null(dataOC)){
    oc_coef <- get_coef_info(sdr, "oc", colnames(gd_oc$V),
                                        pind_mat[4,1]:pind_mat[4,2])
    oc_est <- unmarkedEstimate(name="Presence/absence detection", short.name="oc",
      estimates = oc_coef$ests, covMat = oc_coef$cov, fixed=1:ncol(gd_oc$V),
      invlink = "exp", invlinkGrad = "exp")
    est_list <- c(est_list, list(oc=oc_est))
  }

  if(has_avail){
    avail_coef <- get_coef_info(sdr, "avail", colnames(Xavail_ds),
                                pind_mat[5,1]:pind_mat[5,2])
    avail_est <- unmarkedEstimate(name="Availability", short.name="phi",
      estimates=avail_coef$ests, covMat=avail_coef$cov, fixed=1:ncol(Xavail_ds),
      invlink="exp", invlinkGrad="exp")
    est_list <- c(est_list, list(phi=avail_est))
  }

  est_list <- unmarkedEstimateList(est_list)

  new("unmarkedFitIDS", fitType = "IDS", call = match.call(),
    opt = opt, formula = lambdaformula, formlist=formlist,
    data = dataDS, dataPC=dataPC, dataOC=dataOC,
    surveyDurations=surveyDurations,
    maxDist = list(pc=maxDistPC, oc=maxDistOC),
    keyfun=keyfun,
    # this needs to be fixed
    sitesRemoved = gd_hds$removed.sites,
    unitsOut=unitsOut,
    estimates = est_list, AIC = fmAIC, negLogLike = opt$value,
    nllFun = tmb_obj$fn, TMB=tmb_obj)

}

setMethod("summary", "unmarkedFitIDS", function(object)
{
    cat("\nCall:\n")
    print(object@call)
    cat("\n")
    summaryOut <- summary(object@estimates)
    cat("AIC:", object@AIC,"\n")
    cat("Number of distance sampling sites:", numSites(object@data))
    if(!is.null(object@dataPC)){
      cat("\nNumber of point count sites:", numSites(object@dataPC))
    }
    if(!is.null(object@dataOC)){
      cat("\nNumber of presence/absence sites:", numSites(object@dataOC))
    }
    cat("\noptim convergence code:", object@opt$convergence)
    cat("\noptim iterations:", object@opt$counts[1], "\n")
    if(!identical(object@opt$convergence, 0L))
    warning("Model did not converge. Try providing starting values or increasing maxit control argment.")
    cat("Bootstrap iterations:", length(object@bootstrapSamples), "\n\n")
    invisible(summaryOut)
})

# Need special method since an included dataset may not have a corresponding
# unique submodel in the unmarked estimates
setMethod("names", "unmarkedFitIDS", function(x){
  out <- c("lam","ds")
  if(!is.null(x@dataPC)) out <- c(out, "pc")
  if(!is.null(x@dataOC)) out <- c(out, "oc")
  if("phi" %in% names(x@estimates)) out <- c(out, "phi")
  out
})

# This function converts IDS objects into simpler distsamp objects
# so we can re-use distsamp methods
# the abundance model (lam) is combined with one of the detection models
# to yield a distsamp model
IDS_convert_class <- function(inp, type, ds_type=NULL){
  stopifnot(type %in% names(inp))
  if(is.null(ds_type)) ds_type <- type
  if(type == "lam") type <- "ds"
  data <- inp@data
  if(ds_type == "pc"){
    tempdat <- inp@dataPC
    maxDist <- inp@maxDist$pc
  }
  if(ds_type == "oc"){
    tempdat <- inp@dataOC
    maxDist <- inp@maxDist$oc
  }
  if(ds_type %in% c("pc", "oc")){
    if(is.null(maxDist)) maxDist <- max(data@dist.breaks)
    data <- unmarkedFrameDS(y=tempdat@y, siteCovs=siteCovs(tempdat),
                            dist.breaks=c(0, maxDist), survey="point",
                            unitsIn=data@unitsIn)
  }

  # Select appropriate model estimates; if sigma was not estimated
  # separately for this model, pick the DS model for detection
  det_mod <- type
  if(! det_mod %in% names(inp@estimates)) det_mod <- "ds"
  est <- inp@estimates@estimates[c("lam", det_mod)]
  names(est) <- c("state", "det")

  form <- inp@formlist[[type]]
  if(type=="phi") form <- as.formula(paste(c(as.character(form), "~1"), collapse=""))

  new("unmarkedFitDS", fitType="IDS", opt=inp@opt, formula=form,
      data=data, keyfun=inp@keyfun, unitsOut=inp@unitsOut,
      estimates=unmarkedEstimateList(est),
      AIC=inp@AIC, output="density", TMB=inp@TMB)
}

# This predict method uses IDS_convert_class to allow pass-through to
# distsamp predict method
setMethod("predict", "unmarkedFitIDS", function(object, type, newdata,
          backTransform=TRUE, appendData=FALSE, level=0.95, ...){
  stopifnot(type %in% names(object))

  # Special case of phi and  no newdata
  # We need a separate prediction for each detection dataset
  if(type == "phi" & missing(newdata)){

    dists <- names(object)[names(object) %in% c("ds", "pc", "oc")]
    out <- lapply(dists, function(x){
      conv <- IDS_convert_class(object, "phi", ds_type=x)
      predict(conv, "det", backTransform=backTransform, appendData=appendData,
              level=level, ...)
    })
    names(out) <- dists

  } else { # Regular situation
    conv <- IDS_convert_class(object, type)
    type <- switch(type, lam="state", ds="det", pc="det", oc="det", phi="det")
    out <- predict(conv, type, newdata, backTransform=backTransform, appendData=appendData,
                   level=level, ...)
  }
  out
})

# Get availability probability
setGeneric("getAvail", function(object, ...) standardGeneric("getAvail"))

# Get availability for each data type and site as a probability
setMethod("getAvail", "unmarkedFitIDS", function(object, ...){
  stopifnot("phi" %in% names(object))
  phi <- predict(object, "phi")
  dur <- object@surveyDurations
  out <- lapply(names(phi), function(x){
    1 - exp(-1 * dur[[x]] * phi[[x]]$Predicted)
  })
  names(out) <- names(phi)
  out
})

# Fitted method returns a list of matrices, one per data type
setMethod("fitted", "unmarkedFitIDS", function(object, na.rm=FALSE){

  dists <- names(object)[names(object) %in% c("ds", "pc")]

  # If there is an availability model, get availability
  # Otherwise set it to 1
  avail <- list(ds=1, pc=1, oc=1)
  if("phi" %in% names(object)){
    avail <- getAvail(object)
  }

  # fitted for distance and N-mix data components
  out <- lapply(dists, function(x){
    conv <- IDS_convert_class(object, type=x)
    fitted(conv) * avail[[x]]
  })
  names(out) <- dists

  # fitted for occupancy data
  if("oc" %in% names(object)){
    conv <- IDS_convert_class(object, type="oc")
    lam <- predict(conv, 'state')$Predicted
    A <- pi*max(conv@data@dist.breaks)^2
    switch(conv@data@unitsIn,
            m = A <- A / 1e6,
            km = A <- A)
    switch(conv@unitsOut,
            m = A <- A * 1e6,
            ha = A <- A * 100,
            kmsq = A <- A)
    lam <- lam * A

    p <- getP(conv) * avail$oc
    out$oc <- 1 - exp(-lam*p) ## analytical integration.
  }

  out
})

# getP returns detection probability WITHOUT availability
setMethod("getP", "unmarkedFitIDS", function(object, ...){

  dets <- names(object)[! names(object) %in% c("lam","phi")]

  out <- lapply(dets, function(x){
    conv <- IDS_convert_class(object, type=x)
    getP(conv)
  })
  names(out) <- dets
  out
})


setMethod("residuals", "unmarkedFitIDS", function(object, ...){

  dists <- names(object)[names(object) %in% c("ds", "pc")]

  # distance and N-mix data
  out <- lapply(dists, function(x){
    conv <- IDS_convert_class(object, type=x)
    residuals(conv)
  })
  names(out) <- dists

  # occupancy data
  if("oc" %in% names(object)){
    y <- object@dataOC@y
    ft <- fitted(object)$oc
    out$oc <- y - ft
  }

  out
})

setMethod("hist", "unmarkedFitIDS", function(x, lwd=1, lty=1, ...){

  conv <- IDS_convert_class(x, type='ds')
  hist(conv, lwd=lwd, lty=lty, ...)

})

setMethod("plot", c(x="unmarkedFitIDS", y="missing"), function(x, y, ...){

  r <- residuals(x)
  f <- fitted(x)
  nr <- length(r)
  long_names <- sapply(x@estimates@estimates, function(x) x@name)
  long_names <- long_names[long_names != "Density"]

  old_par <- graphics::par()$mfrow
  graphics::par(mfrow=c(nr,1))

  for (i in 1:nr){
    plot(f[[i]], r[[i]], ylab = "Residuals", xlab = "Predicted values",
         main=long_names[i])
    abline(h = 0, lty = 3, col = "gray")
  }
  graphics::par(mfrow=old_par)

})


setMethod("simulate", "unmarkedFitIDS",
          function(object,  nsim = 1, seed = NULL, na.rm = FALSE){

  dets <- c("ds","pc","oc")

  if("phi" %in% names(object)) avail <- getAvail(object)

  temp <- lapply(dets, function(x){
    if(! x %in% names(object)) return(NULL)
    conv <- IDS_convert_class(object, type=x)
    sims <- simulate(conv, nsim=nsim, na.rm=na.rm)
    # availability process
    if("phi" %in% names(object)){
      sims <- lapply(sims, function(z){
        avail_expand <- matrix(rep(avail[[x]], ncol(z)), ncol=ncol(z))
        out <- rbinom(length(z), z, avail_expand)
        matrix(out, nrow=nrow(z), ncol=ncol(z))
      })
    }
    if(x=="oc"){
      sims <- lapply(sims, function(z){
                    z[z>1] <- 1
                    z
                  })
    }
    sims
  })

  #"permute"
  lapply(1:nsim, function(i){
    sim <- lapply(temp, function(x) x[[i]])
    names(sim) <- c("ds","pc","oc")
    sim
  })

})

setMethod("update", "unmarkedFitIDS",
  function(object, lambdaformula, detformulaDS, detformulaPC, detformulaOC,
           dataDS, dataPC, dataOC, ...){
    call <- object@call

    if(!missing(lambdaformula)){
      call[["lambdaformula"]] <- lambdaformula
    } else {
      call[["lambdaformula"]] <- object@formlist$lam
    }

    if(!missing(detformulaDS)){
      call[["detformulaDS"]] <- detformulaDS
    } else {
      call[["detformulaDS"]] <- split_formula(object@formlist$ds)[[1]]
    }

    if(!missing(detformulaPC)){
      call[["detformulaPC"]] <- detformulaPC
    } else if(!is.null(object@dataPC) & !is.null(call$detformulaPC)){
      call[["detformulaPC"]] <- split_formula(object@formlist$pc)[[1]]
    }

    if(!missing(detformulaOC)){
      call[["detformulaOC"]] <- detformulaOC
    } else if(!is.null(object@dataOC) & !is.null(call$detformulaOC)){
      call[["detformulaOC"]] <- split_formula(object@formlist$oc)[[1]]
    }

    if(!missing(dataDS)){
      call$dataDS <- dataDS
    } else {
      call$dataDS <- object@data
    }

    if(!missing(dataPC)){
      call$dataPC <- dataPC
    } else {
      call$dataPC <- object@dataPC
    }
    
    if(!missing(dataOC)){
      call$dataOC <- dataOC
    } else {
      call$dataOC <- object@dataOC
    }

    extras <- match.call(call=sys.call(-1),
                         expand.dots = FALSE)$...
    if (length(extras) > 0) {
        existing <- !is.na(match(names(extras), names(call)))
        for (a in names(extras)[existing])
            call[[a]] <- extras[[a]]
        if (any(!existing)) {
            call <- c(as.list(call), extras[!existing])
            call <- as.call(call)
            }
        }

    eval(call, parent.frame(2))

})


setMethod("parboot", "unmarkedFitIDS",
    function(object, statistic=SSE, nsim=10, ...)
{
    dots <- list(...)
    statistic <- match.fun(statistic)
    call <- match.call(call = sys.call(-1))
    starts <- as.numeric(coef(object))

    t0 <- statistic(object, ...)
    lt0 <- length(t0)
    t.star <- matrix(NA, nsim, lt0)
    #if(!missing(report))
    #    cat("t0 =", t0, "\n")

    simList <- simulate(object, nsim = nsim, na.rm = FALSE)

    dataDS <- object@data
    dataPC <- object@dataPC
    has_pc <- "pc" %in% names(object)
    dataOC <- object@dataOC
    has_oc <- "oc" %in% names(object)

    t.star <- lapply(1:nsim, function(i){
      dataDS@y <- simList[[i]]$ds
      if(has_pc) dataPC@y <- simList[[i]]$pc
      if(has_oc) dataOC@y <- simList[[i]]$oc
      fit <- update(object, dataDS=dataDS, dataPC=dataPC, dataOC=dataOC, 
                    durationDS = object@surveyDurations$ds,
                    durationPC = object@surveyDurations$pc,
                    durationOC = object@surveyDurations$oc,
                    starts=starts)
      statistic(fit)
    })
    if(lt0 > 1){
      t.star <- t(t.star)
    } else {
      t.star <- matrix(t.star, ncol=lt0)
    }

    if (!is.null(names(t0))){
      colnames(t.star) <- names(t0)
    } else{
      colnames(t.star) <- paste("t*", 1:lt0, sep="")
    }

    out <- new("parboot", call = call, t0 = t0, t.star = t.star)
    return(out)
})

setMethod("SSE", "unmarkedFitIDS", function(fit, ...){
    r <- unlist(residuals(fit))
    return(c(SSE = sum(r^2, na.rm=T)))
})

setMethod("nonparboot", "unmarkedFitIDS",
    function(object, B = 0, keepOldSamples = TRUE, ...)
{
   stop("Not currently supported for unmarkedFitIDS", call.=FALSE)
})

setMethod("ranef", "unmarkedFitIDS",
    function(object, B = 0, keepOldSamples = TRUE, ...)
{
   stop("Not currently supported for unmarkedFitIDS", call.=FALSE)
})