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