diff options
author | Léa Pautrel <54315772+leapautrel@users.noreply.github.com> | 2023-12-07 10:03:43 +0100 |
---|---|---|
committer | GitHub <noreply@github.com> | 2023-12-07 10:03:43 +0100 |
commit | 0fa0a8986f0683df79035e0bafee54bfc3a532f3 (patch) | |
tree | 6f79f6277c6dfcaee002ef5a9203ce9477dd5d57 | |
parent | 591892700c69f343935fee52f6d37510d88bc9c3 (diff) | |
parent | 8111d55a5a42bb05c9ab0caebf54aed75efc1612 (diff) |
Merge pull request #1 from leapautrel/occuCOP2
Add occuCOP - debug package build
-rw-r--r-- | .gitignore | 4 | ||||
-rw-r--r-- | NAMESPACE | 3 | ||||
-rw-r--r-- | R/RcppExports.R | 4 | ||||
-rw-r--r-- | R/occuCOP.R | 170 | ||||
-rw-r--r-- | src/RcppExports.cpp | 18 | ||||
-rw-r--r-- | src/nll_occuCOP.cpp | 13 | ||||
-rw-r--r-- | tests/testthat/test_occuCOP.R | 240 |
7 files changed, 230 insertions, 222 deletions
@@ -27,5 +27,7 @@ runit bc .Rhistory symbols.rds -*.Rproj + +# Rproj .Rproj.user +*.Rproj
\ No newline at end of file @@ -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 |