R/deprecated_sim_power.R

In rbchan/unmarked: Models for Data from Unmarked Animals

get_vars <- function(inp){
  if(is.list(inp)){
    out <- unique(unlist(lapply(inp, all.vars)))
  } else {
    out <- all.vars(inp)
  }
  names(out) <- out
  out
}

var_data <- function(var, guide, n){
  out <- rep(NA, n)
  gv <- guide[[var]]
  if(is.null(gv)){
    out <- stats::rnorm(n, 0, 1)
  } else if(inherits(gv, "factor")){
    levs <- levels(gv)
    out <- factor(sample(levs, n, replace=TRUE), levels=levs)
  } else{
    gv$n <- n
    out <- do.call(gv$dist, gv[!names(gv)=="dist"])
  }
  out
}

generate_data <- function(formulas, guide, n){
  vars <- get_vars(formulas)
  if(length(vars)==0) return(NULL)
  as.data.frame(lapply(vars, var_data, guide=guide, n=n))
}

capitalize <- function(inp){
  paste0(toupper(substring(inp,1,1)),
         substring(inp,2,nchar(inp)))
}

parse_func_name <- function(inp){
  if(!is.character(inp)){
    stop("Argument must be a character string", call.=FALSE)
  }
  capitalize(inp)
}

blank_umFit <- function(fit_function){
  type <- parse_func_name(fit_function)
  type <- ifelse(type=="Pcount", "PCount", type)
  type <- ifelse(type=="MultinomPois", "MPois", type)
  type <- ifelse(type=="Distsamp", "DS", type)
  type <- ifelse(type=="Colext", "ColExt", type)
  type <- ifelse(type=="Gdistsamp", "GDS", type)
  type <- ifelse(type=="Gpcount", "GPC", type)
  type <- ifelse(type=="Gmultmix", "GMM", type)
  type <- ifelse(type=="PcountOpen", "PCO", type)
  type <- ifelse(type=="DistsampOpen", "DSO", type)
  type <- ifelse(type=="MultmixOpen", "MMO", type)
  type <- ifelse(type=="Gdistremoval", "GDR", type)
  type <- paste0("unmarkedFit", type)
  new(type)
}


setMethod("simulate", "character",
  function(object, nsim=1, seed=NULL, formulas, coefs=NULL, design, guide=NULL, ...){

  .Deprecated("simulate", package=NULL, 
              msg = paste("Supplying the name of an unmarked fitting function to simulate will soon be removed from unmarked.\n",
                          "Use the simulate method for an unmarkedFrame instead. See the simulation vignette for more."),
             old = as.character(sys.call(sys.parent()))[1L])
  model <- blank_umFit(object)
  fit <- suppressWarnings(simulate_fit(model, formulas, guide, design, ...))
  coefs <- check_coefs_old(coefs, fit)
  #fit <- replace_sigma(coefs, fit)
  coefs <- generate_random_effects(coefs, fit)
  fit <- replace_estimates(fit, coefs)
  ysims <- suppressWarnings(simulate(fit, nsim))
  umf <- fit@data
  # fix this
  umfs <- lapply(ysims, function(x){
    if(object=="occuMulti"){
      umf@ylist <- x
    } else if(object=="gdistremoval"){
      umf@yDistance=x$yDistance
      umf@yRemoval=x$yRemoval
    } else if(object == "IDS"){
       out <- list()
       out$ds <- fit@data
       out$ds@y <- x$ds
       if("pc" %in% names(fit)){
         out$pc <- fit@dataPC
         out$pc@y <- x$pc
       }
       if("oc" %in% names(fit)){
         out$oc <- fit@dataOC
         out$oc@y <- x$oc
       }
       umf <- out
    } else {
      umf@y <- x
    }
    umf
  })
  if(length(umfs)==1) umfs <- umfs[[1]]
  umfs
})

# Insert specified random effects SD into proper S4 slot in model object
# This is mostly needed by GDR which uses the SD to calculate
# N with E_loglam (this is currently disabled so the function is not needed)
#replace_sigma <- function(coefs, fit){
#  required_subs <- names(fit@estimates@estimates)
#  formulas <- sapply(names(fit), function(x) get_formula(fit, x))
#  rand <- lapply(formulas, lme4::findbars)
#  if(!all(sapply(rand, is.null))){
#    rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars)))
#    for (i in required_subs){
#      if(!is.null(rand[[i]][[1]])){
#        signame <- rvar[[i]]
#        old_coefs <- coefs[[i]]
#        fit@estimates@estimates[[i]]@randomVarInfo$estimates <- coefs[[i]][[signame]]
#      }
#    }
#  }
#  fit
#}



setGeneric("get_umf_components", function(object, ...) standardGeneric("get_umf_components"))

setMethod("get_umf_components", "unmarkedFit",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$state, guide, design$M)
  oc <- generate_data(formulas$det, guide, design$J*design$M)
  yblank <- matrix(0, design$M, design$J)
  list(y=yblank, siteCovs=sc, obsCovs=oc)
})


setGeneric("simulate_fit", function(object, ...) standardGeneric("simulate_fit"))

setMethod("simulate_fit", "unmarkedFitOccu",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  umf <- unmarkedFrameOccu(y=parts$y, siteCovs=parts$siteCovs,
                           obsCovs=parts$obsCovs)
  occu(as.formula(paste(deparse(formulas$det), deparse(formulas$state))),
       data=umf, se=FALSE, control=list(maxit=1))
})

setMethod("simulate_fit", "unmarkedFitPCount",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  umf <- unmarkedFramePCount(y=parts$y, siteCovs=parts$siteCovs,
                             obsCovs=parts$obsCovs)
  args <- list(...)
  K <- ifelse(is.null(args$K), 100, args$K)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  pcount(as.formula(paste(deparse(formulas$det), deparse(formulas$state))),
       data=umf, mixture=mixture, K=K, se=FALSE, control=list(maxit=1))
})

setMethod("simulate_fit", "unmarkedFitOccuRN",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  umf <- unmarkedFrameOccu(y=parts$y, siteCovs=parts$siteCovs,
                           obsCovs=parts$obsCovs)
  occuRN(as.formula(paste(deparse(formulas$det), deparse(formulas$state))),
       data=umf, se=FALSE, control=list(maxit=1))
})

setMethod("get_umf_components", "unmarkedFitMPois",
          function(object, formulas, guide, design, ...){
  args <- list(...)
  sc <- generate_data(formulas$state, guide, design$M)
  oc <- generate_data(formulas$det, guide, design$J*design$M)
  J <- ifelse(args$type=="double", 3, design$J)
  yblank <- matrix(0, design$M, design$J)
  list(y=yblank, siteCovs=sc, obsCovs=oc)
})

setMethod("simulate_fit", "unmarkedFitMPois",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  type <- ifelse(is.null(args$type), "removal", args$type)
  umf <- unmarkedFrameMPois(y=parts$y, siteCovs=parts$siteCovs,
                           obsCovs=parts$obsCovs, type=type)
  multinomPois(as.formula(paste(deparse(formulas$det), deparse(formulas$state))),
               data=umf, se=FALSE, control=list(maxit=1))
})

setMethod("get_umf_components", "unmarkedFitDS",
          function(object, formulas, guide, design, ...){
  #args <- list(...)
  sc <- generate_data(formulas$state, guide, design$M)
  sc2 <- generate_data(formulas$det, guide, design$M)
  dat <- list(sc, sc2)
  keep <- sapply(dat, function(x) !is.null(x))
  dat <- dat[keep]
  sc <- do.call(cbind, dat)
  yblank <- matrix(0, design$M, design$J)
  list(y=yblank, siteCovs=sc)
})

setMethod("simulate_fit", "unmarkedFitDS",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(is.null(args$tlength)) args$tlength <- 0
  umf <- unmarkedFrameDS(y=parts$y, siteCovs=parts$siteCovs,
                         tlength=args$tlength, survey=args$survey, unitsIn=args$unitsIn,
                         dist.breaks=args$dist.breaks)
  keyfun <- ifelse(is.null(args$keyfun), "halfnorm", args$keyfun)
  output <- ifelse(is.null(args$output), "density", args$output)
  unitsOut <- ifelse(is.null(args$unitsOut), "ha", args$unitsOut)

  distsamp(as.formula(paste(deparse(formulas$det), deparse(formulas$state))),
               data=umf, se=FALSE, control=list(maxit=1), keyfun=keyfun,
          output=output, unitsOut=unitsOut)
})


setMethod("get_umf_components", "unmarkedFitColExt",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$psi, guide, design$M)
  ysc <- generate_data(list(formulas$col, formulas$ext), guide, design$M*design$T)
  oc <- generate_data(formulas$det, guide, design$J*design$M*design$T)
  yblank <- matrix(0, design$M, design$T*design$J)
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc, obsCovs=oc)
})


setMethod("simulate_fit", "unmarkedFitColExt",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  umf <- unmarkedMultFrame(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           obsCovs=parts$obsCovs, numPrimary=design$T)
  colext(psiformula=formulas$psi, gammaformula=formulas$col,
         epsilonformula=formulas$ext,pformula=formulas$det,
         data=umf, se=FALSE, control=list(maxit=1))
})

setMethod("get_umf_components", "unmarkedFitOccuTTD",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$psi, guide, design$M)
  ysc <- NULL
  if(design$T>1){
    ysc <- generate_data(list(formulas$col, formulas$ext), guide, design$M*design$T)
  }
  oc <- generate_data(formulas$det, guide, design$J*design$M*design$T)
  yblank <- matrix(0, design$M, design$T*design$J)
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc, obsCovs=oc)
})


setMethod("simulate_fit", "unmarkedFitOccuTTD",
  function(object, formulas, guide, design, ...){
  if(is.null(design$T)) design$T <- 1
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  umf <- unmarkedFrameOccuTTD(y=parts$y,
                              surveyLength=args$surveyLength,
                              siteCovs=parts$siteCovs,
                              yearlySiteCovs=parts$yearlySiteCovs,
                              obsCovs=parts$obsCovs, numPrimary=design$T)
  linkPsi <- ifelse(is.null(args$linkPsi), "logit", args$linkPsi)
  ttdDist <- ifelse(is.null(args$ttdDist), "exp", args$ttdDist)
  occuTTD(psiformula=formulas$psi, gammaformula=formulas$col,
         epsilonformula=formulas$ext,detformula=formulas$det,
         linkPsi=linkPsi, ttdDist=ttdDist,
         data=umf, se=FALSE, control=list(maxit=1))
})


setMethod("get_umf_components", "unmarkedFitGMM",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$lambda, guide, design$M)
  ysc <- generate_data(formulas$phi, guide, design$M*design$T)
  yblank <- matrix(0, design$M, design$T*design$J)
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc)
})


setMethod("simulate_fit", "unmarkedFitGDS",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(args$survey=="line"){
    umf <- unmarkedFrameGDS(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           numPrimary=design$T,
                           tlength=args$tlength, survey=args$survey,
                           unitsIn=args$unitsIn, dist.breaks=args$dist.breaks)
  } else if(args$survey=="point"){
    umf <- unmarkedFrameGDS(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           numPrimary=design$T, survey=args$survey,
                           unitsIn=args$unitsIn, dist.breaks=args$dist.breaks)
  }

  keyfun <- ifelse(is.null(args$keyfun), "halfnorm", args$keyfun)
  output <- ifelse(is.null(args$output), "density", args$output)
  unitsOut <- ifelse(is.null(args$unitsOut), "ha", args$unitsOut)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  K <- ifelse(is.null(args$K), 100, args$K)

  gdistsamp(lambdaformula=formulas$lambda, phiformula=formulas$phi,
            pformula=formulas$det, data=umf, keyfun=keyfun, output=output,
            unitsOut=unitsOut, mixture=mixture, K=K,
            se=FALSE, control=list(maxit=1))
})

setMethod("simulate_fit", "unmarkedFitGPC",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  umf <- unmarkedFrameGPC(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           numPrimary=design$T)
  K <- ifelse(is.null(args$K), 100, args$K)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)

  gpcount(lambdaformula=formulas$lambda, phiformula=formulas$phi,
          pformula=formulas$det, data=umf, mixture=mixture, K=K,
          se=FALSE, control=list(maxit=1))
})


setMethod("simulate_fit", "unmarkedFitGMM",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  umf <- unmarkedFrameGMM(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           numPrimary=design$T, type=args$type)
  K <- ifelse(is.null(args$K), 100, args$K)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)

  gmultmix(lambdaformula=formulas$lambda, phiformula=formulas$phi,
           pformula=formulas$det, data=umf, mixture=mixture, K=K,
           se=FALSE, control=list(maxit=1))
})


setMethod("get_umf_components", "unmarkedFitDailMadsen",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$lambda, guide, design$M)
  ysc <- generate_data(list(formulas$gamma, formulas$omega), guide, design$M*design$T)
  oc <- generate_data(formulas$det, guide, design$M*design$T*design$J)
  yblank <- matrix(0, design$M, design$T*design$J)
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc, obsCovs=oc)
})

setMethod("simulate_fit", "unmarkedFitPCO",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(is.null(args$primaryPeriod)){
    args$primaryPeriod <- matrix(1:design$T, design$M, design$T, byrow=TRUE)
  }
  umf <- unmarkedFramePCO(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           numPrimary=design$T, primaryPeriod=args$primaryPeriod)
  K <- ifelse(is.null(args$K), 100, args$K)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  dynamics <- ifelse(is.null(args$dynamics), "constant", args$dynamics)
  fix <- ifelse(is.null(args$fix), "none", args$fix)
  immigration <- ifelse(is.null(args$immigration), FALSE, args$immigration)
  iotaformula <- args$iotaformula
  if(is.null(iotaformula)) iotaformula <- ~1

  pcountOpen(lambdaformula=formulas$lambda, gammaformula=formulas$gamma,
             omegaformula=formulas$omega, pformula=formulas$det,
             data=umf, mixture=mixture, K=K, dynamics=dynamics, fix=fix,
             se=FALSE, method='SANN', control=list(maxit=1), immigration=immigration,
             iotaformula=iotaformula)
})

setMethod("simulate_fit", "unmarkedFitMMO",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(is.null(args$primaryPeriod)){
    args$primaryPeriod <- matrix(1:design$T, design$M, design$T, byrow=TRUE)
  }
  umf <- unmarkedFrameMMO(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           type=args$type,
                           numPrimary=design$T, primaryPeriod=args$primaryPeriod)
  K <- ifelse(is.null(args$K), 100, args$K)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  dynamics <- ifelse(is.null(args$dynamics), "constant", args$dynamics)
  fix <- ifelse(is.null(args$fix), "none", args$fix)
  immigration <- ifelse(is.null(args$immigration), FALSE, args$immigration)
  iotaformula <- args$iotaformula
  if(is.null(iotaformula)) iotaformula <- ~1

  multmixOpen(lambdaformula=formulas$lambda, gammaformula=formulas$gamma,
             omegaformula=formulas$omega, pformula=formulas$det,
             data=umf, mixture=mixture, K=K, dynamics=dynamics, fix=fix,
             se=FALSE, method='SANN', control=list(maxit=1), immigration=immigration,
             iotaformula=iotaformula)
})

setMethod("get_umf_components", "unmarkedFitDSO",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(formulas$lambda, guide, design$M)
  ysc <- generate_data(list(formulas$gamma, formulas$omega, formulas$det),
                       guide, design$M*design$T)
  yblank <- matrix(0, design$M, design$T*design$J)
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc)
})

setMethod("simulate_fit", "unmarkedFitDSO",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(is.null(args$primaryPeriod)){
    args$primaryPeriod <- matrix(1:design$T, design$M, design$T, byrow=TRUE)
  }
  if(args$survey=="line"){
  umf <- unmarkedFrameDSO(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,
                           tlength=args$tlength, survey=args$survey,
                           unitsIn=args$unitsIn, dist.breaks=args$dist.breaks,
                           numPrimary=design$T, primaryPeriod=args$primaryPeriod)
  } else if(args$survey == "point"){
  umf <- unmarkedFrameDSO(y=parts$y, siteCovs=parts$siteCovs,
                           yearlySiteCovs=parts$yearlySiteCovs,survey=args$survey,
                           unitsIn=args$unitsIn, dist.breaks=args$dist.breaks,
                           numPrimary=design$T, primaryPeriod=args$primaryPeriod)
  }
  K <- ifelse(is.null(args$K), 100, args$K)
  keyfun <- ifelse(is.null(args$keyfun), "halfnorm", args$keyfun)
  output <- ifelse(is.null(args$output), "density", args$output)
  unitsOut <- ifelse(is.null(args$unitsOut), "ha", args$unitsOut)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  dynamics <- ifelse(is.null(args$dynamics), "constant", args$dynamics)
  fix <- ifelse(is.null(args$fix), "none", args$fix)
  immigration <- ifelse(is.null(args$immigration), FALSE, args$immigration)
  iotaformula <- args$iotaformula
  if(is.null(iotaformula)) iotaformula <- ~1
  distsampOpen(lambdaformula=formulas$lambda, gammaformula=formulas$gamma,
             omegaformula=formulas$omega, pformula=formulas$det,
             keyfun=keyfun, unitsOut=unitsOut, output=output,
             data=umf, mixture=mixture, K=K, dynamics=dynamics, fix=fix,
             se=FALSE, method='SANN', control=list(maxit=1), immigration=immigration,
             iotaformula=iotaformula)
})


setMethod("get_umf_components", "unmarkedFitOccuMulti",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(lapply(formulas$state, as.formula), guide, design$M)
  oc <- generate_data(lapply(formulas$det, as.formula), guide, design$J*design$M)
  nspecies <- length(formulas$det)
  yblank <- lapply(1:nspecies, function(x) matrix(0, design$M, design$J))
  list(y=yblank, siteCovs=sc, obsCovs=oc)
})

setMethod("simulate_fit", "unmarkedFitOccuMulti",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  if(is.null(args$maxOrder)) args$maxOrder <- length(parts$y)
  umf <- unmarkedFrameOccuMulti(y=parts$y, siteCovs=parts$siteCovs,
                                obsCovs=parts$obsCovs, maxOrder=args$maxOrder)
  occuMulti(formulas$det, formulas$state, data=umf, maxOrder=args$maxOrder,
            se=FALSE, control=list(maxit=1))
})

setMethod("get_umf_components", "unmarkedFitOccuMS",
          function(object, formulas, guide, design, ...){
  sc <- generate_data(lapply(formulas$state, as.formula), guide, design$M)
  ysc <- NULL
  if(!is.null(formulas$phi)){
    ysc <- generate_data(lapply(formulas$phi, as.formula), guide, design$M*design$T*design$J)
  }
  oc <- generate_data(lapply(formulas$det, as.formula), guide, design$J*design$M)
  nspecies <- length(formulas$det)
  yblank <- matrix(0, design$M, design$T*design$J)
  yblank[1,1] <- 2 # To bypass sanity checker in unmarkedFrameOccuMS
  list(y=yblank, siteCovs=sc, yearlySiteCovs=ysc, obsCovs=oc)
})

setMethod("simulate_fit", "unmarkedFitOccuMS",
  function(object, formulas, guide, design, ...){
  if(is.null(design$T)) design$T <- 1
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  umf <- unmarkedFrameOccuMS(y=parts$y, siteCovs=parts$siteCovs,
                                yearlySiteCovs=parts$yearlySiteCovs,
                                obsCovs=parts$obsCovs, numPrimary=design$T)
  if(is.null(args$parameterization)) args$parameterization <- "multinomial"
  occuMS(formulas$det, formulas$state, formulas$phi, data=umf,
         parameterization=args$parameterization,
         se=FALSE, control=list(maxit=1))
})

setMethod("get_umf_components", "unmarkedFitGDR",
          function(object, formulas, guide, design, ...){
  if(any(! c("M","Jdist","Jrem") %in% names(design))){
    stop("Required design components are M, Jdist, and Jrem")
  }
  sc <- generate_data(list(formulas$lambda, formulas$dist), guide, design$M)
  ysc <- NULL
  if(design$T > 1){
    ysc <- generate_data(formulas$phi, guide, design$M*design$T)
  }
  oc <- generate_data(formulas$rem, guide, design$M*design$T*design$Jrem)

  list(yDistance=matrix(0, design$M, design$T*design$Jdist),
       yRemoval=matrix(0, design$M, design$T*design$Jrem),
       siteCovs=sc, yearlySiteCovs=ysc, obsCovs=oc)
})

setMethod("simulate_fit", "unmarkedFitGDR",
  function(object, formulas, guide, design, ...){
  if(is.null(design$T)) design$T <- 1
  if(is.null(formulas$phi)) formulas$phi <- ~1
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)
  umf <- unmarkedFrameGDR(yDistance=parts$yDistance, yRemoval=parts$yRemoval,
                          numPrimary=design$T, siteCovs=parts$siteCovs,
                          obsCovs=parts$obsCovs, yearlySiteCovs=parts$yearlySiteCovs,
                          dist.breaks=args$dist.breaks, unitsIn=args$unitsIn,
                          period.lengths=args$period.lengths)

  keyfun <- ifelse(is.null(args$keyfun), "halfnorm", args$keyfun)
  output <- ifelse(is.null(args$output), "density", args$output)
  unitsOut <- ifelse(is.null(args$unitsOut), "ha", args$unitsOut)
  mixture <- ifelse(is.null(args$mixture), "P", args$mixture)
  K <- ifelse(is.null(args$K), 100, args$K)

  gdistremoval(lambdaformula=formulas$lambda, phiformula=formulas$phi,
               removalformula=formulas$rem, distanceformula=formulas$dist,
               data=umf, keyfun=keyfun, output=output, unitsOut=unitsOut,
               mixture=mixture, K=K, se=FALSE, control=list(maxit=1), method='L-BFGS-B')
})

# For simulating entirely new datasets
setMethod("get_umf_components", "unmarkedFitIDS",
          function(object, formulas, guide, design, ...){

  # Distance sampling dataset
  sc_ds_lam <- generate_data(formulas$lam, guide, design$Mds)
  sc_ds_det <- generate_data(formulas$ds, guide, design$Mds)
  dat_ds <- list(sc_ds_lam, sc_ds_det)
  if(!is.null(formulas$phi)){
    sc_ds_phi <- generate_data(formulas$phi, guide, design$Mds)
    dat_ds <- c(dat_ds, list(sc_ds_phi))
  }
  keep <- sapply(dat_ds, function(x) !is.null(x))
  dat_ds <- dat_ds[keep]
  sc_ds <- do.call(cbind, dat_ds)
  yblank_ds <- matrix(1, design$Mds, design$J)

  # Point count dataset
  sc_pc <- yblank_pc <- NULL
  if(!is.null(design$Mpc) && design$Mpc > 0){
    if(is.null(formulas$pc)) form_pc <- formulas$ds
    sc_pc_lam <- generate_data(formulas$lam, guide, design$Mpc)
    sc_pc_det <- generate_data(form_pc, guide, design$Mpc)
    sc_pc <- list(sc_pc_lam, sc_pc_det)
    if(!is.null(formulas$phi)){
      sc_pc_phi <- generate_data(formulas$phi, guide, design$Mpc)
      sc_pc <- c(sc_pc, list(sc_pc_phi))
    }
    keep <- sapply(sc_pc, function(x) !is.null(x))
    sc_pc <- sc_pc[keep]
    sc_pc <- do.call(cbind, sc_pc)
    yblank_pc <- matrix(1, design$Mpc, 1)
  }

  # Presence/absence dataset
  sc_oc <- yblank_oc <- NULL
  if(!is.null(design$Moc) && design$Moc > 0){
    if(is.null(formulas$oc)){
      form_oc <- formulas$ds
    } else {
      form_oc <- formulas$oc
    }
    sc_oc_lam <- generate_data(formulas$lam, guide, design$Moc)
    sc_oc_det <- generate_data(form_oc, guide, design$Moc)
    sc_oc <- list(sc_oc_lam, sc_oc_det)
    if(!is.null(formulas$phi)){
      sc_oc_phi <- generate_data(formulas$phi, guide, design$Moc)
      sc_oc <- c(sc_oc, list(sc_oc_phi))
    }
    keep <- sapply(sc_oc, function(x) !is.null(x))
    sc_oc <- sc_oc[keep]
    sc_oc <- do.call(cbind, sc_oc)
    yblank_oc <- matrix(1, design$Moc, 1)
  }

  mget(c("yblank_ds", "sc_ds", "yblank_pc", "sc_pc", "yblank_oc", "sc_oc"))
})


setMethod("simulate_fit", "unmarkedFitIDS",
  function(object, formulas, guide, design, ...){
  parts <- get_umf_components(object, formulas, guide, design, ...)
  args <- list(...)

  args$tlength <- 0
  args$survey <- "point"

  # Distance sampling dataset
  umf_ds <- unmarkedFrameDS(y=parts$yblank_ds, siteCovs=parts$sc_ds,
                            tlength=args$tlength, survey=args$survey,
                            unitsIn=args$unitsIn,
                            dist.breaks=args$dist.breaks)
  # Point count dataset
  umf_pc <- NULL
  if(!is.null(design$Mpc) && design$Mpc > 0){
    umf_pc <- unmarkedFramePCount(y=parts$yblank_pc, siteCovs=parts$sc_pc)
  }

  # Occupancy dataset
  umf_oc <- NULL
  if(!is.null(design$Moc) && design$Moc > 0){
    umf_oc <- unmarkedFrameOccu(y=parts$yblank_oc, siteCovs=parts$sc_oc)
  }

  keyfun <- ifelse(is.null(args$keyfun), "halfnorm", args$keyfun)
  unitsOut <- ifelse(is.null(args$unitsOut), "ha", args$unitsOut)
  K <- ifelse(is.null(args$K), 300, args$K)
  if(is.null(args$maxDistPC)) args$maxDistPC <- max(args$dist.breaks)
  if(is.null(args$maxDistOC)) args$maxDistOC <- max(args$dist.breaks)

  IDS(lambdaformula = formulas$lam,
      detformulaDS = formulas$ds,
      detformulaPC = formulas$pc, detformulaOC = formulas$oc,
      dataDS = umf_ds, dataPC = umf_pc, dataOC = umf_oc,
      availformula = formulas$phi,
      durationDS = args$durationDS, durationPC = args$durationPC,
      durationOC = args$durationOC,
      maxDistPC = args$maxDistPC, maxDistOC = args$maxDistOC,
      keyfun=keyfun, unitsOut=unitsOut, K=K ,control=list(maxit=1))
})


# power -----------------------------------------------------------------------

setClass("unmarkedPower_old",
  representation(call="call", data="unmarkedFrame", M="numeric",
                 J="numeric", T="numeric", coefs="list", estimates="list",
                 alpha="numeric", nulls="list")
)

setMethod("powerAnalysis", "unmarkedFit",
  function(object, coefs=NULL, design=NULL, alpha=0.05, nulls=list(),
                          datalist=NULL,
                          nsim=ifelse(is.null(datalist), 100, length(datalist)),
                          parallel=FALSE){
  
  .Deprecated("powerAnalysis", package=NULL, 
              msg = paste("Using an unmarkedFit object to run a power analysis will soon be removed from unmarked.\n",
                          "Use an unmarkedFrame instead. See the power analysis vignette for more."),
             old = as.character(sys.call(sys.parent()))[1L])

  submodels <- names(object@estimates@estimates)
  coefs <- check_coefs_old(coefs, object)
  coefs <- generate_random_effects(coefs, object)
  fit_temp <- replace_estimates(object, coefs)

  T <- 1
  bdata <- NULL
  if(!is.null(datalist)){
    if(length(datalist) != nsim){
      stop("Length of data list must equal value of nsim", call.=FALSE)
    }
    tryCatch({test <- update(object, data=datalist[[1]], se=FALSE,
                        control=list(maxit=1))
    }, error=function(e){
      stop("Incorrect format of entries in datalist", call.=FALSE)
    })
    bdata <- datalist
    M <- numSites(bdata[[1]])
    sims <- lapply(bdata, function(x){
      #fit_temp@data <- x
      #temporary workaround - not necessary??
      #if(methods::.hasSlot(fit_temp, "knownOcc")){
      #  fit_temp@knownOcc <- rep(FALSE, M)
      #}
      #simulate(fit_temp, 1)[[1]]
      if(inherits(x, "unmarkedFrameOccuMulti")){
        return(x@ylist)
      } else if(inherits(x, "unmarkedFrameGDR")){
        return(list(yDistance=x@yDistance, yRemoval=x@yRemoval))
      } else {
        return(x@y)
      }
    })
    if(methods::.hasSlot(bdata[[1]], "numPrimary")){
      T <- bdata[[1]]@numPrimary
    }
    J <- obsNum(bdata[[1]]) / T
  } else if(is.null(design)){
      sims <- simulate(fit_temp, nsim)
      M <- numSites(object@data)
      if(methods::.hasSlot(object@data, "numPrimary")){
        T <- object@data@numPrimary
      }
      J <- obsNum(object@data) / T
  } else {
    bdata <- bootstrap_data(fit_temp@data, nsim, design)
    sims <- lapply(bdata, function(x){
      fit_temp@data <- x
      #temporary workaround
      if(methods::.hasSlot(fit_temp, "knownOcc")){
        fit_temp@knownOcc <- rep(FALSE, design$M)
      }
      simulate(fit_temp, 1)[[1]]
    })
    M <- design$M
    if(methods::.hasSlot(fit_temp@data, "numPrimary")){
      T <- fit_temp@data@numPrimary
    }
    J <- design$J
  }

  cl <- NULL
  if(parallel){
    cl <- parallel::makeCluster(parallel::detectCores()-1)
    on.exit(parallel::stopCluster(cl))
    parallel::clusterEvalQ(cl, library(unmarked))
  }

  if(!is.null(options()$unmarked_shiny)&&options()$unmarked_shiny){
    ses <- options()$unmarked_shiny_session
    ses <- shiny::getDefaultReactiveDomain()
    pb <- shiny::Progress$new(ses, min=0, max=1)
    pb$set(message="Running simulations")
    if(!requireNamespace("pbapply", quietly=TRUE)){
      stop("You need to install the pbapply package", call.=FALSE)
    }
    fits <- pbapply::pblapply(1:nsim, function(i, sims, fit, bdata=NULL){
      if(!is.null(design)) fit@data <- bdata[[i]]
      if(inherits(fit, "unmarkedFitOccuMulti")){
        fit@data@ylist <- sims[[i]]
      } else{
        fit@data@y <- sims[[i]]
      }
      out <- update(fit, data=fit@data, se=TRUE)
      pb$set(value=i/nsim, message=NULL, detail=NULL)
      out
    }, sims=sims, fit=object, bdata=bdata, cl=NULL)
    pb$close()

  } else {

    fits <- lapply2(1:nsim, function(i, sims, fit, bdata=NULL){
      if(!is.null(design)) fit@data <- bdata[[i]]
      if(inherits(fit, "unmarkedFitOccuMulti")){
        fit@data@ylist <- sims[[i]]
      } else if(inherits(fit, "unmarkedFitGDR")){
        fit@data@yDistance <- sims[[i]]$yDistance
        fit@data@yRemoval <- sims[[i]]$yRemoval
      } else {
        fit@data@y <- sims[[i]]
      }
      update(fit, data=fit@data, se=TRUE)
    }, sims=sims, fit=object, bdata=bdata, cl=cl)

  }

  sum_dfs <- lapply(fits, get_summary_df_old)

  new("unmarkedPower_old", call=object@call, data=object@data, M=M,
      J=J, T=T, coefs=coefs, estimates=sum_dfs, alpha=alpha, nulls=nulls)
})

bootstrap_data <- function(data, nsims, design){
  M <- design$M
  J <- design$J
  sites <- 1:numSites(data)
  if(!is.null(J) & methods::.hasSlot(data, "numPrimary")){
    stop("Can't automatically bootstrap observations with > 1 primary period", call.=FALSE)
  }
  if(J > obsNum(data)){
    stop("Can't currently bootstrap more than the actual number of observations", call.=FALSE)
  }
  obs <- 1:obsNum(data)

  if(M > numSites(data)){
    M_samps <- lapply(1:nsims, function(i) sample(sites, M, replace=TRUE))
  } else if(M < numSites(data)){
    M_samps <- lapply(1:nsims, function(i) sample(sites, M, replace=FALSE))
  } else {
    M_samps <- replicate(nsims, sites, simplify=FALSE)
  }

  if(J > obsNum(data)){
    J_samps <- lapply(1:nsims, function(i) sample(obs, J, replace=TRUE))
  } else if(J < obsNum(data)){
    J_samps <- lapply(1:nsims, function(i) sample(obs, J, replace=FALSE))
  } else {
    J_samps <- replicate(nsims, obs, simplify=FALSE)
  }

  lapply(1:nsims, function(i) data[M_samps[[i]], J_samps[[i]]])
}

check_coefs_old <- function(coefs, fit, template=FALSE){
  required_subs <- names(fit@estimates@estimates)
  required_coefs <- lapply(fit@estimates@estimates, function(x) names(x@estimates))
  required_lens <- lapply(required_coefs, length)

  formulas <- sapply(names(fit), function(x) get_formula(fit, x))

  # If there are random effects, adjust the expected coefficient names
  # to remove the b vector and add the grouping covariate name
  rand <- lapply(formulas, lme4::findbars)
  if(!all(sapply(rand, is.null))){
    stopifnot(all(required_subs %in% names(formulas)))
    rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars)))
    if(!all(sapply(rvar, length)<2)){
      stop("Only 1 random effect per parameter is supported", call.=FALSE)
    }
    for (i in required_subs){
      if(!is.null(rand[[i]][[1]])){
        signame <- rvar[[i]]
        old_coefs <- required_coefs[[i]]
        new_coefs <- old_coefs[!grepl("b_", old_coefs, fixed=TRUE)]
        new_coefs <- c(new_coefs, signame)
        required_coefs[[i]] <- new_coefs
      }
    }
  }

  dummy_coefs <- lapply(required_coefs, function(x){
                    out <- rep(0, length(x))
                    x <- gsub("(Intercept)", "intercept", x, fixed=TRUE)
                    names(out) <- x
                    out
                  })

  if(template) return(dummy_coefs)

  if(is.null(coefs)){
    cat("coefs argument should be a named list of named vectors, with the following structure
        (replacing 0s with your desired coefficient values):\n\n")
    print(dummy_coefs)
    stop("Supply coefs argument as specified above", call.=FALSE)
  }

  for (i in 1:length(required_subs)){
    if(!required_subs[i] %in% names(coefs)){
      stop(paste0("Missing required list element '",required_subs[i], "' in coefs list"), call.=FALSE)
    }

    sub_coefs <- coefs[[required_subs[i]]]

    if(is.null(sub_coefs)){
      stop(paste("Required coefficients for the", required_subs[i], "submodel are:",
                  paste(required_coefs[[i]],collapse=", ")))
    }

    is_named <- !is.null(names(sub_coefs)) & !any(names(sub_coefs)=="")

    if(!is_named){
      warning(paste("At least one coefficient in vector for submodel",required_subs[i],
                    "is unnamed; assuming the following order:\n",
                    paste(required_coefs[[i]], collapse=", ")))
      if(length(sub_coefs) != required_lens[i]){
        stop(paste0("Entry '",required_subs[[i]], "' in coefs list must be length ",
        required_lens[[i]]), call.=FALSE)
      }
    } else {
      rsi <- required_subs[i]
      change_int <- names(coefs[[rsi]])%in%c("intercept","Intercept")
      names(coefs[[rsi]])[change_int] <- "(Intercept)"
      change_int <- names(coefs[[rsi]])%in%c("sigmaintercept","sigmaIntercept")
      names(coefs[[rsi]])[change_int] <- "sigma(Intercept)"
      change_int <- names(coefs[[rsi]])%in%c("shapeintercept","shapeIntercept")
      names(coefs[[rsi]])[change_int] <- "shape(Intercept)"
      change_int <- names(coefs[[rsi]])%in%c("rateintercept","rateIntercept")
      names(coefs[[rsi]])[change_int] <- "rate(Intercept)"
      change_int <- grepl(" intercept", names(coefs[[rsi]]))
      names(coefs[[rsi]])[change_int] <- gsub(" intercept", " (Intercept)",
                                              names(coefs[[rsi]])[change_int])
      change_int <- grepl(" Intercept", names(coefs[[rsi]]))
      names(coefs[[rsi]])[change_int] <- gsub(" Intercept", " (Intercept)",
                                              names(coefs[[rsi]])[change_int])
      sub_coefs <- coefs[[rsi]]

      not_inc <- !required_coefs[[i]] %in% names(sub_coefs)
      extra <- !names(sub_coefs) %in% required_coefs[[i]]

      if(any(not_inc)){
        stop(paste("The following required coefficients in the", required_subs[i], "submodel were not found:",
                   paste(required_coefs[[i]][not_inc], collapse=", ")))
      }
      if(any(extra)){
        warning(paste("Ignoring extra coefficients in the", required_subs[i], "submodel:",
                      paste(names(sub_coefs)[extra], collapse=", ")))
      }
      coefs[[rsi]] <- coefs[[rsi]][required_coefs[[i]]]
    }
  }
  coefs[required_subs]
}

setMethod("summary", "unmarkedPower_old", function(object, ...){
  sum_dfs <- object@estimates
  npar <- nrow(sum_dfs[[1]])

  nulls <- object@nulls
  nulls <- lapply(nulls, function(x){
    nm <- names(x)
    nm[nm %in% c("Intercept","intercept")] <- "(Intercept)"
    names(x) <- nm
    x
  })

  coefs_no_rand <- unlist(object@coefs)[!grepl("b_", names(unlist(object@coefs)))]

  pow <- sapply(1:npar, function(ind){
    submod <- sum_dfs[[1]]$submodel[ind]
    param <- sum_dfs[[1]]$param[ind]
    ni <- nulls[[submod]][param]

    pcrit <- sapply(sum_dfs, function(x) wald(x$Estimate[ind], x$SE[ind], ni)) < object@alpha
    direct <- sapply(sum_dfs, function(x) diff_dir(x$Estimate[ind], coefs_no_rand[ind], ni))
    mean(pcrit & direct, na.rm=T)
  })

  all_nulls <- sapply(1:npar, function(ind){
    submod <- sum_dfs[[1]]$submodel[ind]
    param <- sum_dfs[[1]]$param[ind]
    ni <- nulls[[submod]][param]
    if(is.null(ni) || is.na(ni)) ni <- 0
    ni
  })

  effect_no_random <- unlist(object@coefs)[!grepl("b_",names(unlist(object@coefs)))]

  out <- cbind(sum_dfs[[1]][,1:2], effect=effect_no_random, null=all_nulls,  power=pow)
  rownames(out) <- NULL
  names(out) <- c("Submodel", "Parameter", "Effect", "Null", "Power")
  out
})

setMethod("show", "unmarkedPower_old", function(object){
  cat("\nModel:\n")
  print(object@call)
  cat("\n")

  cat("Power Statistics:\n")
  sumtab <- summary(object)
  sumtab$Power <- round(sumtab$Power, 3)
  print(sumtab, row.names=FALSE)
})



get_summary_df_old <- function(fit){
  n_est <- length(fit@estimates@estimates)
  #est_names <- unname(sapply(fit@estimates@estimates, function(x) x@name))
  est_names <- names(fit@estimates@estimates)
  all_est <- lapply(1:n_est, function(i){
    utils::capture.output(out <- summary(fit@estimates@estimates[[i]]))
    out <- cbind(submodel=est_names[i], param=rownames(out), out)
    rownames(out) <- NULL
    out
  })
  do.call(rbind, all_est)
}

setClass("unmarkedPowerList_old", representation(powerAnalyses="list"))

setMethod("unmarkedPowerList", "list", function(object, ...){
  new("unmarkedPowerList_old", powerAnalyses=object)
})

setMethod("unmarkedPowerList", "unmarkedFit",
  function(object, coefs, design, alpha=0.05, nulls=list(),
           nsim=100, parallel=FALSE, ...){

  ndesigns <- nrow(design)
  out <- lapply(1:ndesigns, function(i){
    cat(paste0("M = ",design$M[i],", J = ",design$J[i],"\n"))
    powerAnalysis(object, coefs, as.list(design[i,]), alpha=alpha, nsim=nsim,
                  nulls=nulls, parallel=FALSE)
  })
  new("unmarkedPowerList_old", powerAnalyses=out)
})

setMethod("summary", "unmarkedPowerList_old", function(object, ...){
  out <- lapply(object@powerAnalyses, function(x){
    stats <- summary(x)
    cbind(M=x@M, T=x@T, J=x@J, stats)
  })
  out <- do.call(rbind, out)
  out$M <- factor(out$M)
  out$T <- factor(out$T)
  out$J <- factor(out$J)
  out
})

setMethod("show", "unmarkedPowerList_old", function(object){
  print(summary(object))
})

setMethod("plot", "unmarkedPowerList_old", function(x, power=NULL, param=NULL, ...){
  dat <- summary(x)
  if(is.null(param)) param <- dat$Parameter[1]
  dat <- dat[dat$Parameter==param,,drop=FALSE]
  ylim <- range(dat$Power, na.rm=T)
  if(!is.null(power)) ylim[2] <- max(power, ylim[2])
  xlim <- range(as.numeric(as.character(dat$M)), na.rm=T)
  cols <- palette.colors(length(levels(dat$J)), palette="Dark 2")
  old_par <- graphics::par()[c("mfrow","mar")]
  nT <- length(levels(dat$T))
  mar <- old_par$mar
  if(nT == 1) mar <- c(5.1, 4.1, 2.1, 2.1)
  graphics::par(mfrow=c(length(levels(dat$T)),1), mar=mar)
  for (i in levels(dat$T)){
    plot_title <- ""
    if(nT > 1) plot_title <- paste0("T = ", i)
    tsub <- dat[dat$T==i,,drop=FALSE]
    Jlev <- levels(tsub$J)
    jsub <- tsub[tsub$J==Jlev[1],,drop=FALSE]
    plot(as.numeric(as.character(jsub$M)), jsub$Power, type="o",
        col=cols[1], ylim=ylim, xlim=xlim, xlab="Sites",
        ylab="Power", pch=19, main=plot_title)
    if(!is.null(power)) abline(h=power, lty=2)
    for (j in 2:length(Jlev)){
      jsub <- tsub[tsub$J==Jlev[j],,drop=FALSE]
      graphics::lines(as.numeric(as.character(jsub$M)), jsub$Power, type="o",
            col=cols[j], pch=19)
    }
    graphics::legend('bottomright', lwd=1, pch=19, col=cols, legend=Jlev, title="Observations")
  }
  graphics::par(mfrow=old_par)
})

setMethod("update", "unmarkedPower_old", function(object, ...){
  args <- list(...)
  if(!is.null(args$alpha)) object@alpha <- args$alpha
  if(!is.null(args$coefs)){
    if(!is.list(args$coefs) || all(names(args$coefs) == names(object@coefs))){
      stop("coefs list structure is incorrect", call.=FALSE)
      object@coefs <- args$coefs
    }
  }
  if(!is.null(args$nulls)) object@nulls <- args$nulls
  object
})

shinyPower <- function(object, ...){

  if(!inherits(object, "unmarkedFit")){
    stop("Requires unmarkedFit object", call.=FALSE)
  }

  .Deprecated("shinyPower", package=NULL, 
              msg = paste("shinyPower used on an unmarkedFit object will soon be removed from unmarked."),
             old = as.character(sys.call(sys.parent()))[1L])
  if(!requireNamespace("shiny")){
    stop("Install the shiny library to use this function", call.=FALSE)
  }
  if(!requireNamespace("pbapply")){
    stop("Install the pbapply library to use this function", call.=FALSE)
  }
  options(unmarked_shiny=TRUE)
  on.exit(options(unmarked_shiny=FALSE))
  .shiny_env$.SHINY_MODEL <- object

  shiny::runApp(system.file("shinyPower", package="unmarked"))

}