diff --git a/R/All_classes.R b/R/All_classes.R
index dae86cde..544bbf01 100644
--- a/R/All_classes.R
+++ b/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,13 +350,13 @@ 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"
 )
 
@@ -362,14 +364,11 @@ setClass("unmarkedFitPCount",
 ### 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",
diff --git a/R/IDS.R b/R/IDS.R
index 3ff994cc..069000c9 100644
--- a/R/IDS.R
+++ b/R/IDS.R
@@ -96,13 +96,19 @@ 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))
@@ -110,9 +116,11 @@ IDS <- function(lambdaformula = ~1,
   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)
diff --git a/R/colext.R b/R/colext.R
index 4560154e..d089ddc8 100644
--- a/R/colext.R
+++ b/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,
diff --git a/R/distsamp.R b/R/distsamp.R
index 6de81329..d0db5fa6 100644
--- a/R/distsamp.R
+++ b/R/distsamp.R
@@ -8,7 +8,9 @@ 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)
@@ -16,7 +18,7 @@ distsamp <- function(formula, data,
     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)
diff --git a/R/distsampOpen.R b/R/distsampOpen.R
index d253d810..39e36051 100644
--- a/R/distsampOpen.R
+++ b/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,
diff --git a/R/fitted.R b/R/fitted.R
index b485c079..a6db1b0e 100644
--- a/R/fitted.R
+++ b/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)
diff --git a/R/gdistremoval.R b/R/gdistremoval.R
index 3bb91a8d..8ec2213c 100644
--- a/R/gdistremoval.R
+++ b/R/gdistremoval.R
@@ -65,8 +65,9 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1,
   mixture <- match.arg(mixture)
   engine <- match.arg(engine)
 
-  formlist <- mget(c("lambdaformula", "phiformula", "distanceformula", "removalformula"))
-  if(any(sapply(formlist, has_random))) engine <- "TMB"
+  formulas <- list(lambda = lambdaformula, phi = phiformula, dist = distanceformula,
+                   rem = removalformula)
+  if(any(sapply(formulas, has_random))) engine <- "TMB"
 
   M <- numSites(data)
   T <- data@numPrimary
@@ -74,7 +75,7 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1,
   Rrem <- ncol(data@yRemoval)
   mixture_code <- switch(mixture, P={1}, NB={2}, ZIP={3})
 
-  gd <- getDesign(data, formlist)
+  gd <- getDesign(data, formulas)
 
   Jdist <- Rdist / T
   ysum <- array(t(gd$yDist), c(Jdist, T, M))
@@ -174,7 +175,7 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1,
 
     dlist <- list(lambda=siteCovs(data), phi=yearlySiteCovs(data),
                   dist=siteCovs(data), rem=obsCovs(data))
-    inps <- get_ranef_inputs(formlist, dlist,
+    inps <- get_ranef_inputs(formulas, dlist,
                              gd[c("X_lambda","X_phi","X_dist","X_rem")],
                              gd[c("Z_lambda","Z_phi","Z_dist","Z_rem")])
 
@@ -264,8 +265,8 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1,
       randomVarInfo=rem_rand_info)
 
   new("unmarkedFitGDR", fitType = "gdistremoval",
-    call = match.call(), formula = as.formula(paste(formlist, collapse="")),
-    formlist = formlist, data = data, estimates = estimateList, sitesRemoved = numeric(0),
+    call = match.call(), formlist = formulas, data = data, 
+    estimates = estimateList, sitesRemoved = numeric(0),
     AIC = fmAIC, opt = opt, negLogLike = opt$value, nllFun = nll,
     mixture=mixture, K=K, keyfun=keyfun, unitsOut=unitsOut, output=output, TMB=tmb_mod)
 
@@ -572,10 +573,10 @@ setMethod("y_to_zeros", "unmarkedFrameGDR", function(object, ...){
 setMethod("rebuild_call", "unmarkedFitGDR", function(object){           
   cl <- object@call
   cl[["data"]] <- quote(object@data)
-  cl[["lambdaformula"]] <- object@formlist$lambdaformula
-  cl[["phiformula"]] <- object@formlist$phiformula
-  cl[["removalformula"]] <- object@formlist$removalformula
-  cl[["distanceformula"]] <- object@formlist$distanceformula
+  cl[["lambdaformula"]] <- object@formlist$lambda
+  cl[["phiformula"]] <- object@formlist$phi
+  cl[["removalformula"]] <- object@formlist$rem
+  cl[["distanceformula"]] <- object@formlist$dist
   cl[["mixture"]] <- object@mixture
   cl[["K"]] <- object@K
   cl[["keyfun"]] <- object@keyfun
diff --git a/R/gdistsamp.R b/R/gdistsamp.R
index 9f7ca562..34c679b4 100644
--- a/R/gdistsamp.R
+++ b/R/gdistsamp.R
@@ -21,14 +21,9 @@ survey <- data@survey
 unitsIn <- data@unitsIn
 mixture <- match.arg(mixture)
 
-formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula,
-    pformula = pformula)
-check_no_support(formlist)
-form <- as.formula(paste(unlist(formlist), collapse=" "))
-D <- getDesign(data, formula = form)
-X_lambda <- D$X_state
-offset_lambda <- D$offset_state
-
+formulas <- list(lambda = lambdaformula, phi = phiformula, det = pformula)
+check_no_support(formulas)
+D <- getDesign(data, formulas)
 y <- D$y  # MxJT
 
 M <- nrow(y)
@@ -82,7 +77,7 @@ switch(unitsOut,
     kmsq = A <- A)
 
 
-lamPars <- colnames(X_lambda)
+lamPars <- colnames(D$X_lambda)
 if(T==1) {
     phiPars <- character(0)
     nPP <- 0
@@ -118,7 +113,7 @@ if(identical(mixture, "NB")) {
     nbPar <- character(0)
 }
 
-nLP <- ncol(X_lambda)
+nLP <- ncol(D$X_lambda)
 nP  <- nLP + nPP + nDP + nSP + nOP
 
 cp <- array(as.numeric(NA), c(M, T, J+1))
@@ -149,7 +144,7 @@ halfnorm = {
         }
 
     nll_R <- function(pars) {
-        lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+        lambda <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
         if(identical(output, "density"))
             lambda <- lambda * A
 
@@ -219,7 +214,7 @@ exp = {
         }
 
     nll_R <- function(pars) {
-        lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+        lambda <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
         if(identical(output, "density"))
             lambda <- lambda * A
         if(T==1)
@@ -287,7 +282,7 @@ hazard = {
         starts <- rep(0, nP)
         }
     nll_R <- function(pars) {
-        lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+        lambda <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
         if(identical(output, "density"))
             lambda <- lambda * A
 
@@ -356,7 +351,7 @@ uniform = {
         starts <- rep(0, nP)
         }
     nll_R <- function(pars) {
-        lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+        lambda <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
         if(identical(output, "density"))
             lambda <- lambda * A
         if(T==1)
@@ -410,7 +405,7 @@ if(engine =="C"){
 
   nll <- function(params){
     nll_gdistsamp(params, n_param, y_long, mixture_code, keyfun, survey,
-                  X_lambda, offset_lambda, A, D$X_phi, D$offset_phi, D$X_det, D$offset_det,
+                  D$X_lambda, D$offset_lambda, A, D$X_phi, D$offset_phi, D$X_det, D$offset_det,
                   db, a, t(u), w, k, lfac.k, lfac.kmytC, kmytC, Kmin, threads)
   }
 
@@ -467,7 +462,7 @@ if(identical(mixture,"ZIP")) {
 }
 
 umfit <- new("unmarkedFitGDS", fitType = "gdistsamp",
-    call = match.call(), formula = form, formlist = formlist,
+    call = match.call(), formlist = formulas,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
     AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll,
     mixture=mixture, K=K, keyfun=keyfun, unitsOut=unitsOut, output=output)
diff --git a/R/getDesign.R b/R/getDesign.R
index 0d875312..269fd521 100644
--- a/R/getDesign.R
+++ b/R/getDesign.R
@@ -5,28 +5,24 @@
 # generic method for unmarkedFrame
 # used by distsamp, multinomPois, occu, occuPEN, occuRN, pcount, pcount.spHDS, IDS
 setMethod("getDesign", "unmarkedFrame", 
-  function(umf, formula, na.rm = TRUE, ...){
+  function(umf, formulas, na.rm = TRUE, ...){
 
   M <- numSites(umf)
   J <- obsNum(umf)
   y <- getY(umf)
 
-  formulas <- split_formula(formula)
-  stateformula <- formulas[[2]]
-  detformula <- formulas[[1]]
-
   # Process covariates
   covs <- clean_up_covs(umf)
   
   # Model matrices and offset for state submodel
-  X_state <- get_model_matrix(stateformula, covs$site_covs)
-  offset_state <- get_offset(stateformula, covs$site_covs)
-  Z_state <- get_Z(stateformula, covs$site_covs)
+  X_state <- get_model_matrix(formulas$state, covs$site_covs)
+  offset_state <- get_offset(formulas$state, covs$site_covs)
+  Z_state <- get_Z(formulas$state, covs$site_covs)
 
   # Model matrices and offset for detection submodel
-  X_det <- get_model_matrix(detformula, covs$obs_covs)
-  offset_det <- get_offset(detformula, covs$obs_covs)
-  Z_det <- get_Z(detformula, covs$obs_covs)
+  X_det <- get_model_matrix(formulas$det, covs$obs_covs)
+  offset_det <- get_offset(formulas$det, covs$obs_covs)
+  Z_det <- get_Z(formulas$det, covs$obs_covs)
 
   # Identify missing values in state covs
   has_na_site <- row_has_na(cbind(X_state, Z_state))
@@ -37,7 +33,7 @@ setMethod("getDesign", "unmarkedFrame",
   has_na_obs <- matrix(has_na_obs, M, J, byrow=TRUE)
   # Multiplying by obsToY handles models where the number of estimated parameters
   # is not the same as the number of observations, such as double-observer models
-  has_na_obs <- has_na_obs %*% umf@obsToY > 0
+  has_na_obs <- has_na_obs %*% obsToY(umf) > 0
 
   # Combine missing value information
   has_na <- has_na_site | has_na_obs
@@ -70,17 +66,13 @@ setMethod("getDesign", "unmarkedFrame",
 # UnmarkedMultFrame
 # used by colext, occuTTD, nmixTTD and base class for G3 and DailMadsen
 setMethod("getDesign", "unmarkedMultFrame",
-  function(umf, formula, na.rm = TRUE){
+  function(umf, formulas, na.rm = TRUE){
 
-  stateformula <- as.formula(formula[[2]][[2]][[2]])  
-  gamformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[2]][[3]])))
-  epsformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[3]])))
-  detformula <- as.formula(as.call(list(as.name("~"), formula[[3]])))
- 
-  # Process state and detection with generic umf method  
-  comb_form <- list(as.name("~"), detformula, stateformula[[2]])
-  comb_form <- as.formula(as.call(comb_form))
-  out <- methods::callNextMethod(umf, formula = comb_form, na.rm = FALSE)
+  # Handle varying name of state variable in models that use this method
+  state_type <- names(formulas)[1]
+  # Process state and detection with generic umf method
+  names(formulas)[1] <- "state"
+  out <- methods::callNextMethod(umf, formulas = formulas, na.rm = FALSE)
  
   M <- numSites(umf)
   R <- obsNum(umf)
@@ -95,12 +87,12 @@ setMethod("getDesign", "unmarkedMultFrame",
   covs <- clean_up_covs(umf, drop_final = TRUE)
 
   # Model matrix for colonization
-  X_col <- get_model_matrix(gamformula, covs$yearly_site_covs)
-  offset_col <- get_offset(gamformula, covs$yearly_site_covs)
+  X_col <- get_model_matrix(formulas$col, covs$yearly_site_covs)
+  offset_col <- get_offset(formulas$col, covs$yearly_site_covs)
 
   # Model matrix for extinction
-  X_ext <- get_model_matrix(epsformula, covs$yearly_site_covs)
-  offset_ext <- get_offset(epsformula, covs$yearly_site_covs)
+  X_ext <- get_model_matrix(formulas$ext, covs$yearly_site_covs)
+  offset_ext <- get_offset(formulas$ext, covs$yearly_site_covs)
 
   # Error if any offsets specified
   if(any(c(offset_col, offset_ext, out$X.offset, out$V.offset) != 0)){
@@ -132,25 +124,26 @@ setMethod("getDesign", "unmarkedMultFrame",
   }
 
   # Combine outputs
-  list(y = y, X_state = out$X_state, X_col = X_col, X_ext = X_ext, X_det = out$X_det,
-       removed.sites = which(drop_sites))
+  out <- list(y = y, X_state = out$X_state, X_col = X_col, X_ext = X_ext, X_det = out$X_det,
+              removed.sites = which(drop_sites))
+  # Rename design matrix for state variable to correct type
+  names(out)[names(out) == "X_state"] <- paste0("X_", state_type)
+  out
 })
 
 
 # unmarkedFrameG3 (gpcount, gmultmix, gdistsamp, goccu)
 setMethod("getDesign", "unmarkedFrameG3",
-  function(umf, formula, na.rm = TRUE){
+  function(umf, formulas, na.rm = TRUE){
 
-  stateformula <- as.formula(formula[[2]][[2]])  
-  phiformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[3]])))
-  detformula <- as.formula(as.call(list(as.name("~"), formula[[3]])))
- 
+  # Handle varying name of state variable in models that use this method
+  state_type <- names(formulas)[1]
+  
   # Process state and detection with generic umf method  
-  comb_form <- list(as.name("~"), detformula, stateformula[[2]])
-  comb_form <- as.formula(as.call(comb_form))
   # Have to use getMethod because this inherits from unmarkedMultFrame
+  names(formulas)[1] <- "state"
   getDesign_generic <- methods::getMethod("getDesign", "unmarkedFrame")
-  out <- getDesign_generic(umf, formula = comb_form, na.rm = FALSE)
+  out <- getDesign_generic(umf, formulas = formulas, na.rm = FALSE)
  
   M <- numSites(umf)
   R <- obsNum(umf)
@@ -163,8 +156,8 @@ setMethod("getDesign", "unmarkedFrameG3",
   covs <- clean_up_covs(umf)
 
   # Model matrix for availability
-  X_phi <- get_model_matrix(phiformula, covs$yearly_site_covs)
-  offset_phi <- get_offset(phiformula, covs$yearly_site_covs)
+  X_phi <- get_model_matrix(formulas$phi, covs$yearly_site_covs)
+  offset_phi <- get_offset(formulas$phi, covs$yearly_site_covs)
 
   # Check missing values in availability
   has_na <- row_has_na(X_phi)
@@ -190,22 +183,20 @@ setMethod("getDesign", "unmarkedFrameG3",
   }
 
   # Combine outputs
-  list(y = y, X_state = out$X_state, offset_state = out$offset_state,
-       X_phi = X_phi, offset_phi = offset_phi,
-       X_det = out$X_det, offset_det = out$offset_det,
-       removed.sites = which(drop_sites))
+  out <- list(y = y, X_state = out$X_state, offset_state = out$offset_state,
+              X_phi = X_phi, offset_phi = offset_phi,
+              X_det = out$X_det, offset_det = out$offset_det,
+              removed.sites = which(drop_sites))
+  # Rename design matrix and offset for state variable to correct type
+  names(out)[names(out) == "X_state"] <- paste0("X_", state_type)
+  names(out)[names(out) == "offset_state"] <- paste0("offset_", state_type)
+  out
 })
 
 
 # unmarkedFrameDailMadsen (pcountOpen, multmixOpen, distsampOpen)
 setMethod("getDesign", "unmarkedFrameDailMadsen",
-  function(umf, formula, na.rm = TRUE){
-
-  lamformula <- as.formula(formula[[2]][[2]][[2]][[2]])  
-  gamformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[2]][[2]][[3]])))
-  omformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[2]][[3]])))
-  pformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[3]])))
-  iotaformula <- as.formula(as.call(list(as.name("~"), formula[[3]])))
+  function(umf, formulas, na.rm = TRUE){
 
   M <- numSites(umf)
   T <- umf@numPrimary
@@ -219,22 +210,22 @@ setMethod("getDesign", "unmarkedFrameDailMadsen",
   covs <- clean_up_covs(umf, drop_final = TRUE)
 
   # Model matrix for abundance
-  X_lambda <- get_model_matrix(lamformula, covs$site_covs)
-  offset_lambda <- get_offset(lamformula, covs$site_covs)
+  X_lambda <- get_model_matrix(formulas$lambda, covs$site_covs)
+  offset_lambda <- get_offset(formulas$lambda, covs$site_covs)
 
   # Transition probs
   # Drop last period of transition prob design matrices
   # NOTE: this is done outside getDesign for colext
   drop_periods <- rep(1:T, M) == T
   ysc_drop <- covs$yearly_site_covs[!drop_periods,,drop=FALSE]
-  X_gamma <- get_model_matrix(gamformula, ysc_drop)
-  offset_gamma <- get_offset(gamformula, ysc_drop)
+  X_gamma <- get_model_matrix(formulas$gamma, ysc_drop)
+  offset_gamma <- get_offset(formulas$gamma, ysc_drop)
 
-  X_omega <- get_model_matrix(omformula, ysc_drop)
-  offset_omega <- get_offset(omformula, ysc_drop)
+  X_omega <- get_model_matrix(formulas$omega, ysc_drop)
+  offset_omega <- get_offset(formulas$omega, ysc_drop)
 
-  X_iota <- get_model_matrix(iotaformula, ysc_drop)
-  offset_iota <- get_offset(iotaformula, ysc_drop)
+  X_iota <- get_model_matrix(formulas$iota, ysc_drop)
+  offset_iota <- get_offset(formulas$iota, ysc_drop)
 
   # Detection, which differs by model type
   det_covs <- covs$obs_covs
@@ -244,8 +235,8 @@ setMethod("getDesign", "unmarkedFrameDailMadsen",
     covs_dso <- clean_up_covs(umf, drop_final = FALSE)
     det_covs <- covs_dso$yearly_site_covs
   }
-  X_det <-  get_model_matrix(pformula, det_covs)
-  offset_det <- get_offset(pformula, det_covs)
+  X_det <-  get_model_matrix(formulas$det, det_covs)
+  offset_det <- get_offset(formulas$det, det_covs)
 
   # Identify missing values in lambda covs
   has_na_site <- row_has_na(X_lambda)
@@ -343,8 +334,8 @@ setMethod("getDesign", "unmarkedFrameDailMadsen",
     } else
       return("matrix")
   }
-  if(length(all.vars(gamformula)) == 0 & length(all.vars(omformula)) == 0 & 
-     length(all.vars(iotaformula)) == 0){
+  if(length(all.vars(formulas$gamma)) == 0 & length(all.vars(formulas$omega)) == 0 & 
+     length(all.vars(formulas$iota)) == 0){
     go.dims <- "scalar"
   } else {
     go.dims.vec <- apply(Xgo, 2, getGOdims)
@@ -402,7 +393,7 @@ formatDelta <- function(d, yna)
 
 # gdistremoval
 setMethod("getDesign", "unmarkedFrameGDR",
-  function(umf, formula, na.rm=TRUE, return.frames=FALSE){
+  function(umf, formulas, na.rm=TRUE){
 
   M <- numSites(umf)
   T <- umf@numPrimary
@@ -424,28 +415,26 @@ setMethod("getDesign", "unmarkedFrameGDR",
   ysc <- cbind(ysc, sc[rep(1:M, each=T),,drop=FALSE])
   oc <- cbind(oc, ysc[rep(1:nrow(ysc), each=Jrem),,drop=FALSE])
 
-  if(return.frames) return(list(sc=sc, ysc=ysc, oc=oc))
-
-  lam_fixed <- reformulas::nobars(formula$lambdaformula)
+  lam_fixed <- reformulas::nobars(formulas$lambda)
   X_lambda <- model.matrix(lam_fixed,
             model.frame(lam_fixed, sc, na.action=NULL))
 
-  phi_fixed <- reformulas::nobars(formula$phiformula)
+  phi_fixed <- reformulas::nobars(formulas$phi)
   X_phi<- model.matrix(phi_fixed,
             model.frame(phi_fixed, ysc, na.action=NULL))
 
-  dist_fixed <- reformulas::nobars(formula$distanceformula)
+  dist_fixed <- reformulas::nobars(formulas$dist)
   X_dist <- model.matrix(dist_fixed,
             model.frame(dist_fixed, ysc, na.action=NULL))
 
-  rem_fixed <- reformulas::nobars(formula$removalformula)
+  rem_fixed <- reformulas::nobars(formulas$rem)
   X_rem <- model.matrix(rem_fixed,
             model.frame(rem_fixed, oc, na.action=NULL))
 
-  Z_lambda <- get_Z(formula$lambdaformula, sc)
-  Z_phi <- get_Z(formula$phiformula, ysc)
-  Z_dist <- get_Z(formula$distanceformula, ysc)
-  Z_rem <- get_Z(formula$removalformula, oc)
+  Z_lambda <- get_Z(formulas$lambda, sc)
+  Z_phi <- get_Z(formulas$phi, ysc)
+  Z_dist <- get_Z(formulas$dist, ysc)
+  Z_rem <- get_Z(formulas$rem, oc)
  
   # Check if there are missing yearlySiteCovs
   ydist_mat <- apply(matrix(yDist, nrow=M*T, byrow=TRUE), 1, function(x) any(is.na(x)))
@@ -484,13 +473,10 @@ setMethod("getDesign", "unmarkedFrameGDR",
 # there are 3 observation formula which are stored in V (true positive detections),
 # U (false positive detections), and W (b or probability detetion is certain)
 setMethod("getDesign", "unmarkedFrameOccuFP",
-  function(umf, detformula, FPformula, Bformula, stateformula, na.rm = TRUE){
+  function(umf, formulas, na.rm = TRUE){
 
   # Process state and true detection with generic umf method
-  # Combine detection and state formulas
-  comb_form <- list(as.name("~"), detformula, stateformula[[2]])
-  comb_form <- as.formula(as.call(comb_form))
-  out <- methods::callNextMethod(umf, formula = comb_form, na.rm = FALSE)
+  out <- methods::callNextMethod(umf, formulas = formulas, na.rm = FALSE)
  
   M <- numSites(umf)
   J <- obsNum(umf)
@@ -500,12 +486,12 @@ setMethod("getDesign", "unmarkedFrameOccuFP",
   covs <- clean_up_covs(umf)
   
   # Model matrix and offset for false positives
-  X_fp <- get_model_matrix(FPformula, covs$obs_covs)
-  offset_fp <- get_offset(FPformula, covs$obs_covs)
+  X_fp <- get_model_matrix(formulas$fp, covs$obs_covs)
+  offset_fp <- get_offset(formulas$fp, covs$obs_covs)
 
   # Model matrices and offset for b (prob detection is certain)
-  X_b <- get_model_matrix(Bformula, covs$obs_covs)
-  offset_b <- get_offset(Bformula, covs$obs_covs)
+  X_b <- get_model_matrix(formulas$b, covs$obs_covs)
+  offset_b <- get_offset(formulas$b, covs$obs_covs)
 
   # Check missing values in FP and b
   has_na <- row_has_na(cbind(X_fp, X_b))
diff --git a/R/gmultmix.R b/R/gmultmix.R
index 29406d7e..c077edd4 100644
--- a/R/gmultmix.R
+++ b/R/gmultmix.R
@@ -11,13 +11,9 @@ engine <- match.arg(engine, c("C", "R"))
 
 mixture <- match.arg(mixture)
 
-formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula,
-    pformula = pformula)
-check_no_support(formlist)
-form <- as.formula(paste(unlist(formlist), collapse=" "))
-D <- getDesign(data, formula = form)
-X_lambda <- D$X_state
-offset_lambda <- D$offset_state
+formulas <- list(lambda = lambdaformula, phi = phiformula, det = pformula)
+check_no_support(formulas)
+D <- getDesign(data, formulas)
 y <- D$y  # MxJT
 
 K <- check_K_multinomial(K, K_adjust = 100, y, data@numPrimary)
@@ -39,9 +35,9 @@ yt <- apply(y, 1:2, function(x) {
 
 piFun <- data@piFun
 
-lamPars <- colnames(X_lambda)
+lamPars <- colnames(D$X_lambda)
 detPars <- colnames(D$X_det)
-nLP <- ncol(X_lambda)
+nLP <- ncol(D$X_lambda)
 if(T==1) {
     nPP <- 0
     phiPars <- character(0)
@@ -72,7 +68,7 @@ for(i in 1:M) {
     }
 
 nll_R <- function(pars) {
-    lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+    lambda <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
     if(T==1)
         phi <- 1
     else if(T>1)
@@ -134,7 +130,7 @@ if(engine=="R"){
   Kmin <- apply(yt, 1, max, na.rm=TRUE)
 
   nll <- function(params) {
-    nll_gmultmix(params, n_param, y_long, mixture_code, piFun, X_lambda, offset_lambda,
+    nll_gmultmix(params, n_param, y_long, mixture_code, piFun, D$X_lambda, D$offset_lambda,
                  D$X_phi, D$offset_phi, D$X_det, D$offset_det, k, lfac.k, lfac.kmytC,
                  kmytC, Kmin, threads)
   }
@@ -194,7 +190,7 @@ if(identical(mixture,"ZIP")) {
 }
 
 umfit <- new("unmarkedFitGMM", fitType = "gmn",
-    call = match.call(), formula = form, formlist = formlist,
+    call = match.call(), formlist = formulas,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
     AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll,
     mixture=mixture, K=K)
diff --git a/R/goccu.R b/R/goccu.R
index 75c922de..c73bafc0 100644
--- a/R/goccu.R
+++ b/R/goccu.R
@@ -17,21 +17,13 @@ goccu <- function(psiformula, phiformula, pformula, data,
   psiInvLink <- ifelse(linkPsi=="cloglog", "cloglog", "logistic")
   psiLinkGrad <- ifelse(linkPsi=="cloglog", "cloglog.grad", "logistic.grad")
 
-  # Pass phiformula as gamma/eps formula so it will be applied to
-  # yearlySiteCovs in getDesign
-  formlist <- list(psiformula=psiformula, phiformula=phiformula,
-                  pformula=pformula)
-
-  formula <- as.formula(paste(unlist(formlist), collapse=" "))
+  formulas <- list(psi=psiformula, phi=phiformula, det=pformula)
 
   data@y[data@y > 1] <- 1
  
   class(data) <- "unmarkedFrameGOccu"
 
-  # handle offsets
-
-  gd <- getDesign(data, formula = formula)
-  X_psi <- gd$X_state
+  gd <- getDesign(data, formulas)
   y <- gd$y
 
   M <- nrow(y)
@@ -78,7 +70,7 @@ goccu <- function(psiformula, phiformula, pformula, data,
 
   # Bundle data for TMB
   dataList <- list(y=y, T=T, link=ifelse(linkPsi=='cloglog', 1, 0), 
-                   X_psi=X_psi, X_phi=gd$X_phi, X_det=gd$X_det,
+                   X_psi=gd$X_psi, X_phi=gd$X_phi, X_det=gd$X_det,
                    n_possible=n_possible,
                    alpha_potential=alpha_potential, alpha_drop = alpha_drop,
                    known_present=known_present, known_available=known_available, 
@@ -86,7 +78,7 @@ goccu <- function(psiformula, phiformula, pformula, data,
 
   # Provide dimensions and starting values for parameters
   # This part should change to be more like occu() if we add random effects
-  psi_ind <- 1:ncol(X_psi)
+  psi_ind <- 1:ncol(gd$X_psi)
   phi_ind <- 1:ncol(gd$X_phi) + max(psi_ind)
   det_ind <- 1:ncol(gd$X_det) + max(phi_ind)
   nP <- max(det_ind)
@@ -106,14 +98,14 @@ goccu <- function(psiformula, phiformula, pformula, data,
 
   covMat <- invertHessian(opt, nP, se)
   ests <- opt$par
-  names(ests) <- c(colnames(X_psi), colnames(gd$X_phi), colnames(gd$X_det))
+  names(ests) <- c(colnames(gd$X_psi), colnames(gd$X_phi), colnames(gd$X_det))
   fmAIC <- 2 * opt$value + 2 * nP
 
 
   psi_est <- unmarkedEstimate(name = "Occupancy", short.name = "psi",
                               estimates = ests[psi_ind],
                               covMat = covMat[psi_ind, psi_ind, drop=FALSE],
-                              fixed = 1:ncol(X_psi),
+                              fixed = 1:ncol(gd$X_psi),
                               invlink = psiInvLink,
                               invlinkGrad = psiLinkGrad,
                               randomVarInfo=list()
@@ -141,7 +133,7 @@ goccu <- function(psiformula, phiformula, pformula, data,
 
   # Create unmarkedFit object--------------------------------------------------
   umfit <- new("unmarkedFitGOccu", fitType = "goccu", call = match.call(),
-                 formula = formula, formlist=formlist, data = data,
+                 formlist=formulas, data = data,
                  sitesRemoved = gd$removed.sites,
                  estimates = estimate_list, AIC = fmAIC, opt = opt,
                  negLogLike = opt$value,
@@ -153,18 +145,6 @@ goccu <- function(psiformula, phiformula, pformula, data,
 
 # Methods
 
-setMethod("predict_inputs_from_umf", "unmarkedFitGOccu",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm=na.rm)
-  X_idx <- switch(type, psi="X_state", phi="X_phi", det="X_det")
-  list(X=designMats[[X_idx]], offset=NULL)
-})
-
-setMethod("get_formula", "unmarkedFitGOccu", function(object, type, ...){
-  fl <- object@formlist
-  switch(type, psi=fl$psiformula, phi=fl$phiformula, det=fl$pformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitGOccu", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
   datatype <- switch(type, psi='site_covs', phi='yearly_site_covs',
@@ -183,11 +163,10 @@ setMethod("getP_internal", "unmarkedFitGOccu", function(object){
 setMethod("fitted_internal", "unmarkedFitGOccu", function(object){
   # TODO: Use predict here
   M <- numSites(object@data)
-  JT <- obsNum(object@data)  
-  gd <- getDesign(object@data, object@formula, na.rm=FALSE)
-  X_psi <- gd$X_state
+  JT <- obsNum(object@data) 
+  gd <- getDesign(object@data, object@formlist, na.rm=FALSE)
 
-  psi <- drop(plogis(X_psi %*% coef(object, "psi")))
+  psi <- drop(plogis(gd$X_psi %*% coef(object, "psi")))
   psi <- matrix(psi, nrow=M, ncol=JT)
 
   phi <- drop(plogis(gd$X_phi %*% coef(object, "phi")))
@@ -208,11 +187,10 @@ setMethod("ranef_internal", "unmarkedFitGOccu", function(object, ...){
   T <- object@data@numPrimary
   J <- JT / T
 
-  gd <- getDesign(object@data, object@formula, na.rm=FALSE)
+  gd <- getDesign(object@data, object@formlist, na.rm=FALSE)
   y_array <- array(t(gd$y), c(J, T, M))
-  X_psi <- gd$X_state
 
-  psi <- drop(plogis(X_psi %*% coef(object, "psi")))
+  psi <- drop(plogis(gd$X_psi %*% coef(object, "psi")))
   phi <- drop(plogis(gd$X_phi %*% coef(object, "phi")))
   phi <- matrix(phi, nrow=M, ncol=T, byrow=TRUE)
   p <- getP(object)
@@ -302,8 +280,8 @@ setMethod("get_fitting_function", "unmarkedFrameGOccu",
 setMethod("rebuild_call", "unmarkedFitGOccu", function(object){           
   cl <- object@call
   cl[["data"]] <- quote(object@data)
-  cl[["psiformula"]] <- object@formlist$psiformula
-  cl[["phiformula"]] <- object@formlist$phiformula
-  cl[["pformula"]] <- object@formlist$pformula
+  cl[["psiformula"]] <- object@formlist$psi
+  cl[["phiformula"]] <- object@formlist$phi
+  cl[["pformula"]] <- object@formlist$det
   cl
 })
diff --git a/R/gpcount.R b/R/gpcount.R
index d0b97483..54a706b3 100644
--- a/R/gpcount.R
+++ b/R/gpcount.R
@@ -9,13 +9,10 @@ if(!is(data, "unmarkedFrameGPC"))
 mixture <- match.arg(mixture)
 engine <- match.arg(engine)
 
-formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula,
-    pformula = pformula)
-check_no_support(formlist)
-form <- as.formula(paste(unlist(formlist), collapse=" "))
-D <- getDesign(data, formula = form)
-X_lambda <- D$X_state
-offset_lambda <- D$offset_state
+formulas <- list(lambda = lambdaformula, phi = phiformula, det = pformula)
+check_no_support(formulas)
+
+D <- getDesign(data, formulas)
 ym <- D$y  # MxJT
 
 if(missing(K) || is.null(K)) {
@@ -32,9 +29,9 @@ J <- numY(data) / T
 
 y <- array(ym, c(I, J, T))
 
-lamPars <- colnames(X_lambda)
+lamPars <- colnames(D$X_lambda)
 detPars <- colnames(D$X_det)
-nLP <- ncol(X_lambda)
+nLP <- ncol(D$X_lambda)
 nPP <- ncol(D$X_phi)
 phiPars <- colnames(D$X_phi)
 nDP <- ncol(D$X_det)
@@ -45,7 +42,7 @@ if(!missing(starts) && length(starts) != nP)
 if(identical(engine, "R")) {
 # Minus negative log-likelihood
 nll <- function(pars) {
-    lam <- exp(X_lambda %*% pars[1:nLP] + offset_lambda)
+    lam <- exp(D$X_lambda %*% pars[1:nLP] + D$offset_lambda)
     phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi)
     phi <- matrix(phi, I, T, byrow=TRUE)
     p <- plogis(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + D$offset_det)
@@ -96,8 +93,8 @@ if(identical(engine, "C")) {
         log.alpha <- 1
         if(mixture %in% c("NB", "ZIP"))
             log.alpha <- pars[nP]
-        nll_gpcount(ym, X_lambda, D$X_phi, D$X_det, beta_lambda, beta_phi, beta_det,
-                    log.alpha, offset_lambda, D$offset_phi, D$offset_det,
+        nll_gpcount(ym, D$X_lambda, D$X_phi, D$X_det, beta_lambda, beta_phi, beta_det,
+                    log.alpha, D$offset_lambda, D$offset_phi, D$offset_det,
                     as.integer(K), mixture, T, threads)
     }
 }
@@ -147,7 +144,7 @@ if(identical(mixture,"ZIP")) {
 }
 
 umfit <- new("unmarkedFitGPC", fitType = "gpcount",
-    call = match.call(), formula = form, formlist = formlist,
+    call = match.call(), formlist = formulas,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
     AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll,
     mixture=mixture, K=K)
diff --git a/R/multinomPois.R b/R/multinomPois.R
index 7c144df3..fbd789d7 100644
--- a/R/multinomPois.R
+++ b/R/multinomPois.R
@@ -7,8 +7,10 @@ multinomPois <- function(formula, data, starts, method = "BFGS",
     if(!is(data, "unmarkedFrameMPois"))
 		    stop("Data is not a data frame or unmarkedFrame.")
     engine <- match.arg(engine, c("C", "R", "TMB"))
-    if(any(sapply(split_formula(formula), has_random))) engine <- "TMB"
-    dm <- getDesign(data, formula)
+    formulas <- split_formula(formula)
+    names(formulas) <- c("det", "state")
+    if(any(sapply(formulas, has_random))) engine <- "TMB"
+    dm <- getDesign(data, formulas)
     y <- dm$y
 
     J <- ncol(y)
@@ -77,9 +79,8 @@ multinomPois <- function(formula, data, starts, method = "BFGS",
 
     } else if(engine == "TMB"){
 
-      forms <- split_formula(formula)
       obs_all <- add_covariates(obsCovs(data), siteCovs(data), numSites(data)*obsNum(data))
-      inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)),
+      inps <- get_ranef_inputs(formulas, list(det=obs_all, state=siteCovs(data)),
                                list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")])
 
       if(!piFun%in%c('doublePiFun','removalPiFun','depDoublePiFun')){
@@ -104,8 +105,8 @@ multinomPois <- function(formula, data, starts, method = "BFGS",
       det_coef <- get_coef_info(tmb_out$sdr, "det", detParms, pIdx)
 
       # 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]], obs_all)
+      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, obs_all)
 
     }
 
@@ -123,7 +124,7 @@ multinomPois <- function(formula, data, starts, method = "BFGS",
         det=detEstimates))
 
     umfit <- new("unmarkedFitMPois", fitType = "multinomPois",
-        call = match.call(), formula = formula, data = data,
+        call = match.call(), formula = formula, formlist = formulas, data = data,
         estimates = estimateList, sitesRemoved = dm$removed.sites,
         AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, TMB=tmb_mod)
 
diff --git a/R/multmixOpen.R b/R/multmixOpen.R
index fca1448c..c9fa67b1 100644
--- a/R/multmixOpen.R
+++ b/R/multmixOpen.R
@@ -24,12 +24,11 @@ multmixOpen <- 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
 
   M <- nrow(y)
@@ -248,7 +247,7 @@ multmixOpen <- function(lambdaformula, gammaformula, omegaformula, pformula,
   }
 
   umfit <- new("unmarkedFitMMO", fitType = "multmixOpen",
-      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)
diff --git a/R/nmixTTD.R b/R/nmixTTD.R
index 2192e8fc..3df5e704 100644
--- a/R/nmixTTD.R
+++ b/R/nmixTTD.R
@@ -19,12 +19,11 @@ nmixTTD <- function(stateformula=~1, detformula=~1, data, K=100,
   mixture <- match.arg(mixture)
   ttdDist <- match.arg(ttdDist)
 
-  formula <- list(stateformula, ~1, ~1, detformula)
-  check_no_support(formula)
-  formula <- as.formula(paste(unlist(formula),collapse=" "))
+  formulas <- list(state = stateformula, col = ~1, ext = ~1, det = detformula)
+  check_no_support(formulas)
 
   #Process input data----------------------------------------------------------
-  dm <- getDesign(data, formula)
+  dm <- getDesign(data, formulas)
   y <- dm$y
   removed <- dm$removed.sites
 
@@ -171,9 +170,7 @@ nmixTTD <- function(stateformula=~1, detformula=~1, data, K=100,
 
   umfit <- new("unmarkedFitNmixTTD", fitType = "nmixTTD",
                call = match.call(),
-               formula = formula,
-               stateformula = stateformula,
-               detformula = detformula,
+               formlist = formulas,
                K = K,
                data = data, sitesRemoved = removed,
                estimates = estimateList,
diff --git a/R/occu.R b/R/occu.R
index 3f65a0ab..ed3b990a 100644
--- a/R/occu.R
+++ b/R/occu.R
@@ -10,7 +10,9 @@ occu <- function(formula, data, knownOcc = numeric(0),
     stop("Data is not an unmarkedFrameOccu object.")
 
   engine <- match.arg(engine, c("C", "R", "TMB"))
-  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"
   if(length(knownOcc)>0 & engine == "TMB"){
     stop("TMB engine does not support knownOcc argument", call.=FALSE)
   }
@@ -21,7 +23,7 @@ occu <- function(formula, data, knownOcc = numeric(0),
   psiLinkGrad <- ifelse(linkPsi=="cloglog", "cloglog.grad", "logistic.grad")
 
   # Format input data----------------------------------------------------------
-  dm <- getDesign(data, formula)
+  dm <- getDesign(data, formulas)
   y <- truncateToBinary(dm$y)
 
   # Re-format some variables for C and R engines
@@ -95,9 +97,8 @@ occu <- function(formula, data, knownOcc = numeric(0),
   } else if(identical(engine, "TMB")){
 
     # Set up TMB input data
-    forms <- split_formula(formula)
     obs_all <- add_covariates(obsCovs(data), siteCovs(data), length(getY(data)))
-    inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)),
+    inps <- get_ranef_inputs(formulas, list(det=obs_all, state=siteCovs(data)),
                              list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")])
 
     tmb_dat <- c(list(y=y, no_detect=nd, link=ifelse(linkPsi=="cloglog",1,0),
@@ -117,9 +118,9 @@ occu <- function(formula, data, knownOcc = numeric(0),
     det_coef <- get_coef_info(tmb_out$sdr, "det", detParms, (nOP+1):nP)
 
     # Organize random effect estimates
-    state_rand_info <- get_randvar_info(tmb_out$sdr, "state", forms[[2]],
+    state_rand_info <- get_randvar_info(tmb_out$sdr, "state", formulas$state,
                                         siteCovs(data))
-    det_rand_info <- get_randvar_info(tmb_out$sdr, "det", forms[[1]],
+    det_rand_info <- get_randvar_info(tmb_out$sdr, "det", formulas$det,
                                       obs_all)
 
     }
@@ -147,7 +148,7 @@ occu <- function(formula, data, knownOcc = numeric(0),
 
   # Create unmarkedFit object--------------------------------------------------
   umfit <- new("unmarkedFitOccu", fitType = "occu", call = match.call(),
-                 formula = formula, data = data,
+                 formula = formula, formlist = formulas, data = data,
                  sitesRemoved = dm$removed.sites,
                  estimates = estimateList, AIC = fmAIC, opt = fm,
                  negLogLike = fm$value,
diff --git a/R/occuCOP.R b/R/occuCOP.R
index 0b223dbb..d9165955 100644
--- a/R/occuCOP.R
+++ b/R/occuCOP.R
@@ -22,6 +22,19 @@
 
 # Methods ----------------------------------------------------------------------
 
+## getDesign method ----
+# All this method does is translate the submodel names
+setMethod("getDesign", "unmarkedFrameOccuCOP", function(umf, formulas, na.rm = TRUE){
+  names(formulas)[names(formulas) == "psi"] <- "state"
+  names(formulas)[names(formulas) == "lambda"] <- "det"
+  out <- methods::callNextMethod(umf, formulas = formulas, na.rm = na.rm)
+  nms <- gsub("state", "psi", names(out))
+  nms <- gsub("det", "lambda", nms)
+  names(out) <- nms
+  out
+})
+
+
 ## getL method ----
 setMethod("getL", "unmarkedFrameOccuCOP", function(object) {
   return(object@L)
@@ -80,37 +93,6 @@ setMethod("summary", "unmarkedFrameOccuCOP", function(object,...) {
 })
 
 
-## fl_getY ----
-setMethod("fl_getY", "unmarkedFitOccuCOP", function(fit, ...){
-  fl <- fit@formlist
-  comb_form <- list(as.name("~"), fl$lambdaformula, fl$psiformula[[2]])
-  getDesign(getData(fit), comb_form)$y
-})
-
-
-## predict_inputs_from_umf ----
-setMethod("predict_inputs_from_umf", "unmarkedFitOccuCOP",
-          function(object, type, newdata, na.rm, re.form = NULL) {
-            fl <- object@formlist
-            comb_form <- list(as.name("~"), fl$lambdaformula, fl$psiformula[[2]])
-            designMats = getDesign(umf = newdata,
-                                   formula = comb_form,
-                                   na.rm = na.rm)
-            if (type == "psi") list_els <- c("X_state", "Z_state")
-            if (type == "lambda") list_els <- c("X_det", "Z_det")
-            X <- designMats[[list_els[1]]]
-            if (is.null(re.form)) X <- cbind(X, designMats[[list_els[2]]])
-            return(list(X = X, offset = NULL))
-          })
-
-
-## get_formula ----
-setMethod("get_formula", "unmarkedFitOccuCOP", function(object, type, ...) {
-  fl <- object@formlist
-  switch(type, psi = fl$psiformula, lambda = fl$lambdaformula)
-})
-
-
 ## get_orig_data ----
 setMethod("get_orig_data", "unmarkedFitOccuCOP", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
@@ -135,14 +117,12 @@ setMethod("fitted_internal", "unmarkedFitOccuCOP", function(object) {
   data <- object@data
   M = nrow(getY(data))
   J = ncol(getY(data))
-  fl <- object@formlist
-  comb_form <- list(as.name("~"), fl$lambdaformula, fl$psiformula[[2]])
-  des <- getDesign(data, comb_form, na.rm = FALSE)
+  des <- getDesign(data, object@formlist, na.rm = FALSE)
   estim_psi = as.numeric(do.call(object["psi"]@invlink,
-                                 list(as.matrix(des$X_state %*% coef(object, 'psi')))))
+                                 list(as.matrix(des$X_psi %*% coef(object, 'psi')))))
   estim_lambda = do.call(object["lambda"]@invlink, 
                          list(matrix(
-                           as.numeric(des$X_det %*% coef(object, 'lambda')),
+                           as.numeric(des$X_lambda %*% coef(object, 'lambda')),
                            nrow = M, ncol = J, byrow = T)))
   return(estim_psi * estim_lambda)
 })
@@ -262,8 +242,8 @@ setMethod("ranef_internal", "unmarkedFitOccuCOP", function(object, ...) {
 # Used by update() method
 setMethod("rebuild_call", "unmarkedFitOccuCOP", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["psiformula"]] <- object@formlist$psiformula
-  cl[["lambdaformula"]] <- object@formlist$lambdaformula
+  cl[["psiformula"]] <- object@formlist$psi
+  cl[["lambdaformula"]] <- object@formlist$lambda
   cl
 })
 
@@ -433,12 +413,11 @@ occuCOP <- function(data,
   # Format input data ----------------------------------------------------------
   
   # Retrieve formlist
-  formlist <- mget(c("psiformula", "lambdaformula"))
+  formlist <- list(lambda = lambdaformula, psi = psiformula)
   
   # Get the design matrix (calling the getDesign method for unmarkedFrame)
   # For more informations, see: getMethod("getDesign", "unmarkedFrameOccuCOP")
-  comb_form <- list(as.name("~"), lambdaformula, psiformula[[2]])
-  designMats <- getDesign(umf = data, formula = comb_form, na.rm = TRUE)
+  designMats <- getDesign(umf = data, formulas = formlist, na.rm = TRUE)
   
   # y is the count detection data (matrix of size M sites x J observations)
   y <- designMats$y
@@ -461,16 +440,16 @@ occuCOP <- function(data,
   }
   
   # Xpsi is the fixed effects design matrix for occupancy
-  Xpsi <- designMats$X_state
+  Xpsi <- designMats$X_psi
   
   # Xlambda is the fixed effects design matrix for detection rate
-  Xlambda <- designMats$X_det
+  Xlambda <- designMats$X_lambda
   
   # Zpsi is the random effects design matrix for occupancy
-  Zpsi <- designMats$Z_state
+  Zpsi <- designMats$Z_psi
   
   # Zlambda is the random effects design matrix for detection rate
-  Zlambda <- designMats$Z_det
+  Zlambda <- designMats$Z_lambda
   
   # removed_obs is a M x J matrix of the observations removed from the analysis
   removed_obs <- is.na(y)
@@ -668,9 +647,6 @@ occuCOP <- function(data,
     "unmarkedFitOccuCOP",
     fitType = "occuCOP",
     call = match.call(),
-    formula = as.formula(paste(
-      formlist["lambdaformula"], formlist["psiformula"], collapse = ""
-    )),
     formlist = formlist,
     data = data,
     estimates = estimateList,
diff --git a/R/occuFP.R b/R/occuFP.R
index da9f83e4..cdab9bef 100644
--- a/R/occuFP.R
+++ b/R/occuFP.R
@@ -8,14 +8,16 @@
 
 occuFP <- function(detformula = ~ 1,FPformula = ~ 1,Bformula = ~ 1,stateformula = ~ 1, data, starts,
                  method = "BFGS", se = TRUE, engine = "R", ...) {
-    check_no_support(list(detformula,FPformula,Bformula,stateformula))
+
 
     if(!is(data, "unmarkedFrameOccuFP"))   stop("Data is not an unmarkedFrameOccuFP object.")
 
     type <- data@type
     if(sum(type[2:3])==0)   stop("Only type 1 data. No data types with false positives. Use occu instead.")
 
-    dm <- getDesign(data, detformula,FPformula,Bformula,stateformula)
+    formulas <- list(det = detformula, fp = FPformula, b = Bformula, state = stateformula)
+    check_no_support(formulas)
+    dm <- getDesign(data, formulas)
     y <- dm$y
     if(any(type[1:2]>0)) if(any(y[,1:sum(type[1:2])]>1,na.rm = TRUE))   stop("Values of y for type 1 and type 2 data must be 0 or 1.")
     if(type[3]>0) if(any(y[,1:sum(type[3])]>2,na.rm = TRUE))   stop("Values of y for type 3 data must be 0, 1, or 2.")
@@ -109,8 +111,7 @@ occuFP <- function(detformula = ~ 1,FPformula = ~ 1,Bformula = ~ 1,stateformula
     }
 
     umfit <- new("unmarkedFitOccuFP", fitType = "occuFP", call = match.call(),
-                 detformula = detformula,FPformula = FPformula,Bformula = Bformula,
-                 stateformula = stateformula, formula = ~1, type = type, data = data,
+                 formlist = formulas, type = type, data = data,
                  sitesRemoved = dm$removed.sites,
                  estimates = estimateList, AIC = fmAIC, opt = fm,
                  negLogLike = fm$value,
diff --git a/R/occuMS.R b/R/occuMS.R
index 3d1bf0f4..d01e9be8 100644
--- a/R/occuMS.R
+++ b/R/occuMS.R
@@ -249,7 +249,7 @@ occuMS <- function(detformulas, psiformulas, phiformulas=NULL, data,
                 detformulas = detformulas, psiformulas = psiformulas,
                 phiformulas = phiformulas,
                 parameterization = parameterization,
-                formula = ~1, data = data,
+                data = data,
                 sitesRemoved = gd$removed.sites,
                 estimates = estimateList, AIC = fmAIC, opt = fm,
                 negLogLike = fm$value, nllFun = nll_R)
diff --git a/R/occuMulti.R b/R/occuMulti.R
index 9ce0ed67..525f8c18 100644
--- a/R/occuMulti.R
+++ b/R/occuMulti.R
@@ -134,7 +134,7 @@ occuMulti <- function(detformulas, stateformulas,  data, maxOrder,
 
   umfit <- new("unmarkedFitOccuMulti", fitType = "occuMulti", call = cl,
                 detformulas = detformulas, stateformulas = stateformulas,
-                formula = ~1, data = data,
+                data = data,
                 #sitesRemoved = designMats$removed.sites,
                 estimates = estimateList, AIC = fmAIC, opt = fm,
                 negLogLike = fm$value, nllFun = nll_R)
diff --git a/R/occuPEN.R b/R/occuPEN.R
index ce2b5202..ef567e07 100644
--- a/R/occuPEN.R
+++ b/R/occuPEN.R
@@ -2,7 +2,9 @@
 computeMPLElambda = function(formula, data, knownOcc = numeric(0), starts, method = "BFGS", engine = c("C", "R")){
 
   engine <- match.arg(engine, c("C", "R"))
-  designMats <- getDesign(data, formula)
+  formulas <- split_formula(formula)
+  names(formulas) <- c("det", "state")
+  designMats <- getDesign(data, formulas)
   X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y
   removed <- designMats$removed.sites
   y <- truncateToBinary(y)
@@ -48,7 +50,9 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts,
   if (length(lambdaVec)==1) stop("Must provide more than one lambda for cross-validation.")
 
   engine <- match.arg(engine, c("C", "R"))
-  designMats <- getDesign(data, formula)
+  formulas <- split_formula(formula)
+  names(formulas) <- c("det", "state")
+  designMats <- getDesign(data, formulas)
   X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y
   y <- truncateToBinary(y)
   J <- ncol(y)
@@ -102,7 +106,7 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts,
     occuTrain = data[which(foldAssignments!=fold),] # train on NOT this fold
     occuTest = data[which(foldAssignments==fold),] # test on this fold
 
-    designMats <- getDesign(occuTest, formula)
+    designMats <- getDesign(occuTest, formulas)
     X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y
     removed <- designMats$removed.sites
     offset_state <- designMats$offset_state; offset_det <- designMats$offset_det
@@ -146,7 +150,7 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts,
   occuOut = occuPEN(formula, data, starts=starts, lambda=bestLambda, pen.type=pen.type)
 
   umfit <- new("unmarkedFitOccuPEN_CV", fitType = "occu", call = match.call(),
-                 formula = formula, data = data,
+                 formula = formula, formlist = formulas, data = data,
                  sitesRemoved = designMats$removed.sites,
                  estimates = occuOut@estimates, AIC = occuOut@AIC,
 		 opt = occuOut@opt,
@@ -168,7 +172,9 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts,
 		 ...)
 {
 
-    check_no_support(split_formula(formula))
+    formulas <- split_formula(formula)
+    names(formulas) <- c("det", "state")
+    check_no_support(formulas)
 
     if(!is(data, "unmarkedFrameOccu"))
         stop("Data is not an unmarkedFrameOccu object.")
@@ -179,7 +185,7 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts,
 
     engine <- match.arg(engine, c("C", "R"))
 
-    designMats <- getDesign(data, formula)
+    designMats <- getDesign(data, formulas)
     X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y
 
     if (ncol(X_state)==1 & pen.type=="MPLE") stop("MPLE requires occupancy covariates.")
@@ -296,7 +302,7 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts,
     estimateList <- unmarkedEstimateList(list(state=state, det=det))
 
     umfit <- new("unmarkedFitOccuPEN", fitType = "occu", call = match.call(),
-                 formula = formula, data = data,
+                 formula = formula, formlist = formulas, data = data,
                  sitesRemoved = designMats$removed.sites,
                  estimates = estimateList, AIC = fmAIC, opt = opt,
                  negLogLike = fm$value,
diff --git a/R/occuRN.R b/R/occuRN.R
index 71f2649d..492ef805 100644
--- a/R/occuRN.R
+++ b/R/occuRN.R
@@ -4,13 +4,15 @@
 occuRN <- function(formula, data, K = 25, starts, method = "BFGS",
     se = TRUE, engine=c("C","R"), threads = 1, ...)
 {
-  check_no_support(split_formula(formula))
+  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, formula)
+  dm <- getDesign(data, formulas)
   y <- dm$y
   y <- truncateToBinary(y)
 
@@ -95,7 +97,7 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS",
           det=detEstimates))
 
   umfit <- new("unmarkedFitOccuRN", fitType = "occuRN",
-      call = match.call(), formula = formula, data = data,
+      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)
 
diff --git a/R/occuTTD.R b/R/occuTTD.R
index fde1b346..a98ee9af 100644
--- a/R/occuTTD.R
+++ b/R/occuTTD.R
@@ -16,9 +16,8 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
   ttdDist <- match.arg(ttdDist, c("exp","weibull"))
   linkPsi <- match.arg(linkPsi, c("logit","cloglog"))
 
-  formula <- list(psiformula, gammaformula, epsilonformula, detformula)
-  check_no_support(formula)
-  formula <- as.formula(paste(unlist(formula),collapse=" "))
+  formulas <- list(psi = psiformula, col = gammaformula, ext = epsilonformula, det = detformula)
+  check_no_support(formulas)
 
   #Psi link function
   linkFunc <- plogis
@@ -31,8 +30,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
   }
 
   #Process input data----------------------------------------------------------
-  dm <- getDesign(data, formula)
-  X_psi <- dm$X_state
+  dm <- getDesign(data, formulas)
   X_col <- dm$X_col; X_ext <- dm$X_ext
   y <- dm$y
 
@@ -51,7 +49,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
 
   #Organize parameters---------------------------------------------------------
   detParms <- colnames(dm$X_det); nDP <- ncol(dm$X_det)
-  occParms <- colnames(X_psi); nOP <- ncol(X_psi)
+  occParms <- colnames(dm$X_psi); nOP <- ncol(dm$X_psi)
   psi_inds <- 1:nOP
 
   gamParms <- NULL; nGP <- 0; col_inds <- c(0,0)
@@ -78,7 +76,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
   nll_R <- function(params){
 
     #Get occupancy and detection parameters
-    psi <- linkFunc(X_psi %*% params[psi_inds])
+    psi <- linkFunc(dm$X_psi %*% params[psi_inds])
     psi <- cbind(1-psi, psi)
     lam <- exp(dm$X_det %*% params[det_inds])
 
@@ -134,7 +132,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
 
   nll_C <- function(params){
     nll_occuTTD(
-          params, yvec, delta, X_psi, dm$X_det, X_col, X_ext,
+          params, yvec, delta, dm$X_psi, dm$X_det, X_col, X_ext,
           range(psi_inds)-1, range(det_inds)-1,
           range(col_inds)-1, range(ext_inds)-1,
           linkPsi, ttdDist, N, T, J, naflag
@@ -199,11 +197,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1,
 
   umfit <- new("unmarkedFitOccuTTD", fitType = "occuTTD",
                call = match.call(),
-               formula = formula,
-               psiformula = psiformula,
-               gamformula = gammaformula,
-               epsformula = epsilonformula,
-               detformula = detformula,
+               formlist = formulas,
                data = data, sitesRemoved = dm$removed.sites,
                estimates = estimateList,
                AIC = fmAIC, opt = fm, negLogLike = fm$value,
diff --git a/R/pcount.R b/R/pcount.R
index 955903b1..a7839447 100644
--- a/R/pcount.R
+++ b/R/pcount.R
@@ -14,12 +14,14 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
         stop("Data is not an unmarkedFramePCount object.")
 
     engine <- match.arg(engine, c("C", "R", "TMB"))
-    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"
     if(identical(mixture, "ZIP") & engine == "R")
         stop("ZIP mixture not available for R engine")
 
     # Generate design matrices-------------------------------------------------
-    dm <- getDesign(data, formula)
+    dm <- getDesign(data, formulas)
     y <- dm$y
 
     # Set up parameter names and indices---------------------------------------
@@ -115,9 +117,8 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
     } else if(engine == "TMB"){
 
       # Set up TMB input data
-      forms <- split_formula(formula)
       obs_all <- add_covariates(obsCovs(data), siteCovs(data), length(getY(data)))
-      inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)),
+      inps <- get_ranef_inputs(formulas, list(det=obs_all, state=siteCovs(data)),
                                list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")])
 
       tmb_dat <- c(list(y=y, K=K, Kmin=Kmin, mixture=mixture_code,
@@ -144,8 +145,8 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
       }
 
       # 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]], obs_all)
+      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, obs_all)
 
     }
 
@@ -181,7 +182,7 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
 
     # Create unmarkedFit object------------------------------------------------
     umfit <- new("unmarkedFitPCount", fitType="pcount", call=match.call(),
-                 formula = formula, data = data,
+                 formula = formula, formlist = formulas, data = data,
                  sitesRemoved = dm$removed.sites,
                  estimates = estimateList, AIC = fmAIC, opt = fm,
                  negLogLike = fm$value,
diff --git a/R/pcount.spHDS.R b/R/pcount.spHDS.R
index 489e492f..b60e2302 100644
--- a/R/pcount.spHDS.R
+++ b/R/pcount.spHDS.R
@@ -11,7 +11,9 @@ pcount.spHDS<-
 function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
     method = "BFGS", se = TRUE, ...)
 {
-    check_no_support(split_formula(formula))
+    formulas <- split_formula(formula)
+    names(formulas) <- c("det", "state")
+    check_no_support(formulas)
 
     mixture <- match.arg(mixture, c("P", "NB", "ZIP"))
     if (!is(data, "unmarkedFramePCount"))
@@ -20,7 +22,7 @@ function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
    #engine <- match.arg(engine, c("C", "R"))
     if (identical(mixture, "ZIP") & identical(engine, "R"))
         stop("ZIP mixture not available for engine='R'")
-    dm <- getDesign(data, formula)
+    dm <- getDesign(data, formulas)
     y <- dm$y
 
     NAmat <- is.na(y)
@@ -111,7 +113,7 @@ function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
                 nP]), invlink = "logistic", invlinkGrad = "logistic.grad")
     }
     umfit <- new("unmarkedFitPCount", fitType = "pcount", call = match.call(),
-        formula = formula, data = data, sitesRemoved = dm$removed.sites,
+        formula = formula, formlist = formulas, data = data, sitesRemoved = dm$removed.sites,
         estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value,
         nllFun = nll, K = K, mixture = mixture)
     return(umfit)
diff --git a/R/pcountOpen.R b/R/pcountOpen.R
index 2e22860a..6b075050 100644
--- a/R/pcountOpen.R
+++ b/R/pcountOpen.R
@@ -15,11 +15,10 @@ fix <- match.arg(fix)
 ##     stop("lambdaformula and omegaformula must be identical for notrend model")
 if((identical(dynamics, "constant") || identical(dynamics, "notrend")) & immigration)
     stop("You can not include immigration in the constant or notrend models")
-formlist <- list(lambdaformula=lambdaformula, gammaformula=gammaformula,
-    omegaformula=omegaformula, pformula=pformula, iotaformula=iotaformula)
-check_no_support(formlist)
-formula <- as.formula(paste(unlist(formlist), collapse=" "))
-D <- getDesign(data, formula)
+formulas <- list(lambda=lambdaformula, gamma=gammaformula, omega=omegaformula, 
+                 det=pformula, iota=iotaformula)
+check_no_support(formulas)
+D <- getDesign(data, formulas)
 y <- D$y
 
 delta <- D$delta; go.dims <- D$go.dims
@@ -226,7 +225,7 @@ if(identical(mixture, "ZIP")) {
         invlinkGrad = "logistic.grad")
     }
 umfit <- new("unmarkedFitPCO", fitType = "pcountOpen",
-    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, immigration = immigration, fix = fix)
diff --git a/R/predict.R b/R/predict.R
index 0f08e73c..c3f8b485 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -126,32 +126,32 @@ check_type <- function(mod, type){
   stop("Valid types are ", paste(opts, collapse=", "), call.=FALSE)
 }
 
-# Get X and offset when newdata is umf
 setMethod("predict_inputs_from_umf", "unmarkedFit",
   function(object, type, newdata, na.rm, re.form){
-  designMats <- getDesign(newdata, object@formula, na.rm = na.rm)
-  if(type == "state") list_els <- c("X_state","Z_state","offset_state")
-  if(type == "det") list_els <- c("X_det","Z_det","offset_det")
-  if(type == "scale"){ # no covariates
-    n <- nrow(designMats$X_det)
-    return(list(X = matrix(1, nrow=n, ncol=1), offset = rep(0, n)))
+
+  dm <- getDesign(newdata, object@formlist, na.rm = na.rm)
+  select_names <- paste0(c("X_", "Z_", "offset_"), type)
+  select <- dm[select_names]
+  names(select) <- c("X", "Z", "offset")
+
+  # Scalar parameters like scale, alpha
+  if(type %in% names(object) & !type %in% names(object@formlist)){
+    return(list(X = matrix(1, nrow=1, ncol=1), offset = NULL))
   }
 
-  X <- designMats[[list_els[1]]]
-  if(is.null(re.form)) X <- cbind(X, designMats[[list_els[2]]])
-  offset <- designMats[[list_els[3]]]
+  # Add random effects matrix if needed
+  if(is.null(re.form) & has_random(object@formlist[[type]])){
+    stopifnot(!is.null(select$Z))
+    select$X <- cbind(select$X, select$Z)
+  }
 
-  list(X=X, offset=offset)
+  list(X = select$X, offset = select$offset)
 })
 
+
 # Get correct individual formula based on type
 setMethod("get_formula", "unmarkedFit", function(object, type, ...){
-  if(type == "state"){
-    return(as.formula(paste("~", object@formula[3], sep="")))
-  } else if(type == "det"){
-    return(as.formula(object@formula[[2]]))
-  }
-  NULL
+  object@formlist[[type]]
 })
 
 # When newdata is data.frame/raster, get original dataset
@@ -320,19 +320,6 @@ setMethod("predict_by_chunk", "unmarkedFitPCount",
 
 # colext methods---------------------------------------------------------------
 
-setMethod("predict_inputs_from_umf", "unmarkedFitColExt",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm = na.rm)
-  list_el <- switch(type, psi="X_state", col="X_col", ext="X_ext", det="X_det")
-  # colext doesn't support offsets
-  list(X=designMats[[list_el]], offset=NULL)
-})
-
-setMethod("get_formula", "unmarkedFitColExt", function(object, type, ...){
-  switch(type, psi=object@psiformula, col=object@gamformula,
-                     ext=object@epsformula, det=object@detformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitColExt", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=TRUE)
   datatype <- switch(type, psi='site_covs', col='yearly_site_covs',
@@ -343,20 +330,6 @@ setMethod("get_orig_data", "unmarkedFitColExt", function(object, type, ...){
 
 # occuFP methods---------------------------------------------------------------
 
-setMethod("predict_inputs_from_umf", "unmarkedFitOccuFP",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@detformula, object@FPformula,
-                          object@Bformula, object@stateformula, na.rm=na.rm)
-  X_idx <- switch(type, state="X_state", det="X_det", fp="X_fp", b="X_b")
-  off_idx <- paste0("offset_", type)
-  list(X=designMats[[X_idx]], offset=designMats[[off_idx]])
-})
-
-setMethod("get_formula", "unmarkedFitOccuFP", function(object, type, ...){
-  switch(type, state=object@stateformula, det=object@detformula,
-         b=object@Bformula, fp=object@FPformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitOccuFP", function(object, type, ...){
   # Get obs data if fp, b, or det
   new_type <- ifelse(type %in% c("fp", "b"), "det", type)
@@ -384,21 +357,6 @@ setMethod("check_predict_arguments", "unmarkedFitDailMadsen",
   methods::callNextMethod(object, type, newdata)
 })
 
-setMethod("predict_inputs_from_umf", "unmarkedFitDailMadsen",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm=na.rm)
-  X_idx <- switch(type, lambda="X_lambda", gamma="X_gamma", omega="X_omega",
-                  iota="X_iota", det="X_det")
-  off_idx <- paste0("offset_", type)
-  list(X=designMats[[X_idx]], offset=designMats[[off_idx]])
-})
-
-setMethod("get_formula", "unmarkedFitDailMadsen", function(object, type, ...){
-  fl <- object@formlist
-  switch(type, lambda=fl$lambdaformula, gamma=fl$gammaformula,
-         omega=fl$omegaformula, iota=fl$iotaformula, det=fl$pformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitDailMadsen", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=TRUE)
   datatype <- switch(type, lambda='site_covs', gamma='yearly_site_covs',
@@ -441,19 +399,6 @@ setMethod("predict_by_chunk", "unmarkedFitDailMadsen",
 # All inherit from GMM so only one set of methods is required
 # (except GDR which has its own predict method right now)
 
-setMethod("predict_inputs_from_umf", "unmarkedFitGMM",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm=na.rm)
-  X_idx <- switch(type, lambda="X_state", phi="X_phi", det="X_det")
-  off_idx <- switch(type, lambda="offset_state", phi="offset_phi", det="offset_det")
-  list(X=designMats[[X_idx]], offset=designMats[[off_idx]])
-})
-
-setMethod("get_formula", "unmarkedFitGMM", function(object, type, ...){
-  fl <- object@formlist
-  switch(type, lambda=fl$lambdaformula, phi=fl$phiformula, det=fl$pformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitGMM", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
   datatype <- switch(type, lambda='site_covs', phi='yearly_site_covs',
@@ -466,18 +411,6 @@ setMethod("get_orig_data", "unmarkedFitGMM", function(object, type, ...){
 
 # Identical to colext
 
-setMethod("predict_inputs_from_umf", "unmarkedFitOccuTTD",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm = na.rm)
-  list_el <- switch(type, psi="X_state", col="X_col", ext="X_ext", det="X_det")
-  list(X=designMats[[list_el]], offset=NULL)
-})
-
-setMethod("get_formula", "unmarkedFitOccuTTD", function(object, type, ...){
-  switch(type, psi=object@psiformula, col=object@gamformula,
-                     ext=object@epsformula, det=object@detformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitOccuTTD", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=TRUE)
   datatype <- switch(type, psi='site_covs', col='yearly_site_covs',
@@ -488,17 +421,6 @@ setMethod("get_orig_data", "unmarkedFitOccuTTD", function(object, type, ...){
 
 # nmixTTD----------------------------------------------------------------------
 
-setMethod("predict_inputs_from_umf", "unmarkedFitNmixTTD",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formula, na.rm = na.rm)
-  list_el <- switch(type, state="X_state", det="X_det")
-  list(X=designMats[[list_el]], offset=NULL)
-})
-
-setMethod("get_formula", "unmarkedFitNmixTTD", function(object, type, ...){
-  switch(type, state=object@stateformula, det=object@detformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitNmixTTD", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
   datatype <- switch(type, state='site_covs', det='obs_covs')
@@ -508,24 +430,6 @@ setMethod("get_orig_data", "unmarkedFitNmixTTD", function(object, type, ...){
 
 # gdistremoval-----------------------------------------------------------------
 
-setMethod("predict_inputs_from_umf", "unmarkedFitGDR",
-  function(object, type, newdata, na.rm, re.form=NA){
-  designMats <- getDesign(newdata, object@formlist)
-  if(type == "lambda") list_els <- c("X_lambda","Z_lambda")
-  if(type == "phi") list_els <- c("X_phi","Z_phi")
-  if(type == "dist") list_els <- c("X_dist","Z_dist")
-  if(type == "rem") list_els <- c("X_rem", "Z_rem")
-  X <- designMats[[list_els[1]]]
-  if(is.null(re.form)) X <- cbind(X, designMats[[list_els[2]]])
-  list(X=X, offset=NULL)
-})
-
-setMethod("get_formula", "unmarkedFitGDR", function(object, type, ...){
-  fl <- object@formlist
-  switch(type, lambda=fl$lambdaformula, phi=fl$phiformula,
-         dist=fl$distanceformula, rem=fl$removalformula)
-})
-
 setMethod("get_orig_data", "unmarkedFitGDR", function(object, type, ...){
   clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
   datatype <- switch(type, lambda='site_covs', phi='yearly_site_covs',
diff --git a/R/ranef.R b/R/ranef.R
index c3987d25..d6eebf1c 100644
--- a/R/ranef.R
+++ b/R/ranef.R
@@ -66,9 +66,7 @@ setMethod("ranef_internal", "unmarkedFitColExt", function(object){
 # DSO and MMO
 setMethod("ranef_internal", "unmarkedFitDailMadsen", function(object, ...){
     dyn <- object@dynamics
-    formlist <- object@formlist
-    formula <- as.formula(paste(unlist(formlist), collapse=" "))
-    D <- getDesign(object@data, formula, na.rm=FALSE)
+    D <- getDesign(object@data, object@formlist, na.rm=FALSE)
     delta <- D$delta
     deltamax <- max(delta, na.rm=TRUE)
     if(!.hasSlot(object, "immigration")){ #For backwards compatibility
@@ -804,10 +802,7 @@ setMethod("ranef_internal", "unmarkedFitPCount", function(object, ...){
 
 setMethod("ranef_internal", "unmarkedFitPCO", function(object, ...){
     dyn <- object@dynamics
-    
-    formlist <- object@formlist
-    formula <- as.formula(paste(unlist(formlist), collapse=" "))
-    D <- getDesign(object@data, formula, na.rm=FALSE)
+    D <- getDesign(object@data, object@formlist, na.rm=FALSE)
     delta <- D$delta
     deltamax <- max(delta, na.rm=TRUE)
     if(!.hasSlot(object, "immigration")) #For backwards compatibility
diff --git a/R/simulate.R b/R/simulate.R
index 6ce15d5d..95ad10d5 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -383,7 +383,6 @@ setMethod("simulate_internal", "unmarkedFitGDS",
 
 setMethod("simulate_internal", "unmarkedFitGMM",
   function(object, nsim){
-    formula <- object@formula
     umf <- object@data
     mixture <- object@mixture
     y <- umf@y
@@ -444,7 +443,6 @@ setMethod("simulate_internal", "unmarkedFitGMM",
 
 setMethod("simulate_internal", "unmarkedFitGPC",
     function(object, nsim){
-    formula <- object@formula
     umf <- object@data
     mixture <- object@mixture
     y <- umf@y
@@ -883,7 +881,7 @@ simOpenN <- function(object){
     #To partially handle old saved model objects
     fix <- tryCatch(object@fix, error=function(e) "none")
     immigration <- tryCatch(object@immigration, error=function(e) FALSE)
-    delta <- getDesign(umf, object@formula, na.rm = FALSE)$delta  
+    delta <- getDesign(umf, object@formlist, na.rm = FALSE)$delta  
     y <- umf@y
     M <- nrow(y)
     T <- umf@numPrimary
diff --git a/R/unmarkedFit.R b/R/unmarkedFit.R
index 555133a3..cee738f8 100644
--- a/R/unmarkedFit.R
+++ b/R/unmarkedFit.R
@@ -464,13 +464,7 @@ setMethod("hist", "unmarkedFitDS", function(x, lwd=1, lty=1, ...) {
 # Extract detection probs
 setMethod("getFP", "unmarkedFitOccuFP", function(object, na.rm = TRUE)
 {
-  formula <- object@formula
-  detformula <- object@detformula
-  stateformula <- object@stateformula
-  FPformula <- object@FPformula
-  Bformula <- object@Bformula
-  umf <- object@data
-  designMats <- getDesign(umf, detformula,FPformula,Bformula,stateformula, na.rm = na.rm)
+  designMats <- getDesign(getData(object), object@formlist, na.rm = na.rm)
   type = object@type
   y <- designMats$y
   M <- nrow(y)
@@ -486,13 +480,7 @@ setMethod("getFP", "unmarkedFitOccuFP", function(object, na.rm = TRUE)
 
 setMethod("getB", "unmarkedFitOccuFP", function(object, na.rm = TRUE)
 {
-  formula <- object@formula
-  detformula <- object@detformula
-  stateformula <- object@stateformula
-  FPformula <- object@FPformula
-  Bformula <- object@Bformula
-  umf <- object@data
-  designMats <- getDesign(umf, detformula,FPformula,Bformula,stateformula, na.rm = na.rm)
+  designMats <- getDesign(getData(object), object@formlist, na.rm = na.rm)
   y <- designMats$y
   M <- nrow(y)
   J <- ncol(y)
diff --git a/R/unmarkedFitList.R b/R/unmarkedFitList.R
index 75d9d3d6..3906bd19 100644
--- a/R/unmarkedFitList.R
+++ b/R/unmarkedFitList.R
@@ -1,7 +1,7 @@
 #Get y (corrected for missing values) based on data and formulas
 #Approach for some fit functions is different, so this is now a method
 setMethod("fl_getY", "unmarkedFit", function(fit, ...){
-  getDesign(getData(fit), fit@formula)$y
+  getDesign(getData(fit), fit@formlist)$y
 })
 
 setMethod("fl_getY", "unmarkedFitOccuMulti", function(fit, ...){
@@ -15,11 +15,6 @@ setMethod("fl_getY", "unmarkedFitOccuMS", function(fit, ...){
             fit@detformulas, fit@parameterization)$y
 })
 
-setMethod("fl_getY", "unmarkedFitOccuFP", function(fit, ...){
-  getDesign(getData(fit), fit@detformula, fit@FPformula,
-            fit@Bformula, fit@stateformula)$y
-})
-
 # constructor of unmarkedFitList objects
 fitList <- function(..., fits, autoNames=c("object","formula")) {
   autoNames <- match.arg(autoNames)
@@ -220,9 +215,7 @@ setMethod("modSel", "unmarkedFitList",
     colnames(out) <- cNames
     out$model <- names(fits)
     out$formula <- sapply(fits, function(x) {
-          f <- as.character(x@formula)
-          f <- paste(f[2], "~", f[3])
-          f
+          paste(unlist(x@formlist), collapse = " ")
         })
     for(i in 1:length(eNames)) {
         out[,eNames[i]] <- sapply(estList, function(x) x[eNames[i]])
diff --git a/R/update.R b/R/update.R
index e5c77378..7b31d7d9 100644
--- a/R/update.R
+++ b/R/update.R
@@ -42,21 +42,21 @@ setMethod("rebuild_call", "unmarkedFit", function(object){
 
 setMethod("rebuild_call", "unmarkedFitColExt", function(object){ 
   cl <- methods::callNextMethod(object) # calls base class method above
-  cl[["psiformula"]] <- object@psiformula
-  cl[["gammaformula"]] <- object@gamformula
-  cl[["epsilonformula"]] <- object@epsformula
-  cl[["pformula"]] <- object@detformula
+  cl[["psiformula"]] <- object@formlist$psi
+  cl[["gammaformula"]] <- object@formlist$col
+  cl[["epsilonformula"]] <- object@formlist$ext
+  cl[["pformula"]] <- object@formlist$det
   cl
 })
 
 # Covers MMO and PCO
 setMethod("rebuild_call", "unmarkedFitDailMadsen", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["lambdaformula"]] <- object@formlist$lambdaformula
-  cl[["gammaformula"]] <- object@formlist$gammaformula
-  cl[["omegaformula"]] <- object@formlist$omegaformula
-  cl[["pformula"]] <- object@formlist$pformula
-  cl[["iotaformula"]] <- object@formlist$iotaformula
+  cl[["lambdaformula"]] <- object@formlist$lambda
+  cl[["gammaformula"]] <- object@formlist$gamma
+  cl[["omegaformula"]] <- object@formlist$omega
+  cl[["pformula"]] <- object@formlist$det
+  cl[["iotaformula"]] <- object@formlist$iota
   cl[["mixture"]] <- object@mixture
   cl[["dynamics"]] <- object@dynamics
   cl[["immigration"]] <- object@immigration
@@ -87,9 +87,9 @@ setMethod("rebuild_call", "unmarkedFitDSO", function(object){
 # Also covers unmarkedFitGPC
 setMethod("rebuild_call", "unmarkedFitGMM", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["lambdaformula"]] <- object@formlist$lambdaformula
-  cl[["phiformula"]] <- object@formlist$phiformula
-  cl[["pformula"]] <- object@formlist$pformula
+  cl[["lambdaformula"]] <- object@formlist$lambda
+  cl[["phiformula"]] <- object@formlist$phi
+  cl[["pformula"]] <- object@formlist$det
   cl[["mixture"]] <- object@mixture
   if(methods::.hasSlot(object, "K")){
     cl[["K"]] <- object@K
@@ -113,8 +113,8 @@ setMethod("rebuild_call", "unmarkedFitMPois", function(object){
 
 setMethod("rebuild_call", "unmarkedFitNmixTTD", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["stateformula"]] <- object@stateformula
-  cl[["detformula"]] <- object@detformula
+  cl[["stateformula"]] <- object@formlist$state
+  cl[["detformula"]] <- object@formlist$det
   cl[["K"]] <- object@K
   cl
 })
@@ -130,10 +130,10 @@ setMethod("rebuild_call", "unmarkedFitOccu", function(object){
 
 setMethod("rebuild_call", "unmarkedFitOccuFP", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["detformula"]] <- object@detformula
-  cl[["FPformula"]] <- object@FPformula
-  cl[["Bformula"]] <- object@Bformula
-  cl[["stateformula"]] <- object@stateformula
+  cl[["detformula"]] <- object@formlist$det
+  cl[["FPformula"]] <- object@formlist$fp
+  cl[["Bformula"]] <- object@formlist$b
+  cl[["stateformula"]] <- object@formlist$state
   cl
 })
 
@@ -186,9 +186,9 @@ setMethod("rebuild_call", "unmarkedFitOccuRN", function(object){
 
 setMethod("rebuild_call", "unmarkedFitOccuTTD", function(object){ 
   cl <- methods::callNextMethod(object)
-  cl[["psiformula"]] <- object@psiformula
-  cl[["gammaformula"]] <- object@gamformula
-  cl[["epsilonformula"]] <- object@epsformula
-  cl[["detformula"]] <- object@detformula
+  cl[["psiformula"]] <- object@formlist$psi
+  cl[["gammaformula"]] <- object@formlist$col
+  cl[["epsilonformula"]] <- object@formlist$ext
+  cl[["detformula"]] <- object@formlist$det
   cl
 })