Add occuCOP - debug package build

10 files changed, 1472 insertions, 5 deletions

diff --git a/.Rbuildignore b/.Rbuildignore
index 70e8958..1ac8847 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -12,3 +12,5 @@ README.Rmd
 ^\.github$
 ^_pkgdown\.yml$
 ^vignettes/colext.Rmd.orig
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
index fe18f2a..7640862 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,3 +28,6 @@ bc
 .Rhistory
 symbols.rds
 
+# Rproj
+.Rproj.user
+*.Rproj
diff --git a/DESCRIPTION b/DESCRIPTION
index 5658925..762bbdf 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -54,6 +54,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R'
     'goccu.R'
     'RcppExports.R'
     'zzz.R'
+    'occuCOP.R'
 LinkingTo: 
     Rcpp,
     RcppArmadillo,
diff --git a/NAMESPACE b/NAMESPACE
index 3234918..c3c80af 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,7 +24,7 @@ importFrom(pbapply, pbsapply, pblapply)
 export(occu, occuFP, occuRN, pcount, pcountOpen, multinomPois, distsamp,
        colext, gmultmix, gdistsamp, gpcount, occuPEN, occuPEN_CV, occuMulti,
        occuMS, computeMPLElambda, pcount.spHDS, occuTTD, distsampOpen,
-       multmixOpen, nmixTTD, gdistremoval, goccu)
+       multmixOpen, nmixTTD, gdistremoval, goccu, occuCOP)
 export(removalPiFun, doublePiFun)
 export(makeRemPiFun, makeCrPiFun, makeCrPiFunMb, makeCrPiFunMh)
 
@@ -40,7 +40,8 @@ exportClasses(unmarkedFit, unmarkedFitOccu, unmarkedFitOccuFP, unmarkedFitDS,
               unmarkedFrameGPC, unmarkedEstimate, unmarkedFitList,
               unmarkedModSel, unmarkedRanef, unmarkedFrameOccuMulti,
               unmarkedFrameOccuMS, unmarkedFrameGDR, unmarkedCrossVal,
-              unmarkedPostSamples, unmarkedPower, unmarkedPowerList)
+              unmarkedPostSamples, unmarkedPower, unmarkedPowerList,
+              unmarkedFrameCOP, unmarkedFitCOP)
 
 # Methods
 exportMethods(backTransform, coef, confint, coordinates, fitted, getData,
@@ -61,7 +62,8 @@ export(unmarkedEstimate, fitList, mapInfo, unmarkedFrame,
     unmarkedFrameDS, unmarkedMultFrame, unmarkedFrameGMM,
     unmarkedFramePCO, unmarkedFrameGDS, unmarkedFrameGPC, unmarkedFrameOccuMulti,
     unmarkedFrameOccuMS, unmarkedFrameOccuTTD, unmarkedFrameDSO, 
-    unmarkedFrameMMO, unmarkedFrameGDR, unmarkedFrameGOccu)
+    unmarkedFrameMMO, unmarkedFrameGDR, unmarkedFrameGOccu,
+    unmarkedFrameCOP)
 
 # Formatting
 export(csvToUMF, formatLong, formatWide, formatMult, formatDistData)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 79ab967..70b57ac 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -49,6 +49,10 @@ nll_occu <- function(y, X, V, beta_psi, beta_p, nd, knownOcc, navec, X_offset, V
     .Call(`_unmarked_nll_occu`, y, X, V, beta_psi, beta_p, nd, knownOcc, navec, X_offset, V_offset, link_psi)
 }
 
+nll_occuCOP <- function(y, L, Xpsi, Xlambda, beta_psi, beta_lambda, removed) {
+    .Call(`_unmarked_nll_occuCOP`, y, L, Xpsi, Xlambda, beta_psi, beta_lambda, removed)
+}
+
 nll_occuMS <- function(beta, y, dm_state, dm_phi, dm_det, sind, pind, dind, prm, S, T, J, N, naflag, guide) {
     .Call(`_unmarked_nll_occuMS`, beta, y, dm_state, dm_phi, dm_det, sind, pind, dind, prm, S, T, J, N, naflag, guide)
 }
diff --git a/R/boot.R b/R/boot.R
index 2f8cffc..0e6505c 100644
--- a/R/boot.R
+++ b/R/boot.R
@@ -69,8 +69,13 @@ setMethod("parboot", "unmarkedFit",  function(object, statistic=SSE, nsim=10,
     simdata <- replaceY(object@data, x)
     tryCatch({
       #if(runif(1,0,1) < 0.5) stop("fail") # for testing error trapping
-      fit <- update(object, data=simdata, starts=starts, se=FALSE)
-      statistic(fit, ...)
+      if (class(object) == "unmarkedFitCOP") {
+        fit <- suppressMessages(update(object, data = simdata, se = FALSE))
+        statistic(fit, ...)
+      } else {
+        fit <- update(object, data = simdata, starts = starts, se = FALSE)
+        statistic(fit, ...)
+      }
     }, error=function(e){
       t0[] <- NA
       t0
diff --git a/R/occuCOP.R b/R/occuCOP.R
new file mode 100644
index 0000000..4f3c0f4
--- /dev/null
+++ b/R/occuCOP.R
@@ -0,0 +1,933 @@
+# Classes ----------------------------------------------------------------------
+
+## unmarkedFrameCOP class ----
+setClass(
+  "unmarkedFrameCOP",
+  representation(L = "matrix"),
+  contains = "unmarkedFrame",
+  validity = function(object) {
+    errors <- character(0)
+    M <- nrow(object@y)
+    J <- ncol(object@y)
+    y_integers = sapply(object@y, check.integer, na.ignore = T)
+    if (!all(y_integers)) {
+      errors <- c(errors,
+                  paste(
+                    "Count detection should be integers. Non-integer values:",
+                    paste(object@y[which(!y_integers)], collapse = ', ')
+                  ))
+    }
+    if (!all(all(dim(object@L) == dim(object@y)))){
+      errors <- c( errors, paste(
+        "L should be a matrix of the same dimension as y, with M =", M,
+        "rows (sites) and J =", J, "columns (sampling occasions)."
+      ))}
+    if (length(errors) == 0) TRUE
+    else errors
+  }
+)
+
+
+## unmarkedFitCOP class ----
+setClass("unmarkedFitCOP",
+         representation(removed_obs = "matrix",
+                        formlist = "list",
+                        nll = "optionalNumeric",
+                        convergence="optionalNumeric"),
+         contains = "unmarkedFit")
+
+
+
+# Methods ----------------------------------------------------------------------
+
+## getDesign method ----
+setMethod(
+  "getDesign", "unmarkedFrameCOP",
+  function(umf, formlist, na.rm = TRUE) {
+    
+    "
+    getDesign convert the information in the unmarkedFrame to a format
+    usable by the likelihood function:
+    - It creates model design matrices for fixed effects (X) for each parameter (here lambda, psi) 
+    - It creates model design matrices for random effects (Z) for each parameter (here lambda, psi)
+    - It handles missing data
+    "
+    
+    # Retrieve useful informations from umf
+    M <- numSites(umf)
+    J <- obsNum(umf)
+    y <- getY(umf)
+    L <- getL(umf)
+    
+    # Occupancy submodel -------------------------------------------------------
+    
+    # Retrieve the fixed-effects part of the formula
+    psiformula <- lme4::nobars(as.formula(formlist$psiformula))
+    psiVars <- all.vars(psiformula)
+    
+    # Retrieve the site covariates
+    sc <- siteCovs(umf)
+    if(is.null(sc)) {
+      sc <- data.frame(.dummy = rep(0, M))
+    }
+    
+    # Check for missing variables
+    psiMissingVars <- psiVars[!(psiVars %in% names(sc))]
+    if (length(psiMissingVars) > 0) {
+      stop(paste0(
+        "Variable(s) '",
+        paste(psiMissingVars, collapse = "', '"),
+        "' not found in siteCovs"
+      ), call. = FALSE)
+    }
+    
+    # State model matrix for fixed effects
+    Xpsi <- model.matrix(
+      psiformula,
+      model.frame(psiformula, sc, na.action = NULL)
+    )
+    # State model matrix for random effects
+    Zpsi <- get_Z(formlist$psiformula, sc)
+    
+    # Detection submodel -------------------------------------------------------
+    
+    # Retrieve the fixed-effects part of the formula
+    lambdaformula <- lme4::nobars(as.formula(formlist$lambdaformula))
+    lambdaVars <- all.vars(lambdaformula)
+    
+    # Retrieve the observation covariates
+    oc <- obsCovs(umf)
+    if(is.null(oc)) {
+      oc <- data.frame(.dummy = rep(0, M*J))
+    }
+    
+    # Check for missing variables
+    lambdaMissingVars <- lambdaVars[!(lambdaVars %in% names(oc))]
+    if (length(lambdaMissingVars) > 0) {
+      stop(paste(
+        "Variable(s)",
+        paste(lambdaMissingVars, collapse = ", "),
+        "not found in obsCovs"
+      ), call. = FALSE)
+    }
+    
+    # Detection model matrix for fixed effects
+    Xlambda <- model.matrix(
+      lambdaformula,
+      model.frame(lambdaformula, oc, na.action = NULL)
+    )
+    # Detection model matrix for random effects
+    Zlambda <- get_Z(formlist$lambdaformula, oc)
+    
+    # Missing data -------------------------------------------------------------
+    
+    # Missing count data
+    missing_y <- is.na(y)
+    
+    # Missing site covariates
+    # (TRUE if at least one site covariate is missing in a site)
+    missing_sc <- apply(Xpsi, 1, function(x) any(is.na(x)))
+    
+    # Missing observation covariates
+    # (TRUE if at least one observation covariate is missing in a sampling occasion in a site)
+    missing_oc <- apply(Xlambda, 1, function(x) any(is.na(x)))
+    
+    # Matrix MxJ of values to not use because there is some data missing
+    removed_obs <- 
+      # If there is count data missing in site i and obs j
+      missing_y | 
+      # If there is any site covariate missing in site i
+      replicate(n = J, missing_sc) |
+      # If there is any observation covariate missing in site i and obs j
+      matrix(missing_oc, M, J, byrow = T)
+    
+    if (any(removed_obs)) {
+      if (na.rm) {
+        nb_missing_sites <- sum(rowSums(!removed_obs) == 0)
+        nb_missing_observations <- sum(is.na(removed_obs))
+        warning("There is missing data: ",
+                sum(missing_y), " missing count data, ",
+                sum(missing_sc), " missing site covariate(s), ",
+                sum(missing_oc), " missing observation covariate(s). ",
+                "Data from only ", (M*J)-sum(removed_obs), " observations out of ", (M*J), " are used, ",
+                "from ", M-nb_missing_sites, " sites out of ", M, ".\n\t"
+        )
+      } else {
+        stop("na.rm=FALSE and there is missing data :\n\t",
+             sum(missing_y), " missing count data (y)\n\t",
+             sum(missing_sc), " missing site covariates (siteCovs)\n\t",
+             sum(missing_oc), " missing observation covariates (obsCovs)")
+      }
+    }
+    
+    # Output -------------------------------------------------------------------
+    return(list(
+      y = y,
+      Xpsi = Xpsi,
+      Zpsi = Zpsi,
+      Xlambda = Xlambda,
+      Zlambda = Zlambda,
+      removed_obs = removed_obs
+    ))
+  })
+
+
+## getL method ----
+setGeneric("getL", function(object) standardGeneric("getL"))
+setMethod("getL", "unmarkedFrameCOP", function(object) {
+  return(object@L)
+})
+
+
+## show method ----
+setMethod("show", "unmarkedFrameCOP", function(object) {
+  J <- ncol(object@L)
+  df_unmarkedFrame <- as(object, "data.frame")
+  df_L <- data.frame(object@L)
+  colnames(df_L) <- paste0("L.", 1:J)
+  if (ncol(df_unmarkedFrame) > J) {
+    df <- cbind(df_unmarkedFrame[, 1:J], 
+                df_L, 
+                df_unmarkedFrame[, (J + 1):ncol(df_unmarkedFrame), drop = FALSE])
+  } else {
+    df <- cbind(df_unmarkedFrame[, 1:J], 
+                df_L)
+  }
+  cat("Data frame representation of unmarkedFrame object.\n")
+  print(df)
+})
+# LP note: as is defined in unmarkedFrame.R part "COERCION"
+
+
+## summary method ----
+setMethod("summary", "unmarkedFrameCOP", function(object,...) {
+  cat("unmarkedFrameCOP Object\n\n")
+  
+  cat(nrow(object@y), "sites\n")
+  cat("Maximum number of sampling occasions per site:",obsNum(object),"\n")
+  mean.obs <- mean(rowSums(!is.na(getY(object))))
+  cat("Mean number of sampling occasions per site:",round(mean.obs,2),"\n")
+  cat("Sites with at least one detection:",
+      sum(apply(getY(object), 1, function(x) any(x > 0, na.rm=TRUE))),
+      "\n\n")
+  
+  cat("Tabulation of y observations:")
+  print(table(object@y, exclude=NULL))
+  
+  if(!is.null(object@L)) {
+    cat("\nTabulation of sampling occasions length:")
+    print(table(object@L))
+  }
+  
+  if(!is.null(object@siteCovs)) {
+    cat("\nSite-level covariates:\n")
+    print(summary(object@siteCovs))
+  }
+  
+  if(!is.null(object@obsCovs)) {
+    cat("\nObservation-level covariates:\n")
+    print(summary(object@obsCovs))
+  }
+  
+})
+
+
+
+## umf[i, j] ----
+setMethod("[", c("unmarkedFrameCOP", "numeric", "numeric", "missing"),
+          function(x, i, j) {
+            # Gey dimensions of x
+            M <- numSites(x)
+            J <- obsNum(x)
+            
+            if (length(i) == 0 & length(j) == 0) {
+              return(x)
+            }
+            
+            # Check i
+            if (any(i < 0) &&
+                any(i > 0)) {
+              stop("i must be all positive or all negative indices.")
+            }
+            if (all(i < 0)) {
+              i <- (1:M)[i]
+            }
+            
+            # Check j
+            if (any(j < 0) &&
+                any(j > 0)) {
+              stop("j must be all positive or all negative indices.")
+            }
+            if (all(j < 0)) {
+              j <- (1:J)[j]
+            }
+            
+            # y observation count data subset
+            y <- getY(x)[i, j]
+            if (min(length(i), length(j)) == 1) {
+              y <- t(y)
+            }
+            
+            # L subset
+            L <- x@L[i, j]
+            if (min(length(i), length(j)) == 1) {
+              L <- t(L)
+            }
+            
+            # siteCovs subset
+            siteCovs <- siteCovs(x)
+            if (!is.null(siteCovs)) {
+              siteCovs <- siteCovs(x)[i, , drop = FALSE]
+            }
+            
+            # obsCovs subset
+            obsCovs <- obsCovs(x)
+            if (!is.null(obsCovs)) {
+              MJ_site <- rep(1:M, each = J)
+              MJ_obs <- rep(1:J, times = M)
+              obsCovs <- obsCovs[((MJ_obs %in% j) & (MJ_site %in% i)), , drop = FALSE]
+              rownames(obsCovs) <- NULL
+            }
+            
+            # Recreate umf
+            new(
+              Class = "unmarkedFrameCOP",
+              y = y,
+              L = L,
+              siteCovs = siteCovs,
+              obsCovs = obsCovs,
+              obsToY = diag(length(j)),
+              mapInfo = x@mapInfo
+            )
+          })
+
+
+## umf[i, ] ----
+setMethod("[", c("unmarkedFrameCOP", "numeric", "missing", "missing"),
+          function(x, i) {
+            x[i, 1:obsNum(x)]
+          })
+
+## umf[, j] ----
+setMethod("[", c("unmarkedFrameCOP", "missing", "numeric", "missing"),
+          function(x, j) {
+            x[1:numSites(x), j]
+          })
+
+
+## fl_getY ----
+setMethod("fl_getY", "unmarkedFitCOP", function(fit, ...){
+  getDesign(getData(fit), fit@formlist)$y
+})
+
+
+## predict_inputs_from_umf ----
+setMethod("predict_inputs_from_umf", "unmarkedFitCOP",
+          function(object, type, newdata, na.rm, re.form = NULL) {
+            designMats = getDesign(umf = newdata,
+                                   formlist = object@formlist,
+                                   na.rm = na.rm)
+            if (type == "psi") list_els <- c("Xpsi", "Zpsi")
+            if (type == "lambda") list_els <- c("Xlambda", "Zlambda")
+            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", "unmarkedFitCOP", function(object, type, ...) {
+  fl <- object@formlist
+  switch(type, psi = fl$psiformula, lambda = fl$lambdaformula)
+})
+
+
+## get_orig_data ----
+setMethod("get_orig_data", "unmarkedFitCOP", function(object, type, ...){
+  clean_covs <- clean_up_covs(object@data, drop_final=FALSE)
+  datatype <- switch(type, psi = 'site_covs', lambda = 'obs_covs')
+  clean_covs[[datatype]]
+})
+
+
+## getP ----
+setMethod("getP", "unmarkedFitCOP", function(object) {
+  data <- object@data
+  M = nrow(getY(data))
+  J = ncol(getY(data))
+  des <- getDesign(data, object@formlist, na.rm = na.rm)
+  matLambda =  do.call(object["lambda"]@invlink, 
+                       list(matrix(
+                         as.numeric(des$Xlambda %*% coef(object, 'lambda')),
+                         nrow = M, ncol = J, byrow = T)))
+  return(matLambda)
+})
+
+
+## fitted ----
+setMethod("fitted", "unmarkedFitCOP", function(object) {
+  data <- object@data
+  M = nrow(getY(data))
+  J = ncol(getY(data))
+  des <- getDesign(data, object@formlist, na.rm = na.rm)
+  estim_psi = as.numeric(do.call(object["psi"]@invlink,
+                                 list(as.matrix(des$Xpsi %*% coef(object, 'psi')))))
+  estim_lambda = do.call(object["lambda"]@invlink, 
+                         list(matrix(
+                           as.numeric(des$Xlambda %*% coef(object, 'lambda')),
+                           nrow = M, ncol = J, byrow = T)))
+  return(estim_psi * estim_lambda)
+})
+
+
+## residuals ----
+setMethod("residuals", "unmarkedFitCOP", function(object) {
+  y <- getY(object@data)
+  e <- fitted(object)
+  r <- y - e
+  return(r)
+})
+
+
+## plot ----
+setMethod("plot", c(x = "unmarkedFitCOP", y = "missing"), function(x, y, ...) {
+  y <- getY(x)
+  r <- residuals(x)
+  e <- fitted(x)
+  
+  old_graph <- graphics::par("mfrow", "mar")
+  on.exit(graphics::par(mfrow = old_graph$mfrow, mar = old_graph$mar))
+  
+  {
+    graphics::par(mfrow = c(2, 1))
+    graphics::par(mar = c(0, 5, 2, 2))
+    plot(e, y,
+         ylab = "Observed data",
+         xlab = "Predicted data",
+         xaxt = 'n')
+    abline(a = 0, b = 1, lty = 3, col = "red")
+    title(main = "COP model - detection events count", outer = F)
+    
+    graphics::par(mar = c(4, 5, 0.5, 2))
+    plot(e, r,
+         ylab = "Residuals",
+         xlab = "Predicted data")
+    abline(h = 0, lty = 3, col = "red")
+  }
+})
+
+
+
+## get_umf_components ----
+setMethod("get_umf_components", "unmarkedFitCOP",
+  function(object, formulas, guide, design, ...){
+    sc <- generate_data(formulas$psi, guide, design$M)
+    oc <- generate_data(formulas$lambda, guide, design$J*design$M)
+    yblank <- matrix(0, design$M, design$J)
+    list(y=yblank, siteCovs=sc, obsCovs=oc)
+})
+
+
+## simulate_fit ----
+setMethod("simulate_fit", "unmarkedFitCOP",
+  function(object, formulas, guide, design, ...){
+    # Generate covariates and create a y matrix of zeros
+    parts <- get_umf_components(object, formulas, guide, design, ...)
+    umf <- unmarkedFrameCOP(y = parts$y, siteCovs = parts$siteCovs, obsCovs=parts$obsCovs)
+    fit <- suppressMessages(
+      occuCOP(
+        data = umf,
+        psiformula = formula(formulas$psi),
+        lambdaformula = formula(formulas$lambda),
+        se = FALSE,
+        control = list(maxit = 1)
+      )
+    )
+    return(fit)
+})
+
+
+## simulate ----
+setMethod("simulate", "unmarkedFitCOP",
+  function(object, nsim = 1, seed = NULL, na.rm = TRUE){
+  set.seed(seed)
+  formula <- object@formula
+  umf <- object@data
+  designMats <- getDesign(umf = umf, formlist = object@formlist, na.rm = na.rm)
+  y <- designMats$y
+  M <- nrow(y)
+  J <- ncol(y)
+  
+  # Occupancy probability psi depending on the site covariates
+  psiParms = coef(object, type = "psi", fixedOnly = FALSE)
+  psi <- as.numeric(plogis(as.matrix(designMats$Xpsi %*% psiParms)))
+  
+  # Detection rate lambda depending on the observation covariates
+  lambda = getP(object = object)
+  
+  # Simulations
+  simList <- vector("list", nsim)
+  for(i in 1:nsim) {
+    Z <- rbinom(M, 1, psi)
+    # Z[object@knownOcc] <- 1
+    y = matrix(rpois(n = M * J, lambda = as.numeric(t(lambda))),
+               nrow = M, ncol = J, byrow = T) * Z
+    simList[[i]] <- y
+  }
+  return(simList)
+})
+
+
+## nonparboot ----
+setMethod("nonparboot", "unmarkedFitCOP",
+  function(object, B = 0, keepOldSamples = TRUE, ...) {
+  stop("Not currently supported for unmarkedFitCOP", call.=FALSE)
+})
+
+
+## ranef ----
+setMethod("ranef", "unmarkedFitCOP", function(object, ...) {
+  # Sites removed (srm) and sites kept (sk)
+  srm <- object@sitesRemoved
+  if (length(srm) > 0) {
+    sk = 1:numSites(getData(object))[-srm]
+  } else{
+    sk = 1:numSites(getData(object))
+  }
+  
+  # unmarkedFrame informations
+  M <- length(sk)
+  J <- obsNum(getData(object))
+  y <- getY(getData(object))[sk,]
+  
+  # Estimated parameters
+  psi <- predict(object, type = "psi")[sk, 1]
+  lambda <- getP(object)[sk,]
+  
+  # Estimate posterior distributions
+  z = c(0, 1)
+  post <- array(0, c(M, 2, 1), dimnames = list(NULL, z))
+  for (i in 1:M) {
+    # psi density
+    f <- dbinom(x = z,
+                size = 1,
+                prob = psi[i])
+    
+    # lambda density
+    g <- c(1, 1)
+    for (j in 1:J) {
+      if (is.na(y[i, j]) | is.na(lambda[i, j])) {
+        next
+      }
+      g = g * dpois(x = y[i, j], lambda = lambda[i, j] * z)
+    }
+    
+    # psi*lambda density
+    fudge <- f * g
+    post[i, , 1] <- fudge / sum(fudge)
+  }
+  
+  new("unmarkedRanef", post = post)
+})
+
+
+# Useful functions -------------------------------------------------------------
+
+check.integer <- function(x, na.ignore = F) {
+  if (is.na(x) & na.ignore) {
+    return(T)
+  }
+  x %% 1 == 0
+}
+
+# unmarkedFrame ----------------------------------------------------------------
+
+unmarkedFrameCOP <- function(y, L, siteCovs = NULL, obsCovs = NULL, mapInfo = NULL) {
+  
+  # Verification that they are non-NA data in y
+  if (all(is.na(y))) {
+    stop("y only contains NA. y needs to contain non-NA integers.")
+  }
+  
+  # Verification that these are count data (and not detection/non-detection)
+  if (max(y, na.rm = T) == 1) {
+    warning("unmarkedFrameCOP is for count data. ",
+            "The data furnished appear to be detection/non-detection.")
+  }
+  
+  # Number of sampling occasions
+  J <- ncol(y)
+  
+  # If missing L: replace by matrix of 1
+  # and the lambda will be the detection rate per observation length
+  if (missing(L)) {
+    L <- y * 0 + 1
+    warning("L is missing, replacing it by a matrix of 1.")
+  } else if (is.null(L)) {
+    L <- y * 0 + 1
+    warning("L is missing, replacing it by a matrix of 1.")
+  }
+  
+  # Transformation observation covariates
+  obsCovs <- covsToDF(
+    covs = obsCovs,
+    name = "obsCovs",
+    obsNum = J,
+    numSites = nrow(y)
+  )
+  
+  # Create S4 object of class unmarkedFrameCOP
+  umf <- new(
+    Class = "unmarkedFrameCOP",
+    y = y,
+    L = L,
+    siteCovs = siteCovs,
+    obsCovs = obsCovs,
+    obsToY = diag(J),
+    mapInfo = mapInfo
+  )
+  
+  return(umf)
+}
+
+
+
+# occuCOP ----------------------------------------------------------------------
+
+occuCOP <- function(data,
+                    psiformula = ~1,
+                    lambdaformula = ~1,
+                    psistarts,
+                    lambdastarts,
+                    method = "BFGS",
+                    se = TRUE,
+                    engine = c("C", "R"),
+                    threads = 1L,
+                    na.rm = TRUE,
+                    get.NLL.params = NULL,
+                    L1 = FALSE,
+                    ...) {
+  #TODO: engines
+  #TODO: random effects
+  
+  # Neg loglikelihood COP ------------------------------------------------------
+  R_nll_occuCOP <- function(params) {
+    
+    # Reading and transforming params
+    # Taking into account the covariates
+    
+    # psi as a function of covariates
+    #   psi in params are transformed with a logit transformation (qlogis)
+    #       so they're back-transformed to a proba with inverse logit (plogis)
+    #   Xpsi is the matrix with occupancy covariates
+    #   params is the vector with all the parameters
+    #   psiIdx is the index of Occupancy Parameters in params
+    psi <- plogis(Xpsi %*% params[psiIdx])
+    
+    # lambda as a function of covariates
+    #   lambda in params are transformed with a log-transformation
+    #          so they're back-transformed to a rate with exp here
+    #   Xlambda is the matrix with detection covariates
+    #   params is the vector with all the parameters
+    #   lambdaIdx is the index of Occupancy Parameters in params
+    lambda <- exp(Xlambda %*% params[lambdaIdx])
+    
+    # Listing sites analysed = sites not removed (due to NAs)
+    if (length(sitesRemoved) > 0) {
+      siteAnalysed = (1:M)[-sitesRemoved]
+    } else {
+      siteAnalysed = (1:M)
+    }
+    
+    # Probability for each site (i)
+    iProb <- rep(NA, M)
+    
+    for (i in siteAnalysed) {
+      # iIdx is the index to access the vectorised vectors of all obs in site i
+      iIdxall <- ((i - 1) * J + 1):(i * J)
+      
+      # Removing NAs
+      iIdx = iIdxall[!removed_obsvec[iIdxall]]
+      
+      if (SitesWithDetec[i]) {
+        # If there is at least one detection in site i
+        iProb[i] = psi[i] * (
+          (sum(lambda[iIdx] * Lvec[iIdx])) ^ sum(yvec[iIdx]) / 
+            factorial(sum(yvec[iIdx])) *
+            exp(-sum(lambda[iIdx] * Lvec[iIdx]))
+        )
+        
+      } else {
+        # If there is zero detection in site i
+        iProb[i] = psi[i] * exp(-sum(lambda[iIdx] * Lvec[iIdx])) + (1 - psi[i]) 
+        
+      }
+    }
+    
+    # log-likelihood
+    ll = sum(log(iProb[siteAnalysed]))
+    return(-ll)
+  }
+  
+  # Check arguments ------------------------------------------------------------
+  if (!is(data, "unmarkedFrameCOP")) {
+    stop("Data is not an unmarkedFrameCOP object. See ?unmarkedFrameCOP if necessary.")
+  }
+  stopifnot(class(psiformula) == "formula")
+  stopifnot(class(lambdaformula) == "formula")
+  if(!missing(psistarts)){stopifnot(class(psistarts) %in%  c("numeric", "double", "integer"))}
+  if(!missing(lambdastarts)){stopifnot(class(lambdastarts) %in%  c("numeric", "double", "integer"))}
+  stopifnot(class(threads) %in%  c("numeric", "double", "integer"))
+  se = as.logical(match.arg(
+    arg = as.character(se),
+    choices = c("TRUE", "FALSE", "0", "1")
+  ))
+  na.rm = as.logical(match.arg(
+    arg = as.character(na.rm),
+    choices = c("TRUE", "FALSE", "0", "1")
+  ))
+  engine <- match.arg(engine, c("C", "R"))
+  L1 = as.logical(match.arg(
+    arg = as.character(L1),
+    choices = c("TRUE", "FALSE", "0", "1")
+  ))
+  
+  
+  # Do not yet manage random effects!!!
+  if (has_random(psiformula) | has_random(lambdaformula)) {
+    stop("occuCOP does not currently handle random effects.")
+  }
+  
+  # Format input data ----------------------------------------------------------
+  
+  # Retrieve formlist
+  formlist <- mget(c("psiformula", "lambdaformula"))
+  
+  # Get the design matrix (calling the getDesign method for unmarkedFrame)
+  # For more informations, see: getMethod("getDesign", "unmarkedFrameCOP")
+  designMats <- getDesign(umf = data, formlist = formlist, na.rm = na.rm)
+  
+  # y is the count detection data (matrix of size M sites x J observations)
+  y <- getY(data)
+  
+  # L is the length of observations (matrix of size M sites x J observations)
+  L <- getL(data)
+  if (!L1) {
+    if (!any(is.na(L))) {
+      if (all(L == 1)) {
+        warning(
+          "All observations lengths (L) are set to 1. ",
+          "If they were not user-defined, lambda corresponds to the ",
+          "detection rate multiplied by the observation length, ",
+          "not just the detection rate per time-unit or space-unit.\n",
+          "You can remove this warning by adding 'L1=TRUE' in the function inputs."
+        )
+      }
+    }
+  }
+  
+  # Xpsi is the fixed effects design matrix for occupancy
+  Xpsi <- designMats$Xpsi
+  
+  # Xlambda is the fixed effects design matrix for detection rate
+  Xlambda <- designMats$Xlambda
+  
+  # Zpsi is the random effects design matrix for occupancy
+  Zpsi <- designMats$Zpsi
+  
+  # Zlambda is the random effects design matrix for detection rate
+  Zlambda <- designMats$Zlambda
+  
+  # removed_obs is a M x J matrix of the observations removed from the analysis
+  removed_obs <- designMats$removed_obs
+  sitesRemoved <- unname(which(apply(removed_obs, 1, function(x) all(x))))
+  
+  # Number of sites
+  M <- nrow(y)
+  
+  # Number of sampling occasions
+  J <- ncol(y)
+  
+  # Total number of detection per site
+  NbDetecPerSite = rowSums(y, na.rm=T)
+  
+  # Sites where there was at least one detection
+  SitesWithDetec = NbDetecPerSite > 0
+  
+  # Number of sites where there was at least one detection
+  NbSitesWithDetec = sum(SitesWithDetec)
+  
+  # Set up parameter names and indices-----------------------------------------
+  
+  # ParamPsi Occupancy parameter names
+  ParamPsi <- colnames(Xpsi)
+  
+  # ParamLambda Detection parameter names
+  ParamLambda <- colnames(Xlambda)
+  
+  # NbParamPsi Number of occupancy parameters
+  NbParamPsi <- ncol(Xpsi)
+  
+  # NbParamLambda Number of detection parameters
+  NbParamLambda <- ncol(Xlambda)
+  
+  # nP Total number of parameters
+  nP <- NbParamPsi + NbParamLambda
+  
+  # psiIdx Index of the occupancy parameters
+  psiIdx <- 1:NbParamPsi
+  
+  # lambdaIdx Index of the detection parameters
+  lambdaIdx <- (NbParamPsi+1):nP
+  
+  # Re-format some variables for C and R engines
+  #   From Matrix of dim MxJ to vector of length MxJ:
+  #   c(ySite1Obs1, ySite1Obs2, ..., ySite1ObsJ, ysite2Obs1, ...)
+  yvec <- as.numeric(t(y))
+  Lvec <- as.numeric(t(L))
+  removed_obsvec <- as.logical(t(removed_obs))
+  
+  # get.NLL.params -------------------------------------------------------------
+  if (!is.null(get.NLL.params)) {
+    df_NLL = data.frame(t(as.data.frame(get.NLL.params)))
+    rownames(df_NLL) = NULL
+    colnames(df_NLL) = c(paste0("logit(psi).", ParamPsi),
+                         paste0("log(lambda).", ParamLambda))
+    df_NLL$nll = NA
+    for (i in 1:nrow(df_NLL)) {
+      df_NLL$nll[i] = nll_COP(params = as.numeric(as.vector(df_NLL[i, -ncol(df_NLL)])))
+    }
+    return(df_NLL)
+  }
+  
+  # nll function depending on engine -------------------------------------------
+  if (identical(engine, "C")) {
+    nll <- function(params) {
+      nll_occuCOP(
+        y = yvec,
+        L = Lvec,
+        Xpsi = Xpsi,
+        Xlambda = Xlambda,
+        beta_psi = params[psiIdx],
+        beta_lambda = params[lambdaIdx],
+        removed = removed_obsvec
+      )
+    }
+  } else if (identical(engine, "R")) {
+    nll <- R_nll_occuCOP
+  }
+  
+  
+  # Optimisation ---------------------------------------------------------------
+  
+  ## Checking the starting point for optim
+  if (missing(lambdastarts)) {
+    # If lambda starts argument was not given:
+    # 0 by default 
+    # so lambda = exp(0) = 1 by default
+    lambdastarts = rep(0, NbParamLambda)
+    message(
+      "No lambda initial values provided for optim. Using lambdastarts = c(",
+      paste(lambdastarts, collapse = ", "),
+      "), equivalent to a detection rate",
+      ifelse(length(lambdastarts) == 1, "", "s"),
+      " of 1."
+    )
+  } else if (length(lambdastarts) != NbParamLambda) {
+    stop("lambdastarts (", paste(lambdastarts, collapse = ", "), ") ",
+         "should be of length ", NbParamLambda, " with lambdaformula ", lambdaformula)
+  }
+  
+  if (missing(psistarts)) {
+    # If psi starts argument was not given
+    # 0 by default 
+    # so psi = plogis(0) = 0.5 by default
+    psistarts = rep(0, NbParamPsi)
+    message(
+      "No psi initial values provided for optim. Using psistarts = c(",
+      paste(psistarts, collapse = ", "),
+      "), equivalent to an occupancy probabilit",
+      ifelse(length(psistarts) == 1, "y", "ies"),
+      " of 0.5."
+    )
+  } else if (length(psistarts) != NbParamPsi) {
+    stop("psistarts (", paste(psistarts, collapse = ", "), ") ",
+         "should be of length ", NbParamPsi, " with psiformula ", psiformula)
+  }
+  
+  starts <- c(psistarts, lambdastarts)
+  
+  ## Run optim
+  opt <- optim(
+    starts,
+    nll,
+    method = method,
+    hessian = se,
+    ...
+  )
+  
+  # Get output -----------------------------------------------------------------
+  covMat <- invertHessian(opt, nP, se)
+  ests <- opt$par
+  tmb_mod <- NULL
+  nll_value <- opt$value
+  fmAIC <- 2 * opt$value + 2 * nP
+  
+  # Organize effect estimates
+  names(ests) <- c(ParamPsi, ParamLambda)
+  psi_coef <- list(ests = ests[psiIdx], cov = as.matrix(covMat[psiIdx, psiIdx]))
+  lambda_coef <- list(ests = ests[lambdaIdx],
+                      cov = as.matrix(covMat[lambdaIdx, lambdaIdx]))
+  
+  # No random effects
+  psi_rand_info <- lambda_rand_info <- list()
+  
+  # Create unmarkedEstimates ---------------------------------------------------
+  psi_uE <- unmarkedEstimate(
+    name = "Occupancy probability",
+    short.name = "psi",
+    estimates = psi_coef$ests,
+    covMat = psi_coef$cov,
+    fixed = 1:NbParamPsi,
+    invlink = "logistic",
+    invlinkGrad = "logistic.grad",
+    randomVarInfo = psi_rand_info
+  )
+  
+  lambda_uE <- unmarkedEstimate(
+    name = "Detection rate",
+    short.name = "lambda",
+    estimates = lambda_coef$ests,
+    covMat = lambda_coef$cov,
+    fixed = 1:NbParamLambda,
+    invlink = "exp",
+    invlinkGrad = "exp",
+    randomVarInfo = lambda_rand_info
+  )
+  
+  estimateList <- unmarkedEstimateList(list(psi = psi_uE, lambda = lambda_uE))
+  
+  # Create unmarkedFit object--------------------------------------------------
+  umfit <- new(
+    "unmarkedFitCOP",
+    fitType = "occuCOP",
+    call = match.call(),
+    formula = as.formula(paste(
+      formlist["lambdaformula"], formlist["psiformula"], collapse = ""
+    )),
+    formlist = formlist,
+    data = data,
+    estimates = estimateList,
+    sitesRemoved = sitesRemoved,
+    removed_obs = removed_obs,
+    AIC = fmAIC,
+    opt = opt,
+    negLogLike = opt$value,
+    nllFun = nll,
+    nll = nll_value,
+    convergence = opt$convergence,
+    TMB = tmb_mod
+  )
+  
+  return(umfit)
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 6d89809..0fe45a8 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -338,6 +338,23 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// nll_occuCOP
+double nll_occuCOP(arma::icolvec y, arma::icolvec L, arma::mat Xpsi, arma::mat Xlambda, arma::colvec beta_psi, arma::colvec beta_lambda, Rcpp::LogicalVector removed);
+RcppExport SEXP _unmarked_nll_occuCOP(SEXP ySEXP, SEXP LSEXP, SEXP XpsiSEXP, SEXP XlambdaSEXP, SEXP beta_psiSEXP, SEXP beta_lambdaSEXP, SEXP removedSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::icolvec >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::icolvec >::type L(LSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type Xpsi(XpsiSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type Xlambda(XlambdaSEXP);
+    Rcpp::traits::input_parameter< arma::colvec >::type beta_psi(beta_psiSEXP);
+    Rcpp::traits::input_parameter< arma::colvec >::type beta_lambda(beta_lambdaSEXP);
+    Rcpp::traits::input_parameter< Rcpp::LogicalVector >::type removed(removedSEXP);
+    rcpp_result_gen = Rcpp::wrap(nll_occuCOP(y, L, Xpsi, Xlambda, beta_psi, beta_lambda, removed));
+    return rcpp_result_gen;
+END_RCPP
+}
 // nll_occuMS
 double nll_occuMS(arma::vec beta, arma::mat y, Rcpp::List dm_state, Rcpp::List dm_phi, Rcpp::List dm_det, arma::mat sind, arma::mat pind, arma::mat dind, std::string prm, int S, int T, int J, int N, arma::mat naflag, arma::mat guide);
 RcppExport SEXP _unmarked_nll_occuMS(SEXP betaSEXP, SEXP ySEXP, SEXP dm_stateSEXP, SEXP dm_phiSEXP, SEXP dm_detSEXP, SEXP sindSEXP, SEXP pindSEXP, SEXP dindSEXP, SEXP prmSEXP, SEXP SSEXP, SEXP TSEXP, SEXP JSEXP, SEXP NSEXP, SEXP naflagSEXP, SEXP guideSEXP) {
@@ -563,6 +580,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_unmarked_nll_multmixOpen", (DL_FUNC) &_unmarked_nll_multmixOpen, 39},
     {"_unmarked_nll_nmixTTD", (DL_FUNC) &_unmarked_nll_nmixTTD, 13},
     {"_unmarked_nll_occu", (DL_FUNC) &_unmarked_nll_occu, 11},
+    {"_unmarked_nll_occuCOP", (DL_FUNC) &_unmarked_nll_occuCOP, 7},
     {"_unmarked_nll_occuMS", (DL_FUNC) &_unmarked_nll_occuMS, 15},
     {"_unmarked_nll_occuMulti_loglik", (DL_FUNC) &_unmarked_nll_occuMulti_loglik, 14},
     {"_unmarked_nll_occuMulti", (DL_FUNC) &_unmarked_nll_occuMulti, 15},
diff --git a/src/nll_occuCOP.cpp b/src/nll_occuCOP.cpp
new file mode 100644
index 0000000..abae16c
--- /dev/null
+++ b/src/nll_occuCOP.cpp
@@ -0,0 +1,43 @@
+#include <RcppArmadillo.h>
+#include <float.h>
+
+using namespace Rcpp ;
+
+// [[Rcpp::export]]
+double nll_occuCOP(arma::icolvec y, arma::icolvec L, 
+    arma::mat Xpsi, arma::mat Xlambda,
+    arma::colvec beta_psi, arma::colvec beta_lambda,
+    Rcpp::LogicalVector removed) {
+  
+  // Number of sites M and obs J
+  int M = Xpsi.n_rows;
+  int J = y.n_elem / M;
+
+  //Calculate psi back-transformed from logit
+  arma::colvec psi = 1.0/(1.0+exp(-Xpsi*beta_psi));
+
+  //Calculate lambda back-transformed from log
+  arma::colvec lambda = exp(Xlambda*beta_lambda);
+  
+
+  double ll=0.0;
+  int k=0; // counter
+  // for each site i in 1:M
+  for(int i=0; i<M; i++) {
+    double iLambdaL=0.0; // init sum(lambda_ij * L_ij)
+    double iN=0.0; // init sum(y) = total count of detec at site i
+    for(int j=0; j<J; j++) {
+      if(!removed(k)) {
+        iLambdaL += lambda(k)*L(k);
+        iN += y(k);
+      }
+      k++;
+    }
+    if(iN>0) {
+      ll += log(psi(i) * pow(iLambdaL, iN) / tgamma(iN + 1) * exp(-iLambdaL));
+    } else {
+      ll += log(psi(i) * exp(-iLambdaL) + 1-psi(i));
+    }
+  }
+  return -ll;
+}
diff --git a/tests/testthat/test_occuCOP.R b/tests/testthat/test_occuCOP.R
new file mode 100644
index 0000000..c8b6738
--- /dev/null
+++ b/tests/testthat/test_occuCOP.R
@@ -0,0 +1,456 @@
+context("occuCOP fitting function")
+skip_on_cran()
+
+COPsimul <- function(psi = 0.5,
+                         lambda = 1,
+                         M = 100,
+                         J = 5) {
+  
+  z_i <- sample(
+    x = c(0, 1),
+    size = M,
+    prob = c(1 - psi, psi),
+    replace = T
+  )
+  
+  y = matrix(rpois(n = M * J, lambda = lambda), nrow = M, ncol = J) * z_i
+  
+  return(y)
+}
+
+
+test_that("unmarkedFrameCOP is constructed correctly", {
+  set.seed(123)
+  
+  # Simulate data
+  M = 100
+  J = 5
+  y = COPsimul(psi = .5,
+                   lambda = 1,
+                   M = M,
+                   J = J)
+  L = y * 0 + 1
+  
+  psiCovs = data.frame(
+    "psiNum" = rnorm(M),
+    "psiInt" = as.integer(rpois(n = M, lambda = 5)),
+    "psiBool" = sample(c(T, F), size = M, replace = T),
+    "psiChar" = sample(letters[1:5], size = M, replace = T),
+    "psiFactUnord" = factor(sample(
+      letters[5:10], size = M, replace = T
+    )),
+    "psiFactOrd" = sample(factor(c("A", "B", "C"), ordered = T), size =
+                            M, replace = T)
+  )
+  
+  lambdaCovs = list(
+    "lambdaNum" = matrix(
+      rnorm(M * J), 
+      nrow = M, ncol = J
+    ),
+    "lambdaInt" = matrix(
+      as.integer(rpois(n = M * J, lambda = 1e5)),
+      nrow = M, ncol = J
+    ),
+    "lambdaBool" = matrix(
+      sample(c(T, F), size = M * J, replace = T),
+      nrow = M, ncol = J
+    ),
+    "lambdaChar" = matrix(
+      sample(letters[1:5], size = M * J, replace = T),
+      nrow = M, ncol = J
+    ),
+    "lambdaFactUnord" = matrix(
+      factor(sample(letters[5:10], size = M * J, replace = T)),
+      nrow = M, ncol = J
+    ),
+    "lambdaFactOrd" = matrix(
+      sample(factor(c("A", "B", "C"), ordered = T), size = M * J, replace = T),
+      nrow = M, ncol = J
+    )
+  )
+
+  
+  # Creating a unmarkedFrameCOP object
+  expect_warning(umf <- unmarkedFrameCOP(y = y))
+  expect_no_error(umf <- unmarkedFrameCOP(y = y, L = L))
+  
+  
+  # Create subsets
+  expect_no_error(umf_sub_i <- umf[1:3, ])
+  expect_no_error(umf_sub_j <- umf[, 1:2])
+  expect_no_error(umf_sub_ij <- umf[1:3, 1:2])
+  
+  # unmarkedFrameCOP organisation ----------------------------------------------
+  expect_true(inherits(umf, "unmarkedFrameCOP"))
+  expect_equivalent(numSites(umf_sub_i), 3)
+  expect_equivalent(obsNum(umf_sub_j), 2)
+  expect_equivalent(numSites(umf_sub_ij), 3)
+  expect_equivalent(obsNum(umf_sub_ij), 2)
+  
+  # unmarkedFrameCOP display ---------------------------------------------------
+  
+  # print method
+  expect_output(print(umf), "Data frame representation of unmarkedFrame object")
+  expect_output(print(umf_sub_i), "Data frame representation of unmarkedFrame object")
+  expect_output(print(umf[1,]), "Data frame representation of unmarkedFrame object")
+  expect_output(print(umf[,1]), "Data frame representation of unmarkedFrame object")
+  expect_output(print(umf[1,1]), "Data frame representation of unmarkedFrame object")
+  
+  # summary method for unmarkedFrameCOP
+  expect_output(summary(umf), "unmarkedFrameCOP Object")
+  expect_output(summary(umf_sub_ij), "unmarkedFrameCOP Object")
+
+  # plot method for unmarkedFrameCOP
+  expect_no_error(plot(umf))
+  expect_no_error(plot(umf_sub_ij))
+  
+  
+  # Input handling: covariates -------------------------------------------------
+  expect_no_error(umfCovs <- unmarkedFrameCOP(
+    y = y,
+    L = L,
+    siteCovs = psiCovs,
+    obsCovs = lambdaCovs
+  ))
+  expect_output(print(umfCovs), "Data frame representation of unmarkedFrame object")
+  expect_output(summary(umfCovs), "unmarkedFrameCOP Object")
+  expect_no_error(plot(umfCovs))
+  
+  # Input handling: NA ---------------------------------------------------------
+  
+  # NA should not pose problems when creating the unmarkedFrameCOP object.
+  # The warnings and potential errors will be displayed when fitting the model.
+  # Except when y only contains NA: then there's an error.
+  
+  ## NA in y
+  yNA <- y
+  yNA[1:2,] <- NA
+  yNA[3:5, 3:4] <- NA
+  yNA[,3] <- NA
+  expect_no_error(umfNA <- unmarkedFrameCOP(y = yNA, L = L))
+  expect_output(print(umfNA), "Data frame representation of unmarkedFrame object")
+  expect_output(summary(umfNA), "unmarkedFrameCOP Object")
+  expect_no_error(plot(umfNA))
+  
+  ## NA in L
+  obsLengthNA <- L
+  obsLengthNA[1:2,] <- NA
+  obsLengthNA[3:5, 3:4] <- NA
+  obsLengthNA[,5] <- NA
+  expect_no_error(umfNA <- unmarkedFrameCOP(y = y, L = obsLengthNA))
+  expect_output(print(umfNA), "Data frame representation of unmarkedFrame object")
+  expect_output(summary(umfNA), "unmarkedFrameCOP Object")
+  expect_no_error(plot(umfNA))
+  
+  ## NA also in covariates
+  psiCovsNA <- psiCovs
+  psiCovsNA[1:5,] <- NA
+  psiCovsNA[c(8,10,12), 3] <- NA
+  psiCovsNA[,1] <- NA
+  lambdaCovsNA <- lambdaCovs
+  lambdaCovsNA[[1]][1:5,] <- NA
+  lambdaCovsNA[[1]][,3] <- NA
+  lambdaCovsNA[[3]][,5] <- NA
+  expect_no_error(umfCovsNA <- unmarkedFrameCOP(
+    y = yNA,
+    L = obsLengthNA,
+    siteCovs = psiCovsNA,
+    obsCovs = lambdaCovsNA
+  ))
+  expect_output(print(umfCovsNA), "Data frame representation of unmarkedFrame object")
+  expect_output(summary(umfCovsNA), "unmarkedFrameCOP Object")
+  expect_no_error(plot(umfCovsNA))
+  
+  ## all NA in y
+  yallNA <- y
+  yallNA[1:M, 1:J] <- NA
+  expect_error(unmarkedFrameCOP(y = yallNA, L = L))
+  
+  # Input handling: error and warnings -----------------------------------------
+  # Error when there are decimals in y
+  y_with_decimals = y
+  y_with_decimals[1, 1] = .5
+  expect_error(unmarkedFrameCOP(y = y_with_decimals, L = L))
+  
+  # Warning if y is detection/non-detection instead of count
+  y_detec_nodetec = (y > 0) * 1
+  expect_warning(unmarkedFrameCOP(y = y_detec_nodetec, L = L))
+  
+  # Error if the dimension of L is different than that of y
+  expect_error(unmarkedFrameCOP(y = y, L = L[1:2, 1:2]))
+})
+
+
+test_that("occuCOP can fit simple models", {
+  # Simulate data
+  set.seed(123)
+  M = 100
+  J = 5
+  y = COPsimul(psi = .5,
+                   lambda = 1,
+                   M = M,
+                   J = J)
+  L = y * 0 + 1
+
+  # Create umf
+  umf <- unmarkedFrameCOP(y = y, L = L)
+
+  # Fitting options ----
+
+  # With default parameters
+  expect_message(fit_default <- occuCOP(data = umf, L1=TRUE))
+  expect_warning(occuCOP(data = umf, psiformula = ~ 1, lambdaformula = ~ 1, psistarts = 0, lambdastarts = 0))
+
+  # With chosen starting points
+  expect_no_error(occuCOP(data = umf,
+                          psiformula = ~ 1, lambdaformula = ~ 1,
+                          psistarts = qlogis(.7),
+                          lambdastarts = log(0.1), L1=T))
+  expect_error(occuCOP(data = umf,
+                       psiformula = ~ 1, lambdaformula = ~ 1,
+                       psistarts = qlogis(c(0.7, 0.5), lambdastarts = 0, L1=T)))
+  expect_error(occuCOP(data = umf,
+                       psiformula = ~ 1, lambdaformula = ~ 1,
+                       lambdastarts = log(c(1, 2)), psistarts = 0, L1=T))
+
+  # With different options
+  expect_no_error(occuCOP(data = umf, method = "Nelder-Mead", psistarts = 0, lambdastarts = 0, L1=T))
+  expect_error(occuCOP(data = umf, method = "ABC", psistarts = 0, lambdastarts = 0, L1=T))
+
+  expect_no_error(occuCOP(data = umf, se = F, psistarts = 0, lambdastarts = 0, L1=T))
+  expect_error(occuCOP(data = umf, se = "ABC"))
+
+  expect_no_error(occuCOP(data = umf, engine = "R", psistarts = 0, lambdastarts = 0, L1=T))
+  expect_error(occuCOP(data = umf, engine = "julia", psistarts = 0, lambdastarts = 0, L1=T))
+
+  expect_no_error(occuCOP(data = umf, na.rm = F, psistarts = 0, lambdastarts = 0, L1=T))
+  expect_error(occuCOP(data = umf, na.rm = "no", psistarts = 0, lambdastarts = 0, L1=T))
+
+  # Looking at at COP model outputs ----
+  expect_is(fit_default, "unmarkedFitCOP")
+
+  ## backTransform
+  expect_no_error(backTransform(fit_default, type = "psi"))
+  expect_no_error(backTransform(fit_default, type = "lambda"))
+  expect_error(backTransform(fit_default, type = "state"))
+  expect_error(backTransform(fit_default, type = "det"))
+  expect_is(backTransform(fit_default, type = "psi"), "unmarkedBackTrans")
+
+  ## predict with newdata = fit@data
+  expect_no_error(predict(object = fit_default, type = "psi"))
+  expect_no_error(predict(object = fit_default, type = "lambda"))
+  expect_error(predict(object = fit_default, type = "state"))
+
+  ## predict with newdata = 1
+  expect_no_error(predict(
+    object = fit_default,
+    newdata = data.frame(1),
+    type = "psi"
+  ))
+  expect_no_error(predict(
+    object = fit_default,
+    newdata = data.frame(1),
+    type = "lambda"
+  ))
+  expect_no_error(predict(
+    object = fit_default,
+    newdata = data.frame("a"=1:5,"b"=10:14),
+    type = "psi"
+  ))
+
+  # Fitting accurately ----
+
+  ## psi = 0.50, lambda = 1 ----
+  psi_test = .5
+  lambda_test = 1
+  fit_accur <- occuCOP(data = unmarkedFrameCOP(
+    y = COPsimul(
+      psi = psi_test,
+      lambda = lambda_test,
+      M = 1000,
+      J = 10
+    ),
+    L = matrix(1, nrow = 1000, ncol = 10)
+  ), psistarts = 0, lambdastarts = 0, L1=T)
+  psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+  lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+  expect_equivalent(
+    psi_estimate,
+    psi_test,
+    tol = 0.05
+  )
+  expect_equivalent(
+    lambda_estimate,
+    lambda_test,
+    tol = 0.05
+  )
+
+  ## psi = 0.25, lambda = 5 ----
+  psi_test = 0.25
+  lambda_test = 5
+  fit_accur <- occuCOP(data = unmarkedFrameCOP(
+    y = COPsimul(
+      psi = psi_test,
+      lambda = lambda_test,
+      M = 1000,
+      J = 10
+    ),
+    L = matrix(1, nrow = 1000, ncol = 10)
+  ), psistarts = 0, lambdastarts = 0, L1=T)
+  psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+  lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+  expect_equivalent(
+    psi_estimate,
+    psi_test,
+    tol = 0.05
+  )
+  expect_equivalent(
+    lambda_estimate,
+    lambda_test,
+    tol = 0.05
+  )
+
+  ## psi = 0.75, lambda = 0.5 ----
+  psi_test = 0.75
+  lambda_test = 0.5
+  fit_accur <- occuCOP(data = unmarkedFrameCOP(
+    y = COPsimul(
+      psi = psi_test,
+      lambda = lambda_test,
+      M = 1000,
+      J = 10
+    ),
+    L = matrix(1, nrow = 1000, ncol = 10)
+  ), psistarts = 0, lambdastarts = 0, L1=T)
+  psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+  lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+  expect_equivalent(
+    psi_estimate,
+    psi_test,
+    tol = 0.05
+  )
+  expect_equivalent(
+    lambda_estimate,
+    lambda_test,
+    tol = 0.05
+  )
+
+  # With NAs ----
+  yNA <- y
+  yNA[1,] <- NA
+  yNA[3, 1] <- NA
+  yNA[4, 3] <- NA
+  yNA[, 5] <- NA
+  expect_no_error(umfNA <- unmarkedFrameCOP(y = yNA, L = L))
+
+  expect_warning(fit_NA <- occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0, L1=T))
+  expect_error(occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0, na.rm = F))
+})
+
+test_that("We can simulate COP data", {
+  
+  # From scratch ----
+  
+  # With no covariates
+  expect_no_error(simulate(
+    "COP",
+    formulas = list(psi =  ~ 1, lambda =  ~ 1),
+    coefs = list(
+      psi = c(intercept = 0),
+      lambda = c(intercept = 0)
+    ),
+    design = list(M = 100, J = 100)
+  ))
+  
+  # With quantitative covariates
+  expect_no_error(simulate(
+    "COP",
+    formulas = list(psi =  ~ elev, lambda =  ~ rain),
+    coefs = list(
+      psi = c(intercept = qlogis(.5), elev = -0.5),
+      lambda = c(intercept = log(3), rain = -1)
+    ),
+    design = list(M = 100, J = 5)
+  ))
+  
+  # With guides
+  expect_no_error(simulate(
+    "COP",
+    formulas = list(psi =  ~ elev, lambda =  ~ rain),
+    coefs = list(
+      psi = c(intercept = qlogis(.5), elev = -0.5),
+      lambda = c(intercept = log(3), rain = -1)
+    ),
+    design = list(M = 100, J = 5),
+    guide = list(elev=list(dist=rnorm, mean=12, sd=0.5))
+  ))
+  
+  # With qualitative covariates
+  expect_no_error(umf <- simulate(
+    "COP",
+    formulas = list(psi =  ~ elev + habitat, lambda =  ~ 1),
+    coefs = list(
+      psi = c(
+        intercept = qlogis(.2),
+        elev = -0.5,
+        habitatB = .5,
+        habitatC = .8
+      ),
+      lambda = c(intercept = log(3))
+    ),
+    design = list(M = 100, J = 5),
+    guide = list(habitat = factor(levels = c("A", "B", "C")))
+  ))
+  
+  # From unmarkedFitCOP ----
+  expect_no_error(umfit <- occuCOP(
+    umf,
+    psiformula =  ~ habitat,
+    L1 = T,
+    psistarts = c(0,0,0),
+    lambdastarts = 0
+  ))
+  expect_no_error(simulate(umfit))
+})
+
+test_that("occuCOP can fit models with covariates", {
+  # Simulate data with covariates ----
+  expect_no_error(umf <- simulate(
+    "COP",
+    formulas = list(psi =  ~ elev + habitat, lambda =  ~ rain),
+    coefs = list(
+      psi = c(
+        intercept = qlogis(.2),
+        elev = -0.5,
+        habitatB = .5,
+        habitatC = .8
+      ),
+      lambda = c(intercept = log(3), rain = -1)
+    ),
+    design = list(M = 100, J = 5),
+    guide = list(habitat = factor(levels = c("A", "B", "C")))
+  ))
+  
+  # Fit
+  expect_no_error(umfit <- occuCOP(
+    umf,
+    psiformula =  ~ habitat + elev,
+    lambdaformula =  ~ rain,
+    L1 = T,
+    psistarts = c(0,0,0,0),
+    lambdastarts = c(0,0)
+  ))
+  
+  expect_error(occuCOP(
+    umf,
+    psiformula =  ~ habitat+elev,
+    lambdaformula =  ~ rain,
+    L1 = T,
+    psistarts = c(0),
+    lambdastarts = c(0,0)
+  ))
+})
+
+
+