diff options
author | lpautrel <lea.pautrel@terroiko.fr> | 2023-12-12 12:26:04 +0100 |
---|---|---|
committer | lpautrel <lea.pautrel@terroiko.fr> | 2023-12-12 12:26:04 +0100 |
commit | 6484766345ebd17c60c214917f0c86c4af569765 (patch) | |
tree | 83ed2af0d19b09314cecc2ae65a5e35032b6bc02 | |
parent | 438a55109809e2bb5000d9f2e487aed18b2d06db (diff) |
Modifications PR #267 review
-rw-r--r-- | DESCRIPTION | 2 | ||||
-rw-r--r-- | R/occuCOP.R | 109 | ||||
-rw-r--r-- | R/simulate.R | 1 | ||||
-rw-r--r-- | man/getP-methods.Rd | 2 | ||||
-rw-r--r-- | man/occuCOP.Rd | 16 | ||||
-rw-r--r-- | man/unmarkedFrameCOP.Rd | 4 | ||||
-rw-r--r-- | tests/testthat/test_occuCOP.R | 48 |
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) }) |