Merge remote-tracking branch 'upstream/master' into gdistsamp_ZIP

27 files changed, 603 insertions, 93 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 5454b14..894773a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.2.5.9010
-Date: 2022-12-29
+Version: 1.2.5.9015
+Date: 2023-06-23
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
@@ -32,7 +32,7 @@ Imports:
     stats,
     TMB (>= 1.7.18),
     utils
-Suggests: knitr, rmarkdown, pkgdown, raster, shiny, testthat
+Suggests: knitr, rmarkdown, pkgdown, raster, shiny, terra, testthat
 Description: Fits hierarchical models of animal abundance and occurrence to data collected using survey methods such as point counts, site occupancy sampling, distance sampling, removal sampling, and double observer sampling. Parameters governing the state and observation processes can be modeled as functions of covariates. Reference: Fiske and Chandler (2011) <doi:10.18637/jss.v043.i10>.
 License: GPL (>=3)
 LazyLoad: yes
@@ -46,6 +46,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R'
     'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R'
     'nmixTTD.R'
     'gdistremoval.R'
+    'plotEffects.R'
     'mixedModelTools.R'
     'power.R' 
     'simulate.R'
diff --git a/NAMESPACE b/NAMESPACE
index 448be4f..7273ce1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,7 +7,7 @@ importFrom(stats, confint, fitted, coef, vcov, predict, update, profile,
            pnorm, qchisq, qnorm, quantile, rbinom,
            reshape, rmultinom, rnbinom, rpois, runif, sd, uniroot,
            update.formula, sigma)
-importFrom(graphics, plot, hist, abline)
+importFrom(graphics, plot, hist, abline, axis, lines, points, polygon, segments)
 importFrom(utils, head, read.csv)
 importFrom(grDevices, devAskNewPage, dev.interactive, palette.colors)
 importFrom(MASS, mvrnorm)
@@ -51,7 +51,7 @@ exportMethods(backTransform, coef, confint, coordinates, fitted, getData,
               "siteCovs<-", summary, update, vcov, yearlySiteCovs,
               "yearlySiteCovs<-", "[", smoothed, projected, nonparboot, logLik,
               LRT, ranef, bup, crossVal, posteriorSamples, sigma, randomTerms, 
-              optimizePenalty, unmarkedPowerList)
+              optimizePenalty, unmarkedPowerList, plotEffectsData, plotEffects)
 
 S3method("print", "unmarkedPostSamples")
 
diff --git a/R/gdistremoval.R b/R/gdistremoval.R
index dc91934..b253625 100644
--- a/R/gdistremoval.R
+++ b/R/gdistremoval.R
@@ -219,9 +219,33 @@ setMethod("getDesign", "unmarkedFrameGDR",
   Zphi <- get_Z(formula$phiformula, ysc)
   Zdist <- get_Z(formula$distanceformula, ysc)
   Zrem <- get_Z(formula$removalformula, oc)
+ 
+  # Check if there are missing yearlySiteCovs
+  ydist_mat <- apply(matrix(yDist, nrow=M*T, byrow=TRUE), 1, function(x) any(is.na(x)))
+  yrem_mat <- apply(matrix(yRem, nrow=M*T, byrow=TRUE), 1, function(x) any(is.na(x)))
+  ok_missing_phi_covs <- ydist_mat | yrem_mat
+  missing_phi_covs <- apply(Xphi, 1, function(x) any(is.na(x)))  
+  if(!all(which(missing_phi_covs) %in% which(ok_missing_phi_covs))){
+    stop("Missing yearlySiteCovs values for some observations that are not missing", call.=FALSE)
+  }
+
+  # Check if there are missing dist covs
+  missing_dist_covs <- apply(Xdist, 1, function(x) any(is.na(x)))
+  ok_missing_dist_covs <- ydist_mat
+  if(!all(which(missing_dist_covs) %in% which(ok_missing_dist_covs))){
+    stop("Missing yearlySiteCovs values for some distance observations that are not missing", call.=FALSE)
+  }
 
-  if(any(is.na(Xlam)) | any(is.na(Xphi)) | any(is.na(Xdist)) | any(is.na(Xrem))){
-    stop("gdistremoval does not currently handle missing values in covariates", call.=FALSE)
+  # Check if there are missing rem covs
+  missing_obs_covs <- apply(Xrem, 1, function(x) any(is.na(x)))
+  missing_obs_covs <- apply(matrix(missing_obs_covs, nrow=M*T, byrow=TRUE), 1, function(x) any(x))
+  ok_missing_obs_covs <- yrem_mat
+  if(!all(which(missing_obs_covs) %in% which(ok_missing_obs_covs))){
+    stop("Missing obsCovs values for some removal observations that are not missing", call.=FALSE)
+  }
+    
+  if(any(is.na(Xlam))){
+    stop("gdistremoval does not currently handle missing values in siteCovs", call.=FALSE)
   }
 
   list(yDist=yDist, yRem=yRem, Xlam=Xlam, Xphi=Xphi, Xdist=Xdist, Xrem=Xrem,
@@ -541,14 +565,14 @@ setMethod("fitted", "unmarkedFitGDR", function(object){
     dist[,,t] <- dist[,,t] * p_rem[,rep(t,ncol(dist[,,t]))]
     if(T > 1){
       rem[,,t] <- rem[,,t] * phi[,rep(t, ncol(rem[,,t]))]
-      det[,,t] <- dist[,,t] * phi[,rep(t, ncol(dist[,,t]))]
+      dist[,,t] <- dist[,,t] * phi[,rep(t, ncol(dist[,,t]))]
     }
   }
 
   if(T > 1){
-    rem_final <- rem[,,t]
-    dist_final <- det[,,t]
-    for (t in 1:T){
+    rem_final <- rem[,,1]
+    dist_final <- dist[,,1]
+    for (t in 2:T){
       rem_final <- cbind(rem_final, rem[,,t])
       dist_final <- cbind(dist_final, dist[,,t])
     }
@@ -615,7 +639,7 @@ setMethod("ranef", "unmarkedFitGDR", function(object){
   pr <- apply(cp, c(1,3), sum)
   prRem <- apply(pif, c(1,3), sum)
 
-  post <- array(0, c(M, K+1, T))
+  post <- array(0, c(M, K+1, 1))
   colnames(post) <- 0:K
   for (i in 1:M){
     if(mixture=="P"){
@@ -625,19 +649,26 @@ setMethod("ranef", "unmarkedFitGDR", function(object){
     } else if(mixture=="ZIP"){
       f <- dzip(0:K, lam[i], alpha)
     }
-    g <- rep(1, K+1)
-    for (t in 1:T){
-      if(has_na[i,t]){
-        g <- rep(NA,K+1)
-        next
-      }
-      for (k in 1:(K+1)){
-        g[k] <- g[k] * dbinom(ysum[i,t], k-1, prob=pr[i,t]*prRem[i,t]*phi[i,t],
-                              log=FALSE)
+
+    # All sampling periods at site i have at least one missing value
+    if(all(has_na[i,])){
+      g <- rep(NA,K+1)
+      next
+    } else {
+      # At least one sampling period wasn't missing
+      g <- rep(1, K+1)
+      for (t in 1:T){
+        if(has_na[i,t]){
+          next
+        }
+        for (k in 1:(K+1)){
+          g[k] <- g[k] * dbinom(ysum[i,t], k-1, prob=pr[i,t]*prRem[i,t]*phi[i,t],
+                                log=FALSE)
+        }
       }
     }
     fg <- f*g
-    post[i,,t] <- fg/sum(fg)
+    post[i,,1] <- fg/sum(fg)
   }
 
   new("unmarkedRanef", post=post)
@@ -702,17 +733,22 @@ setMethod("simulate", "unmarkedFitGDR", function(object, nsim, seed=NULL, na.rm=
     yrem <- matrix(NA, M, T*Jrem)
 
     for (m in 1:M){
-      ysum <- rbinom(T, N[m], p_total[m,])
+      ysum <- suppressWarnings(rbinom(T, N[m], p_total[m,]))
 
       ydist_m <- yrem_m <- c()
 
       for (t in 1:T){
-        rem_class <- sample(1:Jrem, ysum[t], replace=TRUE, prob=rem_scaled[m,,t])
-        rem_class <- factor(rem_class, levels=1:Jrem)
-        yrem_m <- c(yrem_m, as.numeric(table(rem_class)))
-        dist_class <- sample(1:Jdist, ysum[t], replace=TRUE, prob=dist_scaled[m,,t])
-        dist_class <- factor(dist_class, levels=1:Jdist)
-        ydist_m <- c(ydist_m, as.numeric(table(dist_class)))
+        if(is.na(ysum[t])){
+          yrem_m <- c(yrem_m, rep(NA, Jrem))
+          ydist_m <- c(ydist_m, rep(NA, Jdist))
+        } else {
+          rem_class <- sample(1:Jrem, ysum[t], replace=TRUE, prob=rem_scaled[m,,t])
+          rem_class <- factor(rem_class, levels=1:Jrem)
+          yrem_m <- c(yrem_m, as.numeric(table(rem_class)))
+          dist_class <- sample(1:Jdist, ysum[t], replace=TRUE, prob=dist_scaled[m,,t])
+          dist_class <- factor(dist_class, levels=1:Jdist)
+          ydist_m <- c(ydist_m, as.numeric(table(dist_class)))
+        }
       }
       stopifnot(length(ydist_m)==ncol(ydist))
       stopifnot(length(yrem_m)==ncol(yrem))
diff --git a/R/getDesign.R b/R/getDesign.R
index 79f20e6..8597b2d 100644
--- a/R/getDesign.R
+++ b/R/getDesign.R
@@ -23,6 +23,11 @@ setMethod("getDesign", "unmarkedFrame",
     } else {
         siteCovs <- siteCovs(umf)
     }
+    miss_vars <-all.vars(stateformula)[!all.vars(stateformula) %in% names(siteCovs)]
+    if(length(miss_vars) > 0){
+      stop(paste("Variable(s)", paste(miss_vars, collapse=", "), 
+                 "not found in siteCovs"), call.=FALSE)
+    }
     X.mf <- model.frame(stateformula, siteCovs, na.action = NULL)
     X <- model.matrix(stateformula, X.mf)
     X.offset <- as.vector(model.offset(X.mf))
@@ -49,6 +54,11 @@ setMethod("getDesign", "unmarkedFrame",
         obsCovs <- cbind(obsCovs, obsNum = as.factor(rep(1:R, M)))
     }
 
+    miss_vars <-all.vars(detformula)[!all.vars(detformula) %in% names(obsCovs)]
+    if(length(miss_vars) > 0){
+      stop(paste("Variable(s)", paste(miss_vars, collapse=", "), 
+                 "not found in obsCovs"), call.=FALSE)
+    }
     V.mf <- model.frame(detformula, obsCovs, na.action = NULL)
     V <- model.matrix(detformula, V.mf)
     V.offset <- as.vector(model.offset(V.mf))
diff --git a/R/plotEffects.R b/R/plotEffects.R
new file mode 100644
index 0000000..c77f5aa
--- /dev/null
+++ b/R/plotEffects.R
@@ -0,0 +1,109 @@
+#type_to_covs <- function(umf, type){
+#  if(type %in% c("state","psi","lam","lambda","sigma","dist")){
+#    return(methods::slot(umf, "siteCovs"))
+#  } else if(type %in% c("det","rem","fp","b")){
+#    return(methods::slot(umf, "obsCovs"))
+#  } else if(type %in% c("phi","transition","col","ext","gamma","omega","iota")){
+#    return(methods::slot(umf, "yearlySiteCovs"))
+#  }
+#  return(NULL)
+#}
+
+get_base_newdata <- function(object, type){
+  covs <- get_orig_data(object, type)
+  out <- lapply(covs, function(x){
+    if(is.numeric(x)){
+      return(median(x, na.rm=TRUE))
+    } else if(is.factor(x)){
+      return(factor(levels(x)[1], levels=levels(x)))
+    }  else {
+      stop("Unknown column type")
+    }
+  })
+  as.data.frame(out)
+}
+
+get_cov_seq <- function(covariate, object, type){
+  cov_values <- get_orig_data(object, type)[[covariate]]
+  if(is.numeric(cov_values)){
+    rng <- range(cov_values, na.rm=TRUE)
+    return(seq(rng[1], rng[2], length.out=100))
+  } else if(is.factor(cov_values)){
+    return(factor(levels(cov_values), levels=levels(cov_values)))
+  } else {
+    stop("Unknown covariate type")
+  }
+}
+
+setGeneric("plotEffectsData", function(object, ...) standardGeneric("plotEffectsData"))
+
+setMethod("plotEffectsData", "unmarkedFit",
+  function(object, type, covariate, level=0.95, ...){
+
+  #umf <- umf_to_factor(object@data)
+  nd <- get_base_newdata(object, type)
+  if(! covariate %in% names(nd)){
+    stop("Covariate not in this submodel", call.=FALSE)
+  }
+  values <- get_cov_seq(covariate, object, type)
+  nd <- nd[rep(1, length(values)),,drop=FALSE]
+  nd[[covariate]] <- values
+
+  pr <- predict(object, type=type, newdata=nd, level=level, ...)
+  pr$covariate <- covariate
+  pr$covariateValue <- values
+  pr
+})
+
+setGeneric("plotEffects", function(object, ...) standardGeneric("plotEffects"))
+
+setMethod("plotEffects", "unmarkedFit",
+  function(object, type, covariate, level=0.95, ...){
+
+  # Get data for plot
+  plot_data <- plotEffectsData(object, type, covariate, level, ...)
+
+  # Is covariate a factor?
+  is_factor <- is.factor(plot_data$covariateValue)
+
+  # If not a factor (i.e., numeric)
+  if(!is_factor){
+    # Setup basic plot structure
+    plot(plot_data$covariateValue, plot_data$Predicted, type='l',
+         xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)),
+         ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T)))
+    #Draw error ribbon
+    polygon(x=c(plot_data$covariateValue, rev(plot_data$covariateValue)),
+            y=c(plot_data$upper, rev(plot_data$lower)),
+            border=NA, col='gray90')
+    #Draw line on top of ribbon
+    lines(plot_data$covariateValue, plot_data$Predicted)
+
+  # If covariate is a factor
+  } else {
+    # Get number of unique factor levels (groups)
+    ngroup <- nrow(plot_data)
+    # Calculate width of error bar tips based on # of factor levels
+    bw <- ngroup/10
+
+    # Add basic plot of points
+    plot(1:ngroup, plot_data$Predicted,
+         xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)),
+         xlim=c(1-bw/2,ngroup+bw/2),
+         xaxt='n',
+         ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T)))
+
+    # Draw new x-axis with factor labels
+    axis(1, at=1:ngroup, labels=plot_data$covariateValue)
+
+    # Draw error bars
+    for (i in 1:ngroup){
+      segments(i, plot_data$lower[i], i, plot_data$upper[i])
+      segments(i-bw/2, plot_data$lower[i], i+bw/2, plot_data$lower[i])
+      segments(i-bw/2, plot_data$upper[i], i+bw/2, plot_data$upper[i])
+    }
+    #Draw points on top of bars
+    points(1:ngroup, plot_data$Predicted, pch=19)
+  }
+
+})
diff --git a/R/posteriorSamples.R b/R/posteriorSamples.R
index 3b84ef0..9389e14 100644
--- a/R/posteriorSamples.R
+++ b/R/posteriorSamples.R
@@ -57,7 +57,7 @@ print.unmarkedPostSamples <- function(x, ...){
 }
 
 setMethod("[", c("unmarkedPostSamples","ANY","ANY","ANY"), 
-          function(x, i, j, k)
+          function(x, i, j, k, drop = FALSE)
 {
-  x@samples[i,j,k]
+  x@samples[i,j,k, drop = drop]
 })
diff --git a/R/predict.R b/R/predict.R
index 68efbb7..37b9f87 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -25,8 +25,7 @@ setMethod("predict", "unmarkedFit",
     orig_formula <- get_formula(object, type)
 
     # 2. If newdata is raster, get newdata from raster as data.frame
-    if(inherits(newdata, c("RasterLayer","RasterStack"))){
-      if(!require(raster)) stop("raster package required", call.=FALSE)
+    if(inherits(newdata, c("RasterLayer","RasterStack","SpatRaster"))){
       is_raster <- TRUE
       orig_raster <- newdata
       newdata <- newdata_from_raster(newdata, all.vars(orig_formula))
@@ -95,8 +94,8 @@ setMethod("check_predict_arguments", "unmarkedFit",
   check_type(object, type)
 
   # Check newdata class
-  if(!inherits(newdata, c("unmarkedFrame", "data.frame", "RasterLayer", "RasterStack"))){
-    stop("newdata must be unmarkedFrame, data.frame, RasterLayer, or RasterStack", call.=FALSE)
+  if(!inherits(newdata, c("unmarkedFrame", "data.frame", "RasterLayer", "RasterStack", "SpatRaster"))){
+    stop("newdata must be unmarkedFrame, data.frame, RasterLayer, RasterStack, or SpatRaster", call.=FALSE)
   }
   invisible(TRUE)
 })
@@ -261,27 +260,39 @@ setMethod("predict_by_chunk", "unmarkedFit",
 # Raster handling functions----------------------------------------------------
 
 # Convert a raster into a data frame to use as newdata
-newdata_from_raster <- function(rst, vars){
-  nd <- raster::as.data.frame(rst)
-  # Handle factor rasters
-  is_fac <- raster::is.factor(rst)
-  rem_string <- paste(paste0("^",names(rst),"_"), collapse="|")
-  names(nd)[is_fac] <- gsub(rem_string, "", names(nd)[is_fac])
+newdata_from_raster <- function(object, vars){
+  if(inherits(object, "Raster")){
+    if(!requireNamespace("raster", quietly=TRUE)) stop("raster package required", call.=FALSE)
+    nd <- raster::as.data.frame(object)
+    # Handle factor rasters
+    is_fac <- raster::is.factor(object)
+    rem_string <- paste(paste0("^",names(object),"_"), collapse="|")
+    names(nd)[is_fac] <- gsub(rem_string, "", names(nd)[is_fac])
+  } else if(inherits(object, "SpatRaster")){
+    if(!requireNamespace("terra", quietly=TRUE)) stop("terra package required", call.=FALSE)
+    nd <- terra::as.data.frame(object)
+  }
   # Check if variables are missing
   no_match <- vars[! vars %in% names(nd)]
   if(length(no_match) > 0){
-    stop(paste0("Variable(s) ",paste(no_match, collapse=", "), " not found in raster stack"),
+    stop(paste0("Variable(s) ",paste(no_match, collapse=", "), " not found in raster(s)"),
          call.=FALSE)
   }
   return(nd)
 }
 
-# Convert predict output into a raster
-raster_from_predict <- function(pr, orig_rst, appendData){
-  new_rast <- data.frame(raster::coordinates(orig_rst), pr)
-  new_rast <- raster::stack(raster::rasterFromXYZ(new_rast))
-  raster::crs(new_rast) <- raster::crs(orig_rst)
-  if(appendData) new_rast <- raster::stack(new_rast, orig_rst)
+raster_from_predict <- function(pr, object, appendData){
+  if(inherits(object, "Raster")){
+    new_rast <- data.frame(raster::coordinates(object), pr)
+    new_rast <- raster::stack(raster::rasterFromXYZ(new_rast))
+    raster::crs(new_rast) <- raster::crs(object)
+    if(appendData) new_rast <- raster::stack(new_rast, object)
+  } else if(inherits(object, "SpatRaster")){
+    new_rast <- data.frame(terra::crds(object), pr)
+    new_rast <- terra::rast(new_rast, type="xyz")
+    terra::crs(new_rast) <- terra::crs(object)
+    if(appendData) new_rast <- c(new_rast, object)
+  }
   new_rast
 }
 
diff --git a/R/ranef.R b/R/ranef.R
index 350ed6d..5f0928a 100644
--- a/R/ranef.R
+++ b/R/ranef.R
@@ -249,24 +249,35 @@ setMethod("ranef", "unmarkedFitMPois",
     lam <- predict(object, type="state")[,1]
     R <- length(lam)
     cp <- getP(object)
-    cp <- cbind(cp, 1-rowSums(cp))
+    cp <- cbind(cp, 1-rowSums(cp, na.rm=TRUE))
     N <- 0:K
-    post <- array(0, c(R, K+1, 1))
+    post <- array(NA, c(R, K+1, 1))
     colnames(post) <- N
     for(i in 1:R) {
         f <- dpois(N, lam[i])
         g <- rep(1, K+1)
-        if(any(is.na(y[i,])) | any(is.na(cp[i,])))
+        yi <- y[i,]
+        cpi <- cp[i,]
+        if(any(is.na(y[i,])) | any(is.na(cp[i,]))){
+          # This only handles cases when all NAs are at the end of the count vector
+          # Not when they are interspersed. I don't know how to handle that situation
+          if(object@data@piFun == "removalPiFun"){
+            warning("NAs in counts and/or covariates for site ",i, ". Keeping only counts before the first NA", call.=FALSE)
+            notNA <- min(which(is.na(y[i,]))[1], which(is.na(cp[i,]))[1], na.rm=TRUE) - 1
+            yi <- yi[c(1:notNA)]
+            cpi <- cpi[c(1:notNA, length(cpi))]
+          } else {
             next
+          }
+        }
         for(k in 1:(K+1)) {
-            yi <- y[i,]
             ydot <- N[k] - sum(yi)
             if(ydot<0) {
                 g[k] <- 0
                 next
             }
-            yi <- c(yi, ydot)
-            g[k] <- g[k] * dmultinom(yi, size=N[k], prob=cp[i,])
+            yik <- c(yi, ydot)
+            g[k] <- g[k] * dmultinom(yik, size=N[k], prob=cpi)
         }
         fudge <- f*g
         post[i,,1] <- fudge / sum(fudge)
diff --git a/R/unmarkedEstimate.R b/R/unmarkedEstimate.R
index 501808e..86caf40 100644
--- a/R/unmarkedEstimate.R
+++ b/R/unmarkedEstimate.R
@@ -149,6 +149,14 @@ setMethod("summary", signature(object = "unmarkedEstimate"),
                    logistic = "logit",
                    cloglog = "cloglog")
     cat(object@name, " (", link, "-scale)", ":\n", sep="")
+
+    has_random <- methods::.hasSlot(object, "randomVarInfo") && 
+      length(object@randomVarInfo) > 0
+    if(has_random){
+      print_randvar_info(object@randomVarInfo)
+      cat("\nFixed effects:\n")
+    }
+
     outDF <- data.frame(Estimate = ests, SE = SEs, z = Z, "P(>|z|)" = p,
                         check.names = FALSE)
     print(outDF, row.names = printRowNames, digits = 3)
diff --git a/R/unmarkedFitList.R b/R/unmarkedFitList.R
index 2974678..67e94ae 100644
--- a/R/unmarkedFitList.R
+++ b/R/unmarkedFitList.R
@@ -42,32 +42,52 @@ setClass("unmarkedFitList",
 
 
 # constructor of unmarkedFitList objects
-fitList <- function(..., fits) {
-    if(length(list(...)) > 0 & !missing(fits))
-        stop("Do not use both the '...' and 'fits' arguments")
-    if(missing(fits)) {
-        fits <- list(...)
-        isList <- sapply(fits, function(x) is.list(x))
-        if(sum(isList) > 1)
-            stop("Specify models as common-seperated objects, or use fits = 'mylist'")
-        if(isList[1L]) {
-            warning("If supplying a list of fits, use fits = 'mylist'")
-            fits <- fits[[1L]] 	# This is allowed for back-compatability.
-            }
-        if(is.null(names(fits))) {
-            c <- match.call(expand.dots = FALSE)
-            names(fits) <- as.character(c[[2]])
-            warning("Your list was unnamed, so model names were added as object names")
-            }
-        }
+fitList <- function(..., fits, autoNames=c("object","formula")) {
+  autoNames <- match.arg(autoNames)
+  if(length(list(...)) > 0 & !missing(fits))
+    stop("Do not use both the '...' and 'fits' arguments")
+  if(missing(fits)) {
+    fits <- list(...)
+    isList <- sapply(fits, function(x) is.list(x))
+    if(sum(isList) > 1)
+      stop("Specify models as common-seperated objects, or use fits = 'mylist'")
+    if(isList[1L]) {
+      warning("If supplying a list of fits, use fits = 'mylist'")
+      fits <- fits[[1L]] 	# This is allowed for back-compatability.
+    }
     if(is.null(names(fits))) {
-        names(fits) <- as.character(1:(length(fits)))
-        warning("Your list was unnamed, so model names were added as c('1','2',...)")
+      message("Your list was unnamed, so model names were added automatically")
+      if(autoNames=="formula"){
+        if(inherits(fits[[1]], c("unmarkedFitOccuMulti", "unmarkedFitOccuMS"))){
+          warning("simple formula not available, naming based on object instead", call.=FALSE)
+          autoNames <- "object"
+        } else {
+          names(fits) <- lapply(fits, function(x) gsub(" ", "", as.character(deparse(x@formula))))
         }
-    umfl <- new("unmarkedFitList", fits=fits)
-    return(umfl)
+      }
+      if(autoNames=="object"){
+        c <- match.call(expand.dots = FALSE)
+        names(fits) <- as.character(c[[2]])
+      }
     }
-
+  }
+  if(is.null(names(fits))) {
+    message("Your list was unnamed, so model names were added automatically")
+    if(autoNames=="formula"){
+        if(inherits(fits[[1]], c("unmarkedFitOccuMulti", "unmarkedFitOccuMS"))){
+        warning("simple formula not available, naming as numbers instead", call.=FALSE)
+        autoNames <- "object"
+      } else {
+        names(fits) <- lapply(fits, function(x) gsub(" ", "", as.character(deparse(x@formula))))
+      }
+    }
+    if(autoNames=="object"){
+      names(fits) <- as.character(1:(length(fits)))
+    }
+  }
+  umfl <- new("unmarkedFitList", fits=fits)
+  return(umfl)
+}
 
 setMethod("summary", "unmarkedFitList", function(object) {
     fits <- object@fits
diff --git a/R/unmarkedFrame.R b/R/unmarkedFrame.R
index b6a9273..3b8895f 100644
--- a/R/unmarkedFrame.R
+++ b/R/unmarkedFrame.R
@@ -161,7 +161,10 @@ setClass("unmarkedFrameDSO",
 #Convert covs provided as list of matrices/dfs to data frame
 covsToDF <- function(covs, name, obsNum, numSites){
   if(!inherits(covs, "list")) return(covs)
-
+  
+  if(is.null(names(covs)) | any(is.na(names(covs))) | any(names(covs)=="")){
+    stop("All elements of list provided to ", name, " argument must be named", call.=FALSE)
+  }
   lapply(covs, function(x){
     if(!inherits(x, c("matrix", "data.frame")))
       stop(paste("At least one element of", name, "is not a matrix or data frame."))
diff --git a/man/fitList.Rd b/man/fitList.Rd
index e5047bd..e01e315 100644
--- a/man/fitList.Rd
+++ b/man/fitList.Rd
@@ -1,11 +1,12 @@
 \name{fitList}
 \alias{fitList}
 \title{constructor of unmarkedFitList objects}
-\usage{fitList(..., fits)}
+\usage{fitList(..., fits, autoNames=c("object", "formula"))}
 \description{Organize models for model selection or model-averaged prediction.}
 \arguments{
 \item{...}{Fitted models. Preferrably named.}
 \item{fits}{An alternative way of providing the models. A (preferrably named) list of fitted models.}
+\item{autoNames}{Option to change the names \code{unmarked} assigns to models if you don't name them yourself. If \code{autoNames="object"}, models in the \code{fitList} will be named based on their R object names.  If \code{autoNames="formula"}, the models will instead be named based on their formulas. This is not possible for some model types.}
 }
 \note{Two requirements exist to conduct AIC-based model-selection and model-averaging in unmarked. First, the data objects (ie, unmarkedFrames) must be identical among fitted models. Second, the response matrix must be identical among fitted models after missing values have been removed. This means that if a response value was removed in one model due to missingness, it needs to be removed from all models.
 } 
diff --git a/man/plotEffects.Rd b/man/plotEffects.Rd
new file mode 100644
index 0000000..d176d43
--- /dev/null
+++ b/man/plotEffects.Rd
@@ -0,0 +1,100 @@
+\name{plotEffects}
+\alias{plotEffects}
+\alias{plotEffects-methods}
+\alias{plotEffects,unmarkedFit-method}
+\alias{plotEffectsData}
+\alias{plotEffectsData-methods}
+\alias{plotEffectsData,unmarkedFit-method}
+
+\title{Plot marginal effects of covariates in unmarked models}
+
+\description{This function generates a plot visualizing the effects
+ of a single covariate on a parameter (e.g. occupancy, abundance) in an unmarked 
+ model. If the covariate is numeric, the result is a line plot with an error
+ ribbon where the x-axis is the range of the covariate and the y-axis is the 
+ predicted parameter value. If the covariate is an R factor (i.e., categorical),
+ the x-axis instead contains each unique value of the covariate. 
+
+ All covariates in the model besides the one being plotted are held either at
+ their median value (if they are numeric) or at their reference level (if they
+ are factors).
+
+ Some types of unmarked models may require additional arguments, which are passed
+ to the matching \code{predict} method. For example, \code{unmarkedFitOccuMulti} 
+ models require the \code{species} argument to be included in the function 
+ call in order to work properly.
+
+ If you want to customize a plot, the easiest approach is to get data 
+ formatted for plotting using \code{plotEffectsData}, and use that. If
+ you want to see and/or modify the code used by \code{plotEffects} to generate
+ the default plots, run \code{getMethod("plotEffects", "unmarkedFit")} in
+ the R console.
+}
+
+\usage{
+\S4method{plotEffects}{unmarkedFit}(object, type, covariate, level=0.95, ...)
+\S4method{plotEffectsData}{unmarkedFit}(object, type, covariate, level=0.95, ...)
+}
+
+\arguments{
+  \item{object}{A fitted model inheriting class \code{unmarkedFit}}
+  \item{type}{Submodel in which the covariate of interest can be found, for
+    example \code{"state"} or \code{"det"}. This will depend on the fitted model}
+  \item{covariate}{The name of the covariate to be plotted, as a character string}
+  \item{level}{Confidence level for the error ribbons or bars}
+  \item{...}{Other arguments passed to the \code{predict} function, required
+    for some \code{unmarkedFit} types such as \code{unmarkedFitOccuMulti}}
+} 
+
+\value{A plot (\code{plotEffects} or a data frame (\code{plotEffectsData})
+  containing values to be used in a plot.
+}
+
+\author{Ken Kellner \email{contact@kenkellner.com}}
+
+\examples{
+
+\dontrun{
+
+# Simulate data and build an unmarked frame
+set.seed(123)
+dat_occ <- data.frame(x1=rnorm(500))
+dat_p <- data.frame(x2=rnorm(500*5))
+
+y <- matrix(NA, 500, 5)
+z <- rep(NA, 500)
+
+b <- c(0.4, -0.5, 0.3, 0.5)
+
+re_fac <- factor(sample(letters[1:5], 500, replace=T))
+dat_occ$group <- re_fac
+re <- rnorm(5, 0, 1.2)
+re_idx <- as.numeric(re_fac)
+
+idx <- 1
+for (i in 1:500){
+  z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
+  for (j in 1:5){
+    y[i,j] <- z[i]*rbinom(1,1,
+                    plogis(b[3] + b[4]*dat_p$x2[idx]))
+    idx <- idx + 1
+  }
+}
+
+umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
+
+# Fit model
+(fm <- occu(~x2 ~x1 + group, umf))
+
+# Plot marginal effects of various covariates
+plotEffects(fm, "state", "x1")
+plotEffects(fm, "state", "group")
+plotEffects(fm, "det", "x2")
+
+# Get raw data used for a plot
+plotEffectsData(fm, "state", "group")
+
+# See code used by plotEffects so you can edit it yourself and customize the plot
+methods::getMethod("plotEffects", "unmarkedFit")
+}
+}
diff --git a/tests/testthat/test_fitList.R b/tests/testthat/test_fitList.R
index c33c03d..444a3a0 100644
--- a/tests/testthat/test_fitList.R
+++ b/tests/testthat/test_fitList.R
@@ -37,5 +37,9 @@ test_that("fitList operations work",{
 
   se <- SE(mt)
   expect_equal(dim(se), c(2,5))
-
+  
+  fl <- expect_message(fitList(fm, fm2, autoNames='formula'))
+  expect_equal(names(fl@fits), c("~o1+o2~x", "~1~x"))
+  fl <- expect_message(fitList(fits=list(fm, fm2), autoNames='formula'))
+  expect_equal(names(fl@fits), c("~o1+o2~x", "~1~x"))
 })
diff --git a/tests/testthat/test_gdistremoval.R b/tests/testthat/test_gdistremoval.R
index 8f453a5..fdc2ca1 100644
--- a/tests/testthat/test_gdistremoval.R
+++ b/tests/testthat/test_gdistremoval.R
@@ -343,6 +343,12 @@ test_that("gdistremoval can fit models",{
   expect_is(fl, "unmarkedFitList")
   ms <- modSel(fl)
   expect_is(ms, "unmarkedModSel")
+
+  #plotEffects
+  pdf(NULL)
+  plotEffects(fit, 'lambda', "sc1")
+  plotEffects(fit, 'rem', 'oc1')
+  dev.off()
 })
 
 test_that("gdistremoval predict method works",{
@@ -405,9 +411,21 @@ test_that("gdistremoval handles NAs",{
   fit <- gdistremoval(~1,removalformula=~1,distanceformula=~1, data=umf)
   expect_equivalent(coef(fit), c(2.0675,3.908,-2.1433), tol=1e-3)
 
+  # Can't have missing site covs
   umf2 <- umf
   umf2@siteCovs$sc1[1] <- NA
   expect_error(gdistremoval(~sc1,removalformula=~1,distanceformula=~1, data=umf2))
+  
+  # This errors because missing obs cov does not match missing removal data
+  umf2 <- umf
+  umf2@obsCovs$oc1[1] <- NA
+  expect_error(gdistremoval(~1,removalformula=~oc1,distanceformula=~1, data=umf2))
+
+  # This does not error because missing obs cov matches missing removal data
+  umf2 <- umf
+  umf2@obsCovs$oc1[6] <- NA
+  fit <- gdistremoval(~1,removalformula=~oc1,distanceformula=~1, data=umf2)
+  expect_true(is.na(predict(fit, 'rem')$Predicted[6]))
 })
 
 test_that("multi-period data works with gdistremoval",{
@@ -438,8 +456,40 @@ test_that("multi-period data works with gdistremoval",{
 
   # ranef
   r <- ranef(fit)
-  expect_equal(dim(bup(r)), c(30,5))
+  expect_equal(dim(r@post), c(30, 44, 1))
+  expect_equal(length(bup(r)), 30)
 
+  # fitted
+  ft <- fitted(fit)
+  expect_equal(dim(ft$dist), dim(fit@data@yDistance))
+
+  # Entire missing secondary period
+  umf2 <- umf
+  # remove 2nd period at first site
+  umf2@yDistance[1,5:8] <- NA
+  umf2@yRemoval[1, 6:10] <- NA
+  umf2@obsCovs$oc1[6:10] <- NA
+  umf2@yearlySiteCovs$ysc1[2] <- NA
+  
+
+  fit2 <- gdistremoval(~sc1,phiformula=~ysc1, removalformula=~oc1,distanceformula=~1, data=umf2)
+
+  gp <- getP(fit2)
+  expect_true(is.na(gp$phi[1,2]))
+
+  pr <- predict(fit2, 'rem', level=NULL)
+  expect_true(all(is.na(pr$Predicted[6:10])))
+
+  r2 <- ranef(fit2)
+  expect_true(!any(is.na(r2@post)))
+
+  s <- simulate(fit2)
+  expect_true(all(is.na(s[[1]]$yDistance[1,5:8])))
+  
+  res <- residuals(fit2)
+  expect_true(all(is.na(res$rem[1,6:10])))
+  pb <- parboot(fit2, nsim=2)
+  expect_is(pb, "parboot")
 })
 
 test_that("gdistremoval works with random effects",{
diff --git a/tests/testthat/test_modSel.R b/tests/testthat/test_modSel.R
index 14fffc2..b6d1485 100644
--- a/tests/testthat/test_modSel.R
+++ b/tests/testthat/test_modSel.R
@@ -12,11 +12,11 @@ test_that("fitLists can be constructed",{
 
   fits1.1 <- fitList(m1=fm1, m2=fm2)
   expect_equal(names(fits1.1@fits), c("m1","m2"))
-  expect_warning(fits1.2 <- fitList(fm1, fm2))
+  expect_message(fits1.2 <- fitList(fm1, fm2))
   expect_equal(names(fits1.2@fits), c("fm1","fm2"))
   fits2.1 <- fitList(fits = list(m1=fm1, m2=fm2))
   expect_equal(names(fits2.1@fits), c("m1","m2"))
-  expect_warning(fits2.2 <- fitList(fits = list(fm1, fm2)))
+  expect_message(fits2.2 <- fitList(fits = list(fm1, fm2)))
   expect_equal(names(fits2.2@fits), c("1","2"))
 
   expect_equal(fits1.1, fits2.1)
diff --git a/tests/testthat/test_multinomPois.R b/tests/testthat/test_multinomPois.R
index 6b4b693..1103839 100644
--- a/tests/testthat/test_multinomPois.R
+++ b/tests/testthat/test_multinomPois.R
@@ -101,9 +101,14 @@ test_that("multinomPois can fit a removal model",{
     res <- residuals(m2_C)
     expect_equal(dim(res), dim(umf1@y))
 
-    expect_warning(r <- ranef(m2_C))
-    expect_equal(dim(r@post), c(3,56,1))
-    expect_equal(bup(r), c(10.794,0.000,2.655), tol=1e-4)
+    expect_warning(r <- ranef(m2_C, K=50))
+    expect_equal(dim(r@post), c(3,51,1))
+    expect_equal(bup(r), c(10.794,6.9317,2.655), tol=1e-4)
+
+    umf2 <- unmarkedFrameMPois(y=y, siteCovs=data.frame(x1=rnorm(5)), type="removal")
+    m4 <- multinomPois(~1~x1, umf2)
+    r <- ranef(m4, K=30)
+    expect_equal(dim(r@post), c(5,31,1))
 
     expect_warning(s <- simulate(m2_C, 2, na.rm=FALSE))
     expect_equal(length(s), 2)
diff --git a/tests/testthat/test_occu.R b/tests/testthat/test_occu.R
index 7ea3afe..fdcede8 100644
--- a/tests/testthat/test_occu.R
+++ b/tests/testthat/test_occu.R
@@ -126,7 +126,11 @@ test_that("occu can fit models with covariates",{
   stMat <- fm@estimates@estimates$state@covMat
   expect_equivalent(detMat, matrix(rep(NA,9),nrow=3))
   expect_equivalent(stMat, matrix(rep(NA,4),nrow=2))
-
+  
+  # Trap attempts to use a variable in formula that isn't in covariates
+  fake <- rnorm(3)
+  expect_error(fm <- occu(~ o1 + o2 ~ fake, data = umf))
+  expect_error(fm <- occu(~ o1 + fake ~ o1, data = umf))
 })
 
 test_that("occu handles NAs",{
diff --git a/tests/testthat/test_occuFP.R b/tests/testthat/test_occuFP.R
index f1381bb..fc61461 100644
--- a/tests/testthat/test_occuFP.R
+++ b/tests/testthat/test_occuFP.R
@@ -22,7 +22,7 @@ test_that("occuFP model can be fit",{
   m1 <- occuFP(detformula = ~ METH, FPformula = ~1,
                stateformula = ~ habitat, data = umf1)
   expect_equal(names(m1), c("state","det","fp"))
-  expect_warning(fl <- fitList(m1,m1))
+  expect_message(fl <- fitList(m1,m1))
   expect_is(fl,"unmarkedFitList")
   expect_equal(length(fl@fits), 2)
 
diff --git a/tests/testthat/test_occuMS.R b/tests/testthat/test_occuMS.R
index aa0d024..7c25acf 100644
--- a/tests/testthat/test_occuMS.R
+++ b/tests/testthat/test_occuMS.R
@@ -227,7 +227,9 @@ test_that("occuMS can fit the multinomial model",{
   expect_equivalent(r@post[1,,1], c(0,0.5222,0.4778), tol=1e-4)
 
   #Check fitList
-  expect_warning(fl <- fitList(fit_C, fit_C))
+  expect_message(fl <- fitList(fit_C, fit_C))
+  expect_message(expect_warning(fl <- fitList(fit_C, fit_C, autoNames='formula')))
+  expect_message(expect_warning(fl <- fitList(fits=list(fit_C, fit_C), autoNames='formula')))
   expect_is(fl,"unmarkedFitList")
   expect_equivalent(length(fl@fits), 2)
 
diff --git a/tests/testthat/test_occuMulti.R b/tests/testthat/test_occuMulti.R
index fbc264a..d37cbc5 100644
--- a/tests/testthat/test_occuMulti.R
+++ b/tests/testthat/test_occuMulti.R
@@ -42,7 +42,9 @@ test_that("occuMulti can fit simple models",{
   expect_equivalent(det, rep(1,length(detlist)), tolerance= 1e-4)
 
   #Check fitList
-  expect_warning(fl <- fitList(fm, fm))
+  expect_message(fl <- fitList(fm, fm))
+  expect_message(expect_warning(fl <- fitList(fm, fm, autoNames='formula')))
+  expect_message(expect_warning(fl <- fitList(fits=list(fm, fm), autoNames='formula')))
   expect_is(fl,"unmarkedFitList")
   expect_equivalent(length(fl@fits), 2)
 
@@ -335,7 +337,7 @@ test_that("occuMulti predict method works",{
 
   #fitList with maxOrder set
   fm2 <- occuMulti(c("~1","~1"), c("~1","~1"), umf, maxOrder=1)
-  expect_warning(fl2 <- fitList(fm, fm2))
+  expect_message(fl2 <- fitList(fm, fm2))
   expect_is(fl2, "unmarkedFitList")
   ms <- modSel(fl2)
   expect_is(ms, "unmarkedModSel")
diff --git a/tests/testthat/test_occuTTD.R b/tests/testthat/test_occuTTD.R
index f152cde..0fb7d81 100644
--- a/tests/testthat/test_occuTTD.R
+++ b/tests/testthat/test_occuTTD.R
@@ -382,6 +382,12 @@ test_that("occuTTD can fit a dynamic model",{
   expect_equivalent(dim(s[[1]]), c(100,4))
   r <- residuals(fit)
   expect_equivalent(dim(r), c(100,4))
+  pdf(NULL)
+  plotEffects(fit, 'psi', 'elev')
+  plotEffects(fit, 'col', 'forest')
+  plotEffects(fit, 'det', 'wind')
+  plotEffects(fit, 'det', 'obs')
+  dev.off()
 
   #Check ranef
   r <- ranef(fit)
diff --git a/tests/testthat/test_plotEffects.R b/tests/testthat/test_plotEffects.R
new file mode 100644
index 0000000..3aea3ec
--- /dev/null
+++ b/tests/testthat/test_plotEffects.R
@@ -0,0 +1,48 @@
+context("plotEffects")
+
+skip_on_cran()
+
+set.seed(123)
+dat_occ <- data.frame(x1=rnorm(500))
+dat_p <- data.frame(x2=rnorm(500*5))
+
+y <- matrix(NA, 500, 5)
+z <- rep(NA, 500)
+
+b <- c(0.4, -0.5, 0.3, 0.5)
+
+re_fac <- factor(sample(letters[1:5], 500, replace=T))
+dat_occ$group <- re_fac
+re <- rnorm(5, 0, 1.2)
+re_idx <- as.numeric(re_fac)
+
+idx <- 1
+for (i in 1:500){
+  z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
+  for (j in 1:5){
+    y[i,j] <- z[i]*rbinom(1,1,
+                    plogis(b[3] + b[4]*dat_p$x2[idx]))
+    idx <- idx + 1
+  }
+}
+
+umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
+
+fm <- occu(~x2 ~x1 + group, umf)
+
+test_that("plotEffects works", {
+  
+  pdf(NULL)
+  plotEffects(fm, "state", "x1")
+  plotEffects(fm, "state", "group")
+  plotEffects(fm, "det", "x2")
+  dev.off()
+
+  expect_error(plotEffects(fm, "state", "x2"))
+
+  dat <- plotEffectsData(fm, "state", "group")
+
+  expect_true(inherits(dat, "data.frame"))
+  expect_equal(nrow(dat), 5)
+
+})
diff --git a/tests/testthat/test_predict.R b/tests/testthat/test_predict.R
index e92add1..686da9b 100644
--- a/tests/testthat/test_predict.R
+++ b/tests/testthat/test_predict.R
@@ -138,3 +138,41 @@ test_that("predicting from raster works",{
   nd_2 <- nd_raster[[1]]
   expect_error(predict(mod, 'state', newdata=nd_2))
 })
+
+test_that("predicting from terra::rast works",{
+
+  skip_if(!require(terra), "terra package unavailable")
+
+  set.seed(123)
+  # Create rasters
+  # Elevation
+  r_elev <- data.frame(x=rep(1:10, 10), y=rep(1:10, each=10), z=rnorm(100))
+  r_elev <- terra::rast(r_elev, type="xyz")
+
+  #Group
+  r_group <- data.frame(x=rep(1:10, 10), y=rep(1:10, each=10),
+                      z=sample(1:length(levels(umf@siteCovs$group)), 100, replace=T))
+  # Convert to 'factor' raster
+  r_group <- terra::as.factor(terra::rast(r_group, type="xyz"))
+  levels(r_group) <- data.frame(ID=terra::levels(r_group)[[1]]$ID, group=levels(umf@siteCovs$group))
+
+  # Stack
+  nd_raster <- c(r_elev, r_group)
+  names(nd_raster) <- c("elev", "group")
+  terra::crs(nd_raster) <- "epsg:32616"
+
+  pr <- predict(mod, 'state', newdata=nd_raster)
+  expect_is(pr, 'SpatRaster')
+  expect_equal(names(pr), c("Predicted","SE","lower","upper"))
+  expect_equivalent(pr[1,1][1], 0.3675313, tol=1e-5)
+  expect_equal(crs(pr), crs(nd_raster))
+
+  #append data
+  pr <- predict(mod, 'state', newdata=nd_raster, appendData=TRUE)
+  expect_is(pr, 'SpatRaster')
+  expect_equal(names(pr)[5:6], c("elev","group"))
+
+  # Missing levels are handled
+  nd_2 <- nd_raster[[1]]
+  expect_error(predict(mod, 'state', newdata=nd_2))
+})
diff --git a/tests/testthat/test_ranef_predict.R b/tests/testthat/test_ranef_predict.R
index 251ed29..7da3f54 100644
--- a/tests/testthat/test_ranef_predict.R
+++ b/tests/testthat/test_ranef_predict.R
@@ -29,7 +29,7 @@ test_that("ranef predict method works",{
   expect_equivalent(dim(ps@samples), c(9,1,10))
 
   # Brackets
-  expect_equal(ps[1,1,1], ps@samples[1,1,1])
+  expect_equal(ps[1,1,1], ps@samples[1,1,1,drop=FALSE])
 
   # Method for unmarkedFit objects
   set.seed(123)
diff --git a/tests/testthat/test_unmarkedFrame.R b/tests/testthat/test_unmarkedFrame.R
index b5c6044..c3194fe 100644
--- a/tests/testthat/test_unmarkedFrame.R
+++ b/tests/testthat/test_unmarkedFrame.R
@@ -214,4 +214,45 @@ test_that("yearlySiteCovs processing works",{
   expect_equivalent(dim(as_df2),c(50,25))
 })
 
+test_that("lists provided to obsCovs or yearlySiteCovs must be named", {
+    y <- matrix(1:27, 3)
+    sc <- data.frame(x1 = 1:3)
+    ysc <- list(x2 = matrix(1:9, 3))
+    oc <- list(x3 = matrix(1:27, 3))
+
+    umf1 <- unmarkedMultFrame(
+        y = y,
+        siteCovs = sc,
+        yearlySiteCovs = ysc,
+        obsCovs = oc,
+        numPrimary = 3)
+    expect_is(umf1, "unmarkedMultFrame")
 
+
+    oc <- list(matrix(1:27, 3))
+    umf1 <- expect_error(unmarkedMultFrame(
+        y = y,
+        siteCovs = sc,
+        yearlySiteCovs = ysc,
+        obsCovs = oc,
+        numPrimary = 3))
+
+    oc <- list(x3 = matrix(1:27, 3))
+    ysc <- list(matrix(1:9, 3))
+    umf1 <- expect_error(unmarkedMultFrame(
+        y = y,
+        siteCovs = sc,
+        yearlySiteCovs = ysc,
+        obsCovs = oc,
+        numPrimary = 3))
+
+    ysc <- list(x2 = matrix(1:9, 3))
+    oc <- list(x3 = matrix(1:27, 3), matrix(1:27, 3))
+
+    umf1 <- expect_error(unmarkedMultFrame(
+        y = y,
+        siteCovs = sc,
+        yearlySiteCovs = ysc,
+        obsCovs = oc,
+        numPrimary = 3))
+})
diff --git a/vignettes/occuMulti.Rmd b/vignettes/occuMulti.Rmd
index 25eae85..59f602e 100644
--- a/vignettes/occuMulti.Rmd
+++ b/vignettes/occuMulti.Rmd
@@ -64,7 +64,7 @@ The data are a simplified version of the data used in @Rota2016, with the data c
 The dataset is included with `unmarked` and is called `MesoCarnivores`.
 First, we need to load in the dataset, which is a list with several components.
 
-```{r, eval=FALSE}
+```{r}
 library(unmarked)
 data(MesoCarnivores)
 names(MesoCarnivores)