Modifications PR #267 review

7 files changed, 114 insertions, 68 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 762d025..35cdf43 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -53,9 +53,9 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R'
     'simulate.R'
     'predict.R'
     'goccu.R'
+    'occuCOP.R'
     'RcppExports.R'
     'zzz.R'
-    'occuCOP.R'
 LinkingTo: 
     Rcpp,
     RcppArmadillo,
diff --git a/R/occuCOP.R b/R/occuCOP.R
index 04dba22..f5aff19 100644
--- a/R/occuCOP.R
+++ b/R/occuCOP.R
@@ -465,7 +465,8 @@ setMethod("simulate_fit", "unmarkedFitCOP",
 ## simulate ----
 setMethod("simulate", "unmarkedFitCOP",
   function(object, nsim = 1, seed = NULL, na.rm = TRUE){
-  set.seed(seed)
+  # set.seed(seed) 
+  # Purposefully not implemented
   formula <- object@formula
   umf <- object@data
   designMats <- getDesign(umf = umf, formlist = object@formlist, na.rm = na.rm)
@@ -613,11 +614,12 @@ occuCOP <- function(data,
                     lambdaformula = ~1,
                     psistarts,
                     lambdastarts,
+                    starts,
                     method = "BFGS",
                     se = TRUE,
                     engine = c("C", "R"),
                     na.rm = TRUE,
-                    get.NLL.params = NULL,
+                    return.negloglik = NULL,
                     L1 = FALSE,
                     ...) {
   #TODO: random effects
@@ -798,9 +800,9 @@ occuCOP <- function(data,
   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)))
+  # return.negloglik -----------------------------------------------------------
+  if (!is.null(return.negloglik)) {
+    df_NLL = data.frame(t(as.data.frame(return.negloglik)))
     rownames(df_NLL) = NULL
     colnames(df_NLL) = c(paste0("logit(psi).", ParamPsi),
                          paste0("log(lambda).", ParamLambda))
@@ -831,44 +833,67 @@ occuCOP <- function(data,
   
   # 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)
+  ## Checking the starting point for optim ----
+  # Check if either (psistarts AND lambdastarts) OR starts is provided
+  if (!missing(psistarts) & !missing(lambdastarts)) {
+    # Both psistarts and lambdastarts provided
+    if (!missing(starts)){
+      if (!all(c(psistarts, lambdastarts) == starts)) {
+        warning(
+          "You provided psistarts, lambdastarts and starts. ",
+          "Please provide either (psistarts AND lambdastarts) OR starts. ",
+          "Using psistarts and lambdastarts."
+        )
+      }
+    }
+    if (length(lambdastarts) != NbParamLambda) {
+      stop("lambdastarts (", paste(lambdastarts, collapse = ", "), ") ",
+           "should be of length ", NbParamLambda, " with lambdaformula ", lambdaformula)
+    }
+    if (length(psistarts) != NbParamPsi) {
+      stop("psistarts (", paste(psistarts, collapse = ", "), ") ",
+           "should be of length ", NbParamPsi, " with psiformula ", psiformula)
+    }
+    starts <- c(psistarts, lambdastarts)
+  } else if (!missing(starts)) {
+    # starts provided
+    if (length(starts) != nP) {
+      stop("starts (", paste(starts, collapse = ", "), ") ",
+           "should be of length ", nP, 
+           " with psiformula ", psiformula, 
+           " and lambdaformula ", lambdaformula)
+    }
+    
+    psistarts <- starts[1:NbParamPsi]
+    lambdastarts <- starts[(NbParamPsi + 1):(NbParamPsi + NbParamLambda)]
+    
+  } else {
+    # No arguments provided, apply default values
+    
+    if (missing(lambdastarts)) {
+      # If lambda starts argument was not given:
+      # 0 by default 
+      # so lambda = exp(0) = 1 by default
+      lambdastarts = rep(0, NbParamLambda)
+    } 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)
+    } else if (length(psistarts) != NbParamPsi) {
+      stop("psistarts (", paste(psistarts, collapse = ", "), ") ",
+           "should be of length ", NbParamPsi, " with psiformula ", psiformula)
+    }
+    
+    starts <- c(psistarts, lambdastarts)
   }
   
-  starts <- c(psistarts, lambdastarts)
-  
-  ## Run optim
+  ## Run optim ----
   opt <- optim(
     starts,
     nll,
@@ -939,4 +964,4 @@ occuCOP <- function(data,
   )
   
   return(umfit)
-}
\ No newline at end of file
+}
diff --git a/R/simulate.R b/R/simulate.R
index 9d93af9..a8887cb 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -61,7 +61,6 @@ blank_umFit <- function(fit_function){
 
 setMethod("simulate", "character",
   function(object, nsim=1, seed=NULL, formulas, coefs=NULL, design, guide=NULL, ...){
-  set.seed(seed)
   model <- blank_umFit(object)
   fit <- suppressWarnings(simulate_fit(model, formulas, guide, design, ...))
   coefs <- check_coefs(coefs, fit)
diff --git a/man/getP-methods.Rd b/man/getP-methods.Rd
index b7b8638..030c02d 100644
--- a/man/getP-methods.Rd
+++ b/man/getP-methods.Rd
@@ -20,7 +20,7 @@
 \alias{getP,unmarkedFitGOccu-method}
 \title{Methods for Function getP in Package `unmarked'}
 \description{
-Methods for function \code{getP} in Package `unmarked'. These methods return a matrix of dimension MxJ of the back-transformed detection parameter (\eqn{p} the detection probability or \eqn{\lambda} the detection rate, depending on the model), with M the number of sites and J the number of sampling periods.
+Methods for function \code{getP} in Package `unmarked'. These methods return a matrix of the back-transformed detection parameter (\eqn{p} the detection probability or \eqn{\lambda} the detection rate, depending on the model). The matrix is of dimension MxJ, with M the number of sites and J the number of sampling periods; or of dimension MxJT for models with multiple primary periods T.
 }
 \section{Methods}{
 \describe{
diff --git a/man/occuCOP.Rd b/man/occuCOP.Rd
index ca37e4e..aba49e3 100644
--- a/man/occuCOP.Rd
+++ b/man/occuCOP.Rd
@@ -7,10 +7,10 @@
 \usage{
 occuCOP(data, 
         psiformula = ~1, lambdaformula = ~1, 
-        psistarts, lambdastarts,
+        psistarts, lambdastarts, starts,
         method = "BFGS", se = TRUE, 
         engine = c("C", "R"), na.rm = TRUE, 
-        get.NLL.params = NULL, L1 = FALSE, ...)}
+        return.negloglik = NULL, L1 = FALSE, ...)}
 
 \arguments{
 
@@ -20,9 +20,11 @@ occuCOP(data,
     
     \item{lambdaformula}{Formula describing the detection covariates.}
     
-    \item{psistarts}{Vector of starting values for likelihood maximisation with \code{\link{optim}} for occupancy probability \eqn{\psi}{psi}. These values must be logit-transformed (with \code{\link{qlogis}}). See details. By default, optimisation will start at 0, corresponding to an occupancy probability of 0.5 (\code{plogis(0)} is 0.5).}
+    \item{psistarts}{Vector of starting values for likelihood maximisation with \code{\link{optim}} for occupancy probability \eqn{\psi}{psi}. These values must be logit-transformed (with \code{\link{qlogis}}) (see details). By default, optimisation will start at 0, corresponding to an occupancy probability of 0.5 (\code{plogis(0)} is 0.5).}
     
-    \item{lambdastarts}{Vector of starting values for likelihood maximisation with \code{\link{optim}} for detection rate \eqn{\lambda}{lambda}. These values must be log-transformed (with \code{\link{log}}). See details. By default, optimisation will start at 0, corresponding to detection rate of 1 (\code{exp(0)} is 1).}
+    \item{lambdastarts}{Vector of starting values for likelihood maximisation with \code{\link{optim}} for detection rate \eqn{\lambda}{lambda}. These values must be log-transformed (with \code{\link{log}}) (see details). By default, optimisation will start at 0, corresponding to detection rate of 1 (\code{exp(0)} is 1).}
+    
+    \item{starts}{Vector of starting values for likelihood maximisation with \code{\link{optim}}. If \code{psistarts} and \code{lambdastarts} are provided, \code{starts = c(psistarts, lambdastarts)}.}
 
     \item{method}{Optimisation method used by \code{\link{optim}}.}
     
@@ -32,7 +34,7 @@ occuCOP(data,
     
     \item{na.rm}{Logical specifying whether to fit the model (\code{na.rm=TRUE}) or not (\code{na.rm=FALSE}) if there are NAs in the \code{\link{unmarkedFrameCOP}} object.}
     
-    \item{get.NLL.params}{A list of vectors of parameters (\code{c(psiparams, lambdaparams)}). If specified, the function will not maximise likelihood but return the likelihood for the those parameters. See an example below.}
+    \item{return.negloglik}{A list of vectors of parameters (\code{c(psiparams, lambdaparams)}). If specified, the function will not maximise likelihood but return the negative log-likelihood for the those parameters. See an example below.}
 
     \item{L1}{Logical specifying whether the length of observations (\code{L}) are purposefully set to 1 (\code{L1=TRUE}) or not (\code{L1=FALSE}).}
 
@@ -222,11 +224,11 @@ occuCOP(data = umf, method = "Nelder-Mead")
 
 # We can run our model with a C++ or with a R likelihood function.
 ## They give the same result. 
-occuCOP(data = umf,engine = "C", psistarts = 0, lambdastarts = 0)
+occuCOP(data = umf, engine = "C", psistarts = 0, lambdastarts = 0)
 occuCOP(data = umf, engine = "R", psistarts = 0, lambdastarts = 0)
 
 ## The C++ (the default) is faster.
-system.time(occuCOP(data = umf,engine = "C", psistarts = 0, lambdastarts = 0))
+system.time(occuCOP(data = umf, engine = "C", psistarts = 0, lambdastarts = 0))
 system.time(occuCOP(data = umf, engine = "R", psistarts = 0, lambdastarts = 0))
 
 ## However, if you want to understand how the likelihood is calculated,
diff --git a/man/unmarkedFrameCOP.Rd b/man/unmarkedFrameCOP.Rd
index 937dcbe..a901b81 100644
--- a/man/unmarkedFrameCOP.Rd
+++ b/man/unmarkedFrameCOP.Rd
@@ -3,7 +3,7 @@
 
 \title{Organize data for the occupancy model using count data fit by \code{occuCOP}}
 
-\usage{unmarkedFrameCOP(y, L, siteCovs = NULL, obsCovs = NULL, mapInfo = NULL)}
+\usage{unmarkedFrameCOP(y, L, siteCovs = NULL, obsCovs = NULL)}
 
 \description{Organizes count data along with the covariates. The \linkS4class{unmarkedFrame} S4 class required by the \code{data} argument of \code{\link{occuCOP}}.}
 
@@ -17,7 +17,7 @@
   
   \item{obsCovs}{A named list of dataframes of dimension MxJ, with one dataframe per covariate that varies between sites and observation periods}
   
-  \item{mapInfo}{Currently ignored}
+  %% \item{mapInfo}{Currently ignored}
 }
 
 \details{
diff --git a/tests/testthat/test_occuCOP.R b/tests/testthat/test_occuCOP.R
index e698709..eca1130 100644
--- a/tests/testthat/test_occuCOP.R
+++ b/tests/testthat/test_occuCOP.R
@@ -199,21 +199,28 @@ test_that("occuCOP can fit simple models", {
 
   # Fitting options ----
 
-  # With default parameters
-  expect_message(fit_default <- occuCOP(data = umf, L1=TRUE))
+  ## With default parameters ----
+  expect_no_error(fit_default <- occuCOP(data = umf, L1 = TRUE))
   expect_warning(occuCOP(data = umf, psiformula = ~ 1, lambdaformula = ~ 1, psistarts = 0, lambdastarts = 0))
 
-  # With chosen starting points
+  ## 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))
+  expect_no_error(occuCOP(data = umf,
+                          psiformula = ~ 1, lambdaformula = ~ 1,
+                          starts = c(qlogis(.7), log(0.1)), L1 = T))
+  # warning if all starts and psistarts and lambdastarts were furnished
+  # and starts != c(psistarts, lambdastarts)
+  expect_no_error(occuCOP(data = umf, starts = c(0, 0),
+                          psistarts = c(0), lambdastarts = c(0), L1 = T))
+  expect_warning(occuCOP(data = umf, starts = c(0, 1),
+                          psistarts = c(0), lambdastarts = c(0), L1 = T))
+  # errors if starting vectors of the wrong length
+  expect_error(occuCOP(data = umf, starts = c(0), L1 = T))
+  expect_error(occuCOP(data = umf, psistarts = c(0, 0), lambdastarts = 0, L1 = T))
+  expect_error(occuCOP(data = umf, lambdastarts = c(0, 0), L1 = T))
 
   # With different options
   expect_no_error(occuCOP(data = umf, method = "Nelder-Mead", psistarts = 0, lambdastarts = 0, L1=T))
@@ -230,6 +237,7 @@ test_that("occuCOP can fit simple models", {
 
   # Looking at at COP model outputs ----
   expect_is(fit_default, "unmarkedFitCOP")
+  expect_equivalent(coef(fit_default), c(0.13067954, 0.06077929), tol = 1e-5)
   
   ## backTransform
   expect_no_error(backTransform(fit_default, type = "psi"))
@@ -239,8 +247,10 @@ test_that("occuCOP can fit simple models", {
   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_no_error(umpredpsi <- predict(object = fit_default, type = "psi"))
+  expect_equal(umpredpsi$Predicted[1], 0.5326235, tol = 1e-5)
+  expect_no_error(umpredlambda <- predict(object = fit_default, type = "lambda"))
+  expect_equal(umpredlambda$Predicted[1], 1.062664, tol = 1e-5)
   expect_error(predict(object = fit_default, type = "state"))
   
   ## predict with newdata = 1
@@ -416,6 +426,7 @@ test_that("We can simulate COP data", {
 
 test_that("occuCOP can fit and predict models with covariates", {
   # Simulate data with covariates ----
+  set.seed(123)
   expect_no_error(umf <- simulate(
     "COP",
     formulas = list(psi =  ~ elev + habitat, lambda =  ~ rain),
@@ -451,15 +462,24 @@ test_that("occuCOP can fit and predict models with covariates", {
     lambdastarts = c(0,0)
   ))
   
+  expect_equivalent(
+    coef(umfit), 
+    c(-1.5350679, 0.4229763, 0.7398768, -1.0456397, 1.2333424, -0.8344109), 
+    tol = 1e-5
+  )
+  
   # Predict ----
   expect_no_error(predict(umfit, type = "psi"))
-  expect_no_error(predict(umfit, type = "lambda", appendData = TRUE))
-  expect_no_error(predict(umfit, type = "lambda", level = 0.5))
-  expect_no_error(predict(
+  expect_no_error(umpredpsi <- predict(
     umfit,
     type = "psi",
     newdata = data.frame("habitat" = c("A", "B", "C"), "elev" = c(0, 0, 0)),
     appendData = TRUE
   ))
+  expect_equivalent(umpredpsi$Predicted, c(0.1772534, 0.2474811, 0.3110551), tol = 1e-5)
+  
+  expect_no_error(umpredlambda <- predict(umfit, type = "lambda", appendData = TRUE))
+  expect_no_error(predict(umfit, type = "lambda", level = 0.5))
+  expect_equal(umpredlambda$Predicted[1], 1.092008, tol = 1e-5)
 })