Standardize getDesign inputs

R/All_classes.R

@@ -275,7 +275,7 @@ setClass("unmarkedEstimateList",
 setClass("unmarkedFit",
   representation(fitType = "character",
     call = "call",
-    formula = "formula",
+    formlist = "list",
     data = "unmarkedFrame",
     sitesRemoved = "numeric",  # vector of indices of removed sites
     estimates = "unmarkedEstimateList",
@@ -292,22 +292,25 @@ setClass("unmarkedFit",
 
 # Distance sampling
 setClass("unmarkedFitDS",
-  representation(keyfun = "character", unitsOut = "character", output = "character"), 
+  representation(formula = "formula", keyfun = "character", unitsOut = "character", output = "character"), 
   contains = "unmarkedFit"
 )
 
 # Multinomial mixture model
-setClass("unmarkedFitMPois", contains = "unmarkedFit")
+setClass("unmarkedFitMPois", 
+  representation(formula = "formula"),
+  contains = "unmarkedFit"
+)
 
 # N-mixture TTD model
 setClass("unmarkedFitNmixTTD",
-  representation(stateformula = "formula", detformula = "formula", K = "numeric"),
+  representation(K = "numeric"),
   contains = "unmarkedFit"
 )
 
 # Basic occupancy model
 setClass("unmarkedFitOccu",
-  representation(knownOcc = "logical"),
+  representation(formula = "formula", knownOcc = "logical"),
   contains = "unmarkedFit"
 )
 
@@ -316,27 +319,26 @@ setClass("unmarkedFitOccuComm", contains="unmarkedFitOccu")
 
 # Count-based occupancy model
 setClass("unmarkedFitOccuCOP",
-  representation(removed_obs = "matrix", formlist = "list"),
+  representation(removed_obs = "matrix"),
   contains = "unmarkedFit"
 )
 
 # False positive occupancy
 setClass("unmarkedFitOccuFP",
-  representation(knownOcc = "logical", detformula = "formula", FPformula = "formula",
-                 Bformula = "formula", stateformula = "formula", type = "numeric"),
+  representation(knownOcc = "logical", type = "numeric"),
   contains = "unmarkedFit"
 )
 
 # Penalized likelihood occupancy model
 setClass("unmarkedFitOccuPEN",
-  representation(knownOcc = "logical", pen.type = "character", lambda = "numeric"),
+  representation(formula = "formula", knownOcc = "logical", pen.type = "character", lambda = "numeric"),
   contains = "unmarkedFit"
 )
 
 setClass("unmarkedFitOccuPEN_CV",
-  representation(knownOcc = "logical", pen.type = "character", lambdaVec = "numeric",
-	               k = "numeric", foldAssignments = "numeric", lambdaScores = "numeric", 
-                 chosenLambda = "numeric"),
+  representation(formula = "formula", knownOcc = "logical", pen.type = "character", 
+                 lambdaVec = "numeric", k = "numeric", foldAssignments = "numeric", 
+                 lambdaScores = "numeric", chosenLambda = "numeric"),
   contains = "unmarkedFit"
 )
 
@@ -348,28 +350,25 @@ setClass("unmarkedFitOccuMulti",
 
 # Royle-Nichols occupancy model
 setClass("unmarkedFitOccuRN",
-  representation(K = "numeric"),
+  representation(formula = "formula", K = "numeric"),
   contains = "unmarkedFit"
 )
 
 # N-mixture model
 setClass("unmarkedFitPCount",
-  representation(K = "numeric", mixture = "character"),
+  representation(formula = "formula", K = "numeric", mixture = "character"),
   contains = "unmarkedFit"
 )
 
 
 ### Temporary emigration (TE) model types
 
 # TE occupancy
-setClass("unmarkedFitGOccu",
-  representation(formlist = "list"),
-  contains = "unmarkedFit"
-)
+setClass("unmarkedFitGOccu", contains = "unmarkedFit")
 
 # TE multinomial mixture
 setClass("unmarkedFitGMM",
-  representation(formlist = "list", mixture = "character", K = "numeric"),
+  representation(mixture = "character", K = "numeric"),
   contains = "unmarkedFit"
 )
 
@@ -390,8 +389,7 @@ setClass("unmarkedFitGPC", contains = "unmarkedFitGMM")
 
 # Dynamic occupancy model
 setClass("unmarkedFitColExt",
-  representation(phi = "matrix", psiformula = "formula", gamformula = "formula",
-                 epsformula = "formula", detformula = "formula", projected = "array",
+  representation(phi = "matrix", projected = "array",
                  projected.mean = "matrix", smoothed = "array", smoothed.mean = "matrix",
                  projected.mean.bsse = "optionalMatrix", smoothed.mean.bsse = "optionalMatrix"),
   contains = "unmarkedFit"
@@ -405,16 +403,12 @@ setClass("unmarkedFitOccuMS",
 )
 
 # Time-to-detection occupancy model
-setClass("unmarkedFitOccuTTD",
-  representation(psiformula = "formula", gamformula = "formula",
-                 epsformula = "formula", detformula = "formula"),
-  contains = "unmarkedFit"
-)
+setClass("unmarkedFitOccuTTD", contains = "unmarkedFit")
 
 # This class is not used directly, just used as a base for for PCO, MMO, DSO
 setClass("unmarkedFitDailMadsen",
-  representation(K = "numeric", mixture = "character", formlist = "list",
-                 dynamics = "character", immigration = "logical", fix = "character"),
+  representation(K = "numeric", mixture = "character", dynamics = "character",
+                 immigration = "logical", fix = "character"),
   contains = "unmarkedFit"
 )
 
@@ -436,7 +430,6 @@ setClass("unmarkedFitPCO", contains = "unmarkedFitDailMadsen")
 setClassUnion("unmarkedFrameOrNULL", members=c("unmarkedFrame", "NULL"))
 setClass("unmarkedFitIDS",
     representation(
-        formlist = "list",
         keyfun = "character",
         K = "numeric",
         dataPC = "unmarkedFrameOrNULL",

R/IDS.R

@@ -96,23 +96,31 @@ IDS <- function(lambdaformula = ~1,
   # Design matrices------------------------------------------------------------
 
   # Need to add offset support here eventually
-  gd_hds <- getDesign(dataDS, form_hds)
+  form_hds_split <- split_formula(form_hds)
+  names(form_hds_split) <- c("det", "state")
+  gd_hds <- getDesign(dataDS, form_hds_split)
   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_state
+  if(has_avail){
+    form_avail_split <- split_formula(form_avail)
+    names(form_avail_split) <- c("det", "state")
+    Xavail_ds <- getDesign(dataDS, form_avail_split)$X_state
+  }
   if(is.null(durationDS)) durationDS <- rep(0,0)
 
   gd_pc <- list(y=matrix(0,0,0), X_state=matrix(0,0,0), X_det=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)
+    form_pc_split <- split_formula(form_pc)
+    names(form_pc_split) <- c("det", "state")
+    gd_pc <- getDesign(dataPC, form_pc_split)
     ds_pc <- get_ds_info(c(0, maxDistPC))
-    if(has_avail) Xavail_pc <- getDesign(dataPC, form_avail)$X_state
+    if(has_avail) Xavail_pc <- getDesign(dataPC, form_avail_split)$X_state
   }
 
   gd_oc <- list(y=matrix(0,0,0), X_state=matrix(0,0,0), X_det=matrix(0,0,0))
@@ -121,10 +129,12 @@ IDS <- function(lambdaformula = ~1,
   Xavail_oc <- matrix(0,0,0)
   if(is.null(durationOC)) durationOC <- rep(0,0)
   if(!is.null(dataOC)){
-    gd_oc <- getDesign(dataOC, form_oc)
+    form_oc_split <- split_formula(form_pc)
+    names(form_oc_split) <- c("det", "state")
+    gd_oc <- getDesign(dataOC, form_oc_split)
     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_state
+    if(has_avail) Xavail_oc <- getDesign(dataOC, form_avail_split)$X_state
   }
 
   # Density conversion and unequal area correction
@@ -322,7 +332,7 @@ IDS <- function(lambdaformula = ~1,
   est_list <- unmarkedEstimateList(est_list)
 
   new("unmarkedFitIDS", fitType = "IDS", call = match.call(),
-    opt = opt, formula = lambdaformula, formlist=formlist,
+    opt = opt, formlist=formlist,
     data = dataDS, dataPC=dataPC, dataOC=dataOC, K=K,
     surveyDurations=surveyDurations,
     maxDist = list(pc=maxDistPC, oc=maxDistOC),
@@ -400,8 +410,10 @@ IDS_convert_class <- function(inp, type, ds_type=NULL){
 
   form <- inp@formlist[[type]]
   if(type=="phi") form <- as.formula(paste(c(as.character(form), "~1"), collapse=""))
+  formlist <- split_formula(form)
+  names(formlist) <- c("det", "state")
 
-  new("unmarkedFitDS", fitType="IDS", opt=inp@opt, formula=form,
+  new("unmarkedFitDS", fitType="IDS", opt=inp@opt, formula=form, formlist=formlist,
       data=data, keyfun=inp@keyfun, unitsOut=inp@unitsOut,
       estimates=unmarkedEstimateList(est),
       AIC=inp@AIC, output="density", TMB=inp@TMB)

R/colext.R

@@ -5,16 +5,17 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1,
   ## truncate to 1
   data@y <- truncateToBinary(data@y)
 
-  formula <- list(psiformula = psiformula, gammaformula = gammaformula,
-                  epsilonformula = epsilonformula, pformula = pformula)
-  check_no_support(formula)
-  dm <- getDesign(data, formula = as.formula(paste(unlist(formula), collapse=" ")))
-  X_psi <- dm$X_state; X_col <- dm$X_col; X_ext <- dm$X_ext
+  formulas <- list(psi = psiformula, col = gammaformula, ext = epsilonformula,
+                   det = pformula)
+  check_no_support(formulas)
+  dm <- getDesign(data, formulas)
+
+  X_col <- dm$X_col; X_ext <- dm$X_ext
   y <- dm$y
   M <- nrow(y)
   T <- data@numPrimary
   J <- ncol(y) / T
-  psiParms <- colnames(X_psi)
+  psiParms <- colnames(dm$X_psi)
   colParms <- colnames(X_col)
   extParms <- colnames(X_ext)
   detParms <- colnames(dm$X_det)
@@ -56,7 +57,7 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1,
   pind_mat[3,] <- max(pind_mat[2,]) + c(1, length(extParms))
   pind_mat[4,] <- max(pind_mat[3,]) + c(1, length(detParms))
 
-  tmb_dat <- list(y = as.vector(t(y)), X_psi = X_psi, X_col = X_col, 
+  tmb_dat <- list(y = as.vector(t(y)), X_psi = dm$X_psi, X_col = X_col, 
                   X_ext = X_ext, X_det = dm$X_det,
                   M = M, T = T, J = J, 
                   site_sampled = site_sampled, nd = no_detects)
@@ -109,7 +110,7 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1,
   estimateList <- unmarkedEstimateList(list(psi = psi, col = col,
                                             ext = ext, det = det))
 
-  psis <- plogis(X_psi %*% psi_coef$ests)
+  psis <- plogis(dm$X_psi %*% psi_coef$ests)
 
   # Compute projected estimates
   phis <- array(NA,c(2,2,T-1,M))
@@ -140,11 +141,7 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1,
 
   umfit <- new("unmarkedFitColExt", fitType = "colext",
                 call = match.call(),
-                formula = as.formula(paste(unlist(formula),collapse=" ")),
-                psiformula = psiformula,
-                gamformula = gammaformula,
-                epsformula = epsilonformula,
-                detformula = pformula,
+                formlist = formulas,
                 data = data, sitesRemoved = dm$removed.sites,
                 estimates = estimateList,
                 AIC = fmAIC, opt = opt, negLogLike = opt$value,

R/distsamp.R

@@ -8,15 +8,17 @@ distsamp <- function(formula, data,
 
     # Check arguments
     engine <- match.arg(engine)
-    if(any(sapply(split_formula(formula), has_random))) engine <- "TMB"
+    formulas <- split_formula(formula)
+    names(formulas) <- c("det", "state")
+    if(any(sapply(formulas, has_random))) engine <- "TMB"
     keyfun <- match.arg(keyfun)
     output <- match.arg(output)
     unitsOut <- match.arg(unitsOut)
 
     if(missing(starts)) starts <- NULL
 
     #Generate design matrix
-    dm <- getDesign(data, formula)
+    dm <- getDesign(data, formulas)
     y <- dm$y
     M <- nrow(y)
     J <- ncol(y)
@@ -253,8 +255,7 @@ distsamp <- function(formula, data,
 
       # Set up TMB input data
       if(output == "abund") A <- rep(1, length(A))
-      forms <- split_formula(formula)
-      inps <- get_ranef_inputs(forms, list(det=siteCovs(data), state=siteCovs(data)),
+      inps <- get_ranef_inputs(formulas, list(det=siteCovs(data), state=siteCovs(data)),
                                list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")])
 
       keyfun_type <- switch(keyfun, uniform={0}, halfnorm={1}, exp={2},
@@ -290,8 +291,8 @@ distsamp <- function(formula, data,
       }
 
       # Organize random-effect estimates from TMB output
-      state_rand_info <- get_randvar_info(tmb_out$sdr, "state", forms[[2]], siteCovs(data))
-      det_rand_info <- get_randvar_info(tmb_out$sdr, "det", forms[[1]], siteCovs(data))
+      state_rand_info <- get_randvar_info(tmb_out$sdr, "state", formulas$state, siteCovs(data))
+      det_rand_info <- get_randvar_info(tmb_out$sdr, "det", formulas$det, siteCovs(data))
 
     }
 
@@ -317,7 +318,7 @@ distsamp <- function(formula, data,
     }
 
     dsfit <- new("unmarkedFitDS", fitType = "distsamp", call = match.call(),
-        opt = fm, formula = formula, data = data, keyfun=keyfun,
+        opt = fm, formula = formula, formlist = formulas, data = data, keyfun=keyfun,
         sitesRemoved = dm$removed.sites, unitsOut=unitsOut,
         estimates = estimateList, AIC = fmAIC, negLogLike = fm$value,
         nllFun = nll, output=output, TMB=tmb_mod)

R/distsampOpen.R

@@ -44,12 +44,11 @@ distsampOpen <- function(lambdaformula, gammaformula, omegaformula, pformula,
 
   fix <- match.arg(fix)
 
-  formlist <- mget(c("lambdaformula", "gammaformula", "omegaformula",
-                   "pformula", "iotaformula"))
-  check_no_support(formlist)
-  formula <- as.formula(paste(unlist(formlist), collapse=" "))
+  formulas <- list(lambda = lambdaformula, gamma = gammaformula,
+                   omega = omegaformula, det = pformula, iota = iotaformula)
+  check_no_support(formulas)
 
-  D <- getDesign(data, formula)
+  D <- getDesign(data, formulas)
   y <- D$y
 
   deltamax <- max(D$delta, na.rm=TRUE)
@@ -307,7 +306,7 @@ distsampOpen <- function(lambdaformula, gammaformula, omegaformula, pformula,
   }
 
   umfit <- new("unmarkedFitDSO", fitType = "distsampOpen",
-      call = match.call(), formula = formula, formlist = formlist, data = data,
+      call = match.call(), formlist = formulas, data = data,
       sitesRemoved=D$removed.sites, estimates = estimateList, AIC = fmAIC,
       opt = fm, negLogLike = fm$value, nllFun = nll, K = K, mixture = mixture,
       dynamics = dynamics, fix = fix, immigration=immigration, keyfun=keyfun,

R/fitted.R

@@ -319,7 +319,7 @@ setMethod("fitted_internal", "unmarkedFitDailMadsen", function(object){
     fix <- tryCatch(object@fix, error=function(e) "none")
     immigration <- tryCatch(object@immigration, error=function(e) FALSE)
     data <- getData(object)
-    D <- getDesign(data, object@formula, na.rm = FALSE)
+    D <- getDesign(data, object@formlist, na.rm = FALSE)
     delta <- D$delta #FIXME this isn't returned propertly when na.rm=F
 
     M <- numSites(data)