Merge pull request #1 from leapautrel/occuCOP2

Add occuCOP - debug package build

7 files changed, 230 insertions, 222 deletions

diff --git a/.gitignore b/.gitignore
index 0a9ba4c..54223fc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -27,5 +27,7 @@ runit
 bc
 .Rhistory
 symbols.rds
-*.Rproj
+
+# Rproj
 .Rproj.user
+*.Rproj
\ No newline at end of file
diff --git a/NAMESPACE b/NAMESPACE
index 3666be0..b5aaf9f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -63,7 +63,8 @@ export(unmarkedEstimate, fitList, mapInfo, unmarkedFrame,
     unmarkedFrameDS, unmarkedMultFrame, unmarkedFrameGMM,
     unmarkedFramePCO, unmarkedFrameGDS, unmarkedFrameGPC, unmarkedFrameOccuMulti,
     unmarkedFrameOccuMS, unmarkedFrameOccuTTD, unmarkedFrameDSO, 
-    unmarkedFrameMMO, unmarkedFrameGDR, unmarkedFrameGOccu, unmarkedPostSamples)
+    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/occuCOP.R b/R/occuCOP.R
index 40601f3..02d351f 100644
--- a/R/occuCOP.R
+++ b/R/occuCOP.R
@@ -1,43 +1,25 @@
-# 
-# # Prerequisites (while its not integrated within the package)
-# setwd("~/synchros_git/MOBICO_git/unmarked/R")
-# library(unmarked)
-# source("unmarkedFrame.R")
-# source("utils.R")
-# source("unmarkedEstimate.R")
-# source("unmarkedFitList.R")
-# source("predict.R")
-# source("getDesign.R")
-# source("mixedModelTools.R")
-# source("simulate.R")
-# source("power.R")
-# source("boot.R")
-# # source("~/synchros_git/MOBICO_git/unmarked/R/RcppExports.R")
-# # sapply(paste0("./", list.files("./")), source)
-# 
-# library(lattice)
-# library(testthat)
-
 # Fit the occupancy model COP
 # (Counting Occurrences Process)
 
 # Occupancy
-# z_i ~ Bernoulli(psi)
+# z_i ~ Bernoulli(psi_i)
 # 
 # with:
 #   z_i = Occupancy state of site i
 #       = 1 if the site i is occupied
 #       = 0 else
+# with:
+#   psi_i = Occupancy probability of site i
 
 # Detection
-# N_ijs | z_i = 1 ~ Poisson(lambda*T_s)
-# N_ijs | z_i = 0 ~ 0
+# N_ij | z_i = 1 ~ Poisson(lambda_is*L_is)
+# N_ij | z_i = 0 ~ 0
 # 
 # with:
-#   N_is = Number of detections of site i during session s
+#   N_ij = Number of detections of site i during observation j
 #   z_i = Occupancy state of site i
-#   lambda = Detection rate
-#   T_s = Duration of the session s
+#   lambda_ij = Detection rate of the observation j in site i
+#   L_ij = Length/Duration of the observation j in site i
 
 # CLASSES ----------------------------------------------------------------------
 
@@ -60,9 +42,9 @@ setClass(
     }
     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)."
-    ))}
+        "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
   }
@@ -77,11 +59,9 @@ setClass("unmarkedFitCOP",
          contains = "unmarkedFit")
 
 
-# UNMARKED FRAME METHODS -------------------------------------------------------
+# Methods ----------------------------------------------------------------------
 
 ## getDesign method ----
-# Example for occu: https://github.com/rbchan/unmarked/blob/536e32ad7b2526f8fac1b315ed2e99accac0d50f/R/getDesign.R#L10-L97
-# Example for GDR: https://github.com/rbchan/unmarked/blob/c82e63947d7df7dfc896066e51dbf63bda3babf4/R/gdistremoval.R#L177-L253
 setMethod(
   "getDesign", "unmarkedFrameCOP",
   function(umf, formlist, na.rm = TRUE) {
@@ -101,7 +81,6 @@ setMethod(
     L <- getL(umf)
     
     # Occupancy submodel -------------------------------------------------------
-
     # Retrieve the fixed-effects part of the formula
     psiformula <- lme4::nobars(as.formula(formlist$psiformula))
     psiVars <- all.vars(psiformula)
@@ -187,12 +166,12 @@ setMethod(
         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). ",
+                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",
@@ -203,12 +182,12 @@ setMethod(
     
     # Output -------------------------------------------------------------------
     return(list(
-        y = y,
-        Xpsi = Xpsi,
-        Zpsi = Zpsi,
-        Xlambda = Xlambda,
-        Zlambda = Zlambda,
-        removed_obs = removed_obs
+      y = y,
+      Xpsi = Xpsi,
+      Zpsi = Zpsi,
+      Xlambda = Xlambda,
+      Zlambda = Zlambda,
+      removed_obs = removed_obs
     ))
   })
 
@@ -269,9 +248,9 @@ setMethod("summary", "unmarkedFrameCOP", function(object,...) {
     cat("\nObservation-level covariates:\n")
     print(summary(object@obsCovs))
   }
-
 })
 
+
 ## umf[i, j] ----
 setMethod("[", c("unmarkedFrameCOP", "numeric", "numeric", "missing"),
           function(x, i, j) {
@@ -359,10 +338,6 @@ setMethod("fl_getY", "unmarkedFitCOP", function(fit, ...){
   getDesign(getData(fit), fit@formlist)$y
 })
 
-# PREDICT METHODS ----
-# Fit type-specific methods to generate different components of prediction
-# 5. predict_by_chunk(): Take inputs and generate predictions
-# Basic methods are shown below; fit type-specific methods in their own sections
 
 ## predict_inputs_from_umf ----
 setMethod("predict_inputs_from_umf", "unmarkedFitCOP",
@@ -377,12 +352,14 @@ setMethod("predict_inputs_from_umf", "unmarkedFitCOP",
             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)
@@ -390,8 +367,6 @@ setMethod("get_orig_data", "unmarkedFitCOP", function(object, type, ...){
   clean_covs[[datatype]]
 })
 
-# UNMARKED FIT METHOD ----------------------------------------------------------
-# Required methods include: ranef, ,
 
 ## getP ----
 setMethod("getP", "unmarkedFitCOP", function(object) {
@@ -406,6 +381,7 @@ setMethod("getP", "unmarkedFitCOP", function(object) {
   return(matLambda)
 })
 
+
 ## fitted ----
 setMethod("fitted", "unmarkedFitCOP", function(object) {
   data <- object@data
@@ -421,6 +397,7 @@ setMethod("fitted", "unmarkedFitCOP", function(object) {
   return(estim_psi * estim_lambda)
 })
 
+
 ## residuals ----
 setMethod("residuals", "unmarkedFitCOP", function(object) {
   y <- getY(object@data)
@@ -429,9 +406,9 @@ setMethod("residuals", "unmarkedFitCOP", function(object) {
   return(r)
 })
 
+
 ## plot ----
-setMethod("plot", c(x = "unmarkedFitCOP", y = "missing"), function(x, y, ...)
-{
+setMethod("plot", c(x = "unmarkedFitCOP", y = "missing"), function(x, y, ...) {
   y <- getY(x)
   r <- residuals(x)
   e <- fitted(x)
@@ -455,10 +432,8 @@ setMethod("plot", c(x = "unmarkedFitCOP", y = "missing"), function(x, y, ...)
          xlab = "Predicted data")
     abline(h = 0, lty = 3, col = "red")
   }
-  
 })
 
-# SIMULATION METHODS -----------------------------------------------------------
 
 ## get_umf_components ----
 setMethod("get_umf_components", "unmarkedFitCOP",
@@ -469,6 +444,7 @@ setMethod("get_umf_components", "unmarkedFitCOP",
     list(y=yblank, siteCovs=sc, obsCovs=oc)
 })
 
+
 ## simulate_fit ----
 setMethod("simulate_fit", "unmarkedFitCOP",
   function(object, formulas, guide, design, ...){
@@ -487,46 +463,44 @@ setMethod("simulate_fit", "unmarkedFitCOP",
     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)
+  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)
-          })
-
-
+  function(object, B = 0, keepOldSamples = TRUE, ...) {
+  stop("Not currently supported for unmarkedFitCOP", call.=FALSE)
+})
 
-# POSTERIOR DISTRIBUTION METHODS -----------------------------------------------
 
 ## ranef ----
 setMethod("ranef", "unmarkedFitCOP", function(object, ...) {
@@ -574,7 +548,6 @@ setMethod("ranef", "unmarkedFitCOP", function(object, ...) {
 })
 
 
-
 # Useful functions -------------------------------------------------------------
 
 check.integer <- function(x, na.ignore = F) {
@@ -644,7 +617,7 @@ occuCOP <- function(data,
                     lambdastarts,
                     method = "BFGS",
                     se = TRUE,
-                    engine = "R",
+                    engine = c("C", "R"),
                     threads = 1L,
                     na.rm = TRUE,
                     get.NLL.params = NULL,
@@ -754,7 +727,7 @@ occuCOP <- function(data,
   y <- getY(data)
   
   # L is the length of observations (matrix of size M sites x J observations)
-  L = getL(data)
+  L <- getL(data)
   if (!L1) {
     if (!any(is.na(L))) {
       if (all(L == 1)) {
@@ -762,7 +735,7 @@ occuCOP <- function(data,
           "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.\n",
+          "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."
         )
       }
@@ -849,8 +822,8 @@ occuCOP <- function(data,
       nll_occuCOP(
         y = yvec,
         L = Lvec,
-        X = X,
-        V = V,
+        Xpsi = Xpsi,
+        Xlambda = Xlambda,
         beta_psi = params[psiIdx],
         beta_lambda = params[lambdaIdx],
         removed = removed_obsvec
@@ -897,13 +870,13 @@ occuCOP <- function(data,
     stop("psistarts (", paste(psistarts, collapse = ", "), ") ",
          "should be of length ", NbParamPsi, " with psiformula ", psiformula)
   }
-
+  
   starts <- c(psistarts, lambdastarts)
   
   ## Run optim
   opt <- optim(
     starts,
-    nll_occuCOP,
+    nll,
     method = method,
     hessian = se,
     ...
@@ -913,7 +886,7 @@ occuCOP <- function(data,
   covMat <- invertHessian(opt, nP, se)
   ests <- opt$par
   tmb_mod <- NULL
-  nll <- opt$value
+  nll_value <- opt$value
   fmAIC <- 2 * opt$value + 2 * nP
   
   # Organize effect estimates
@@ -966,12 +939,11 @@ occuCOP <- function(data,
     AIC = fmAIC,
     opt = opt,
     negLogLike = opt$value,
-    nllFun = nll_COP,
-    nll = nll,
+    nllFun = nll,
+    nll = nll_value,
     convergence = opt$convergence,
     TMB = tmb_mod
   )
   
   return(umfit)
-}
-
+}
\ No newline at end of file
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
index 6cb6cb7..abae16c 100644
--- a/src/nll_occuCOP.cpp
+++ b/src/nll_occuCOP.cpp
@@ -5,22 +5,19 @@ using namespace Rcpp ;
 
 // [[Rcpp::export]]
 double nll_occuCOP(arma::icolvec y, arma::icolvec L, 
-    arma::mat X, arma::mat V,
+    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 = X.n_rows;
+  int M = Xpsi.n_rows;
   int J = y.n_elem / M;
 
   //Calculate psi back-transformed from logit
-  arma::colvec logit_psi = X*beta_psi;
-  arma::colvec psi(M);
-  psi = 1.0/(1.0+exp(-logit_psi));
+  arma::colvec psi = 1.0/(1.0+exp(-Xpsi*beta_psi));
 
   //Calculate lambda back-transformed from log
-  arma::colvec log_lambda = V*beta_lambda;
-  arma::colvec lambda = exp(log_lambda);
+  arma::colvec lambda = exp(Xlambda*beta_lambda);
   
 
   double ll=0.0;
@@ -36,7 +33,7 @@ double nll_occuCOP(arma::icolvec y, arma::icolvec L,
       }
       k++;
     }
-    if(iN==0) {
+    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));
diff --git a/tests/testthat/test_occuCOP.R b/tests/testthat/test_occuCOP.R
index a4f0adb..ffecae0 100644
--- a/tests/testthat/test_occuCOP.R
+++ b/tests/testthat/test_occuCOP.R
@@ -19,23 +19,6 @@ COPsimul <- function(psi = 0.5,
 }
 
 
-# COPsimulfit <- function(psi = 0.5,
-#                         lambda = 1,
-#                         M = 100,
-#                         J = 5,
-#                         obsLength = NULL,
-#                         ...) {
-#   occuCOP(data = unmarkedFrameCOP(
-#     y = simulDataCOP(
-#       psi = psi,
-#       lambda = lambda,
-#       M = M,
-#       J = J
-#     ),
-#     obsLength = obsLength
-#   ), ...)
-# }
-
 test_that("unmarkedFrameCOP is constructed correctly", {
   set.seed(123)
   
@@ -46,7 +29,7 @@ test_that("unmarkedFrameCOP is constructed correctly", {
                    lambda = 1,
                    M = M,
                    J = J)
-  obsLength = y * 0 + 1
+  L = y * 0 + 1
   
   psiCovs = data.frame(
     "psiNum" = rnorm(M),
@@ -90,7 +73,7 @@ test_that("unmarkedFrameCOP is constructed correctly", {
   
   # Creating a unmarkedFrameCOP object
   expect_warning(umf <- unmarkedFrameCOP(y = y))
-  expect_no_error(umf <- unmarkedFrameCOP(y = y, obsLength = obsLength))
+  expect_no_error(umf <- unmarkedFrameCOP(y = y, L = L))
   
   
   # Create subsets
@@ -108,19 +91,15 @@ test_that("unmarkedFrameCOP is constructed correctly", {
   # unmarkedFrameCOP display ---------------------------------------------------
   
   # print method
-  expect_no_error(print(umf))
-  expect_no_warning(print(umf))
-  expect_no_message(print(umf))
-  expect_no_error(print(umf_sub_i))
-  expect_no_error(print(umf_sub_j))
-  expect_no_error(print(umf_sub_ij))
-  expect_no_error(print(umf[1,]))
-  expect_no_error(print(umf[,1]))
-  expect_no_error(print(umf[1,1]))
+  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_no_error(summary(umf))
-  expect_no_error(summary(umf_sub_ij))
+  expect_output(summary(umf), "unmarkedFrameCOP Object")
+  expect_output(summary(umf_sub_ij), "unmarkedFrameCOP Object")
   
   # plot method for unmarkedFrameCOP
   expect_no_error(plot(umf))
@@ -130,12 +109,12 @@ test_that("unmarkedFrameCOP is constructed correctly", {
   # Input handling: covariates -------------------------------------------------
   expect_no_error(umfCovs <- unmarkedFrameCOP(
     y = y,
-    obsLength = obsLength,
+    L = L,
     siteCovs = psiCovs,
     obsCovs = lambdaCovs
   ))
-  expect_no_error(print(umfCovs))
-  expect_no_error(summary(umfCovs))
+  expect_output(print(umfCovs), "Data frame representation of unmarkedFrame object")
+  expect_output(summary(umfCovs), "unmarkedFrameCOP Object")
   expect_no_error(plot(umfCovs))
   
   # Input handling: NA ---------------------------------------------------------
@@ -149,19 +128,20 @@ test_that("unmarkedFrameCOP is constructed correctly", {
   yNA[1:2,] <- NA
   yNA[3:5, 3:4] <- NA
   yNA[,3] <- NA
-  expect_no_error(umfNA <- unmarkedFrameCOP(y = yNA, obsLength = obsLength))
-  expect_no_error(print(umfNA))
-  expect_no_error(summary(umfNA))
+  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 obsLength
-  obsLengthNA <- obsLength
+  ## 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, obsLength = obsLengthNA))
-  expect_no_error(print(umfNA))
-  expect_no_error(summary(umfNA))
+  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
@@ -175,31 +155,31 @@ test_that("unmarkedFrameCOP is constructed correctly", {
   lambdaCovsNA[[3]][,5] <- NA
   expect_no_error(umfCovsNA <- unmarkedFrameCOP(
     y = yNA,
-    obsLength = obsLengthNA,
+    L = obsLengthNA,
     siteCovs = psiCovsNA,
     obsCovs = lambdaCovsNA
   ))
-  expect_no_error(print(umfCovsNA))
-  expect_no_error(summary(umfCovsNA))
+  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, obsLength = obsLength))
+  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, obsLength = obsLength))
+  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, obsLength = obsLength))
+  expect_warning(unmarkedFrameCOP(y = y_detec_nodetec, L = L))
   
-  # Error if the dimension of obsLength is different than that of y
-  expect_error(unmarkedFrameCOP(y = y, obsLength = obsLength[1:2, 1:2]))
+  # Error if the dimension of L is different than that of y
+  expect_error(unmarkedFrameCOP(y = y, L = L[1:2, 1:2]))
 })
 
 
@@ -212,42 +192,42 @@ test_that("occuCOP can fit simple models", {
                    lambda = 1,
                    M = M,
                    J = J)
-  obsLength = y * 0 + 1
-  
+  L = y * 0 + 1
+
   # Create umf
-  umf <- unmarkedFrameCOP(y = y, obsLength = obsLength)
-  
+  umf <- unmarkedFrameCOP(y = y, L = L)
+
   # Fitting options ----
-  
+
   # With default parameters
-  expect_message(fit_default <- occuCOP(data = umf))
-  expect_no_error(occuCOP(data = umf, psiformula = ~ 1, lambdaformula = ~ 1, psistarts = 0, lambdastarts = 0))
-  
+  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, 
+  expect_no_error(occuCOP(data = umf,
+                          psiformula = ~ 1, lambdaformula = ~ 1,
                           psistarts = qlogis(.7),
-                          lambdastarts = log(0.1)))
-  expect_error(occuCOP(data = umf, 
-                       psiformula = ~ 1, lambdaformula = ~ 1, 
-                       psistarts = qlogis(c(0.7, 0.5), lambdastarts = 0)))
-  expect_error(occuCOP(data = umf, 
-                       psiformula = ~ 1, lambdaformula = ~ 1, 
-                       lambdastarts = log(c(1, 2)), psistarts = 0))
-  
+                          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))
-  expect_error(occuCOP(data = umf, method = "ABC", psistarts = 0, lambdastarts = 0))
-  
-  expect_no_error(occuCOP(data = umf, se = F, psistarts = 0, lambdastarts = 0))
+  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))
-  expect_error(occuCOP(data = umf, engine = "julia", psistarts = 0, lambdastarts = 0))
-  
-  expect_no_error(occuCOP(data = umf, na.rm = F, psistarts = 0, lambdastarts = 0))
-  expect_error(occuCOP(data = umf, na.rm = "no", psistarts = 0, lambdastarts = 0))
-  
+
+  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")
   
@@ -281,7 +261,6 @@ test_that("occuCOP can fit simple models", {
   ))
   
   # Fitting accurately ----
-  
   ## psi = 0.50, lambda = 1 ----
   psi_test = .5
   lambda_test = 1
@@ -292,21 +271,21 @@ test_that("occuCOP can fit simple models", {
       M = 1000,
       J = 10
     ),
-    obsLength = matrix(1, nrow = 1000, ncol = 10)
-  ), psistarts = 0, lambdastarts = 0)
+    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, 
+    psi_estimate,
+    psi_test,
     tol = 0.05
   )
   expect_equivalent(
-    lambda_estimate, 
-    lambda_test, 
+    lambda_estimate,
+    lambda_test,
     tol = 0.05
   )
-  
+
   ## psi = 0.25, lambda = 5 ----
   psi_test = 0.25
   lambda_test = 5
@@ -317,21 +296,21 @@ test_that("occuCOP can fit simple models", {
       M = 1000,
       J = 10
     ),
-    obsLength = matrix(1, nrow = 1000, ncol = 10)
-  ), psistarts = 0, lambdastarts = 0)
+    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, 
+    psi_estimate,
+    psi_test,
     tol = 0.05
   )
   expect_equivalent(
-    lambda_estimate, 
-    lambda_test, 
+    lambda_estimate,
+    lambda_test,
     tol = 0.05
   )
-  
+
   ## psi = 0.75, lambda = 0.5 ----
   psi_test = 0.75
   lambda_test = 0.5
@@ -342,18 +321,18 @@ test_that("occuCOP can fit simple models", {
       M = 1000,
       J = 10
     ),
-    obsLength = matrix(1, nrow = 1000, ncol = 10)
-  ), psistarts = 0, lambdastarts = 0)
+    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, 
+    psi_estimate,
+    psi_test,
     tol = 0.05
   )
   expect_equivalent(
-    lambda_estimate, 
-    lambda_test, 
+    lambda_estimate,
+    lambda_test,
     tol = 0.05
   )
 
@@ -363,21 +342,12 @@ test_that("occuCOP can fit simple models", {
   yNA[3, 1] <- NA
   yNA[4, 3] <- NA
   yNA[, 5] <- NA
-  expect_no_error(umfNA <- unmarkedFrameCOP(y = yNA, obsLength = obsLength))
-  
-  expect_warning(fit_NA <- occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0))
-  expect_error(occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0, na.rm = F))
-  
-  
-  #TODO
-})
+  expect_no_error(umfNA <- unmarkedFrameCOP(y = yNA, L = L))
 
-
-test_that("occuCOP can fit models with covariates", {
-  
+  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 ----
@@ -417,7 +387,7 @@ test_that("We can simulate COP data", {
   ))
   
   # With qualitative covariates
-  expect_no_error(simulate(
+  expect_no_error(umf <- simulate(
     "COP",
     formulas = list(psi =  ~ elev + habitat, lambda =  ~ 1),
     coefs = list(
@@ -434,6 +404,50 @@ test_that("We can simulate COP data", {
   ))
   
   # From unmarkedFitCOP ----
-  #TODO
+  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)
+  ))
+})
\ No newline at end of file