From 7b3c02a11a1383db77ae9480bc8e495f568d3c0d Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 30 Nov 2023 16:14:40 -0500 Subject: [PATCH 01/21] Initial commit of occuRNMulti function --- DESCRIPTION | 1 + NAMESPACE | 2 +- R/RcppExports.R | 4 + R/occuRNMulti.R | 173 +++++++++++++++++++++++++++++++ src/RcppExports.cpp | 22 ++++ src/nll_occuRNMulti.cpp | 218 ++++++++++++++++++++++++++++++++++++++++ 6 files changed, 419 insertions(+), 1 deletion(-) create mode 100644 R/occuRNMulti.R create mode 100644 src/nll_occuRNMulti.cpp diff --git a/DESCRIPTION b/DESCRIPTION index 87d8c4f6..260434d4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -46,6 +46,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R' 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R' 'nmixTTD.R' 'gdistremoval.R' + 'occuRNMulti.R' 'plotEffects.R' 'mixedModelTools.R' 'power.R' diff --git a/NAMESPACE b/NAMESPACE index 32349189..61da06b5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -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, occuRNMulti) export(removalPiFun, doublePiFun) export(makeRemPiFun, makeCrPiFun, makeCrPiFunMb, makeCrPiFunMh) diff --git a/R/RcppExports.R b/R/RcppExports.R index 79ab9671..8b206e5d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -69,6 +69,10 @@ nll_occuRN <- function(beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, thre .Call(`_unmarked_nll_occuRN`, beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, threads) } +nll_occuRNMulti <- function(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads) { + .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads) +} + nll_occuTTD <- function(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) { .Call(`_unmarked_nll_occuTTD`, beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) } diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R new file mode 100644 index 00000000..5653f287 --- /dev/null +++ b/R/occuRNMulti.R @@ -0,0 +1,173 @@ +setClass("unmarkedFrameOccuRNMulti", + contains="unmarkedFrameOccuMulti") + +setClass("unmarkedFitOccuRNMulti", + representation( + detformulas = "list", + stateformulas = "list"), + contains = "unmarkedFit") + +occuRNMulti <- function(detformulas, stateformulas, data, + K = rep(10, length(stateformulas)), + se = TRUE, threads=1){ + + stopifnot(all(names(data@ylist) == names(stateformulas))) + + S <- length(stateformulas) + dep <- create_dep_matrix(stateformulas) + + data <- as(data, "unmarkedFrameOccuRNMulti") + gd <- getDesign(data, detformulas, stateformulas) + + Kmin <- sapply(gd$ylist, function(x) apply(x, 1, max)) + + nP <- max(gd$det_ind) + state_rng <- min(gd$state_ind):max(gd$state_ind) + det_rng <- min(gd$det_ind):max(gd$det_ind) + starts <- rep(0, nP) + names(starts) <- gd$par_names + + nll_C <- function(pars){ + nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, + gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads) + } + + fm <- optim(starts, nll_C, method="BFGS", hessian=TRUE) + covMat <- invertHessian(fm, nP, se) + + fmAIC <- 2 * fm$value + 2 * nP + #---------------------------------------------------------------------------- + + #Format output--------------------------------------------------------------- + ests <- fm$par + names(ests) <- gd$par_names + + state <- unmarkedEstimate(name = "Abundance", short.name = "lam", + estimates = ests[state_rng], + covMat = as.matrix(covMat[state_rng, state_rng]), + invlink = "exp", + invlinkGrad = "exp") + + det <- unmarkedEstimate(name = "Detection", short.name = "p", + estimates = ests[det_rng], + covMat = as.matrix(covMat[det_rng, det_rng]), + invlink = "logistic", + invlinkGrad = "logistic.grad") + + estimateList <- unmarkedEstimateList(list(state=state, det=det)) + + umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(), + detformulas = detformulas, stateformulas = stateformulas, + formula = ~1, data = data, + #sitesRemoved = designMats$removed.sites, + estimates = estimateList, AIC = fmAIC, opt = fm, + negLogLike = fm$value, nllFun = nll_C) + + return(umfit) +} + +setMethod("getDesign", "unmarkedFrameOccuRNMulti", + function(umf, detformulas, stateformulas){ + M <- nrow(umf@ylist[[1]]) + J <- ncol(umf@ylist[[1]]) + + ylist <- umf@ylist + + # Site covs + sc <- umf@siteCovs + if(is.null(sc)) sc <- data.frame(dummy_=rep(1, M)) + + mm <- lapply(stateformulas, function(s){ + if(!is.list(s)){ + return(model.matrix(s, sc)) + } + lapply(s, function(x) model.matrix(x, sc)) + }) + + # Obs covs + sc_long_ind <- rep(1:M, each=J) + sc_long <- sc[sc_long_ind,] + oc <- cbind(sc_long, umf@obsCovs) + + vv <- lapply(detformulas, function(x){ + model.matrix(x, oc) + }) + + # Indices + state_n <- unlist(lapply(mm, function(x){ + if(is.list(x)){ + lapply(x, function(y) ncol(y)) + } else { + ncol(x) + } + })) + state_end <- cumsum(state_n) + state_start <- state_end - state_n + 1 + state_ind <- cbind(state_start, state_end) + + state_prefix <- rep(rownames(state_ind), state_n) + state_prefix <- gsub(".", ":", state_prefix, fixed=TRUE) + state_prefix <- paste0("[",state_prefix,"]") + + state_par <- unlist(lapply(mm, function(x){ + if(is.matrix(x)) return(colnames(x)) + lapply(x, function(y){ + return(colnames(y)) + }) + })) + stopifnot(length(state_prefix) == length(state_par)) + state_par <- paste(state_prefix, state_par) + + det_n <- sapply(vv, ncol) + det_end <- cumsum(det_n) + max(state_end) + det_start <- det_end - det_n + 1 + det_ind <- cbind(det_end, det_start) + + det_rng <- min(det_ind):max(det_ind) + + det_prefix <- paste0("[",rep(rownames(det_ind), det_n),"]") + det_par <- unlist(lapply(vv, colnames)) + stopifnot(length(det_prefix) == length(det_par)) + det_par <- paste(det_prefix, det_par) + par_names <- c(state_par, det_par) + + list(par_names = par_names, state_ind = state_ind, det_ind = det_ind, + ylist = ylist, state_dm = mm, det_dm = vv) +}) + +create_dep_matrix <- function(stateformulas){ + S <- length(stateformulas) + dep <- matrix(0, S, S) + rownames(dep) <- colnames(dep) <- names(stateformulas) + allowed_species <- character(0) + + # Species 1 + stopifnot(is.call(stateformulas[[1]])) + allowed_species <- c(allowed_species, names(stateformulas)[1]) + +# Species 2 + if(is.list(stateformulas[[2]])){ + dep_sp <- names(stateformulas[[2]])[2:length(stateformulas[[2]])] + stopifnot(all(dep_sp %in% allowed_species)) + for (i in dep_sp){ + dep[names(stateformulas)[2], dep_sp] <- 1 + } + } + allowed_species <- c(allowed_species, names(stateformulas)[2]) + + # Species 3 + if(S > 2){ + if(is.list(stateformulas[[3]])){ + dep_sp <- names(stateformulas[[3]])[2:length(stateformulas[[3]])] + stopifnot(all(dep_sp %in% allowed_species)) + for (i in dep_sp){ + dep[names(stateformulas)[3], dep_sp] <- 1 + } + } + allowed_species <- c(allowed_species, names(stateformulas)[3]) + } + + dep +} + + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 6d89809a..5d2eae36 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -453,6 +453,27 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// nll_occuRNMulti +double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads); +RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type state_ind(state_indSEXP); + Rcpp::traits::input_parameter< arma::mat >::type det_ind(det_indSEXP); + Rcpp::traits::input_parameter< int >::type S(SSEXP); + Rcpp::traits::input_parameter< Rcpp::List >::type ylist(ylistSEXP); + Rcpp::traits::input_parameter< Rcpp::List >::type Xlist(XlistSEXP); + Rcpp::traits::input_parameter< Rcpp::List >::type Vlist(VlistSEXP); + Rcpp::traits::input_parameter< arma::imat >::type dep(depSEXP); + Rcpp::traits::input_parameter< arma::ivec >::type K(KSEXP); + Rcpp::traits::input_parameter< arma::imat >::type Kmin(KminSEXP); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads)); + return rcpp_result_gen; +END_RCPP +} // nll_occuTTD double nll_occuTTD(arma::vec beta, arma::vec y, arma::vec delta, arma::mat W, arma::mat V, arma::mat Xgam, arma::mat Xeps, arma::vec pind, arma::vec dind, arma::vec cind, arma::vec eind, std::string lpsi, std::string tdist, int N, int T, int J, arma::vec naflag); RcppExport SEXP _unmarked_nll_occuTTD(SEXP betaSEXP, SEXP ySEXP, SEXP deltaSEXP, SEXP WSEXP, SEXP VSEXP, SEXP XgamSEXP, SEXP XepsSEXP, SEXP pindSEXP, SEXP dindSEXP, SEXP cindSEXP, SEXP eindSEXP, SEXP lpsiSEXP, SEXP tdistSEXP, SEXP NSEXP, SEXP TSEXP, SEXP JSEXP, SEXP naflagSEXP) { @@ -568,6 +589,7 @@ static const R_CallMethodDef CallEntries[] = { {"_unmarked_nll_occuMulti", (DL_FUNC) &_unmarked_nll_occuMulti, 15}, {"_unmarked_nll_occuPEN", (DL_FUNC) &_unmarked_nll_occuPEN, 11}, {"_unmarked_nll_occuRN", (DL_FUNC) &_unmarked_nll_occuRN, 10}, + {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 11}, {"_unmarked_nll_occuTTD", (DL_FUNC) &_unmarked_nll_occuTTD, 17}, {"_unmarked_nll_pcount", (DL_FUNC) &_unmarked_nll_pcount, 11}, {"_unmarked_nll_pcountOpen", (DL_FUNC) &_unmarked_nll_pcountOpen, 35}, diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp new file mode 100644 index 00000000..d4059262 --- /dev/null +++ b/src/nll_occuRNMulti.cpp @@ -0,0 +1,218 @@ +#include +#include "utils.h" + +#ifdef _OPENMP + #include +#endif + +using namespace Rcpp; +using namespace arma; + +//arma::mat inv_logit( arma::mat inp ){ +// return(1 / (1 + exp(-1 * inp))); +//} + +//Calculate part of abundance likelihood to integrate over (abun f[g] * detect [g]) +double calc_fg(int k, double lam, int J, arma::vec r, arma::rowvec y){ + double f = R::dpois(k, lam, 0); + double p; + double g = 0; + for (int j = 0; j(ylist[0]); + unsigned M = y1.n_rows; + unsigned J = y1.n_cols; + + //Parameters for dominant species + mat X1 = as(Xlist[0]); + vec lam1 = exp(X1 * beta.subvec(state_ind(0, 0), state_ind(0, 1))); + + mat V1 = as(Vlist[0]); + vec r1 = inv_logit(V1 * beta.subvec(det_ind(0, 0), det_ind(0,1))); + + //Parameters for 2nd species + mat y2 = as(ylist[1]); + + int sind_row = 1; + + vec lam2; + vec lp2, lp2_sp1; + if(dep(1,0)){ + List Xlist_sp2 = Xlist[1]; + mat X2 = as(Xlist_sp2[0]); + lp2 = X2 * beta.subvec(state_ind(1,0), state_ind(1,1)); + mat X2_sp1 = as(Xlist_sp2[1]); + lp2_sp1 = X2_sp1 * beta.subvec(state_ind(2,0), state_ind(2,1)); + sind_row += 2; + } else { + mat X2 = as(Xlist[1]); + lam2 = exp(X2 * beta.subvec(state_ind(1,0), state_ind(1,1))); + sind_row += 1; + } + + mat V2 = as(Vlist[1]); + vec r2 = inv_logit(V2 * beta.subvec(det_ind(1, 0), det_ind(1,1))); + + //Parameters for 3rd species (if needed) + mat y3; + vec r3; + vec lam3, lp3, lp3_sp1, lp3_sp2; + + //Very janky... + if(S > 2){ + y3 = as(ylist[2]); + + if(dep(2,0) | dep(2,1)){ + List Xlist_sp3 = Xlist[2]; + mat X3 = as(Xlist_sp3[0]); + lp3 = X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1)); + sind_row += 1; + + int mat_idx = 1; + if(dep(2,0)){ + mat X3_sp1 = as(Xlist_sp3[mat_idx]); + mat_idx += 1; + lp3_sp1 = X3_sp1 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1)); + sind_row += 1; + } + if(dep(2,1)){ + mat X3_sp2 = as(Xlist_sp3[mat_idx]); + lp3_sp2 = X3_sp2 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1)); + } + } else { + mat X3 = as(Xlist[2]); + lam3 = exp(X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1))); + } + + mat V3 = as(Vlist[2]); + r3 = inv_logit(V3 * beta.subvec(det_ind(2, 0), det_ind(2,1))); + } + + #pragma omp parallel for reduction(-: nll) if(threads > 1) + for (int m = 0; m sp2 + + lik_sp2 = 0; + + //ALLOW OCCUPANCY HERE + lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m)); + + lik_m += fg1 * lik_sp2; + + //----------------------------------------------------------------------- + } else if((S == 3) & (dep(1,0) & dep(2,1))){ // sp1 --> sp2 --> sp3 + + lik_sp2 = 0; + + lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); + + for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ + fg2 = calc_fg(k2, lam2_m, J, r2_sub, y2.row(m)); + + lik_sp3 = 0; + + //ALLOW OCCUPANCY HERE + lam3_m = exp(lp3(m) + k2 * lp3_sp2(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); + + lik_sp2 += fg2 * lik_sp3; + } + + lik_m += fg1 * lik_sp2; + + //----------------------------------------------------------------------- + } else if((S == 3) & (dep(2,0) & dep(2,1))){ // sp1 --> sp3 <-- sp2 + + lik_sp2 = 0; + + for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ + + fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(m)); + + lik_sp3 = 0; + + //ALLOW OCCUPANCY HERE + lam3_m = exp(lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); + + lik_sp2 += fg2 * lik_sp3; + } + + lik_m += fg1 * lik_sp2; + + //----------------------------------------------------------------------- + } else if((S == 3) & (dep(1,0) & dep(2,0))){ // sp2 <-- sp1 --> sp3 + + lik_sp2 = 0; + //ALLOW OCCUPANCY HERE + lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m)); + + lik_sp3 = 0; + + //ALLOW OCCUPANCY HERE + lam3_m = exp(lp3(m) + k1 * lp3_sp1(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); + + lik_m += fg1 * lik_sp2 * lik_sp3; + } + } + + + nll -= log(lik_m); + } + return nll; +} From 731137db3ae2bdd050c6fbbdc568dfda2dd2e4e2 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 1 Dec 2023 11:46:27 -0500 Subject: [PATCH 02/21] Add occupancy support, model sp1-->sp3<--sp2 doesnt work --- R/RcppExports.R | 4 +- R/occuRNMulti.R | 72 +++++++++++++++++++++++++++++------- src/RcppExports.cpp | 9 +++-- src/nll_occuRNMulti.cpp | 81 ++++++++++++++++++++++++----------------- 4 files changed, 114 insertions(+), 52 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 8b206e5d..3d86707d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -69,8 +69,8 @@ nll_occuRN <- function(beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, thre .Call(`_unmarked_nll_occuRN`, beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, threads) } -nll_occuRNMulti <- function(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads) { - .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads) +nll_occuRNMulti <- function(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads) { + .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads) } nll_occuTTD <- function(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) { diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 5653f287..a8278c1a 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -7,14 +7,31 @@ setClass("unmarkedFitOccuRNMulti", stateformulas = "list"), contains = "unmarkedFit") -occuRNMulti <- function(detformulas, stateformulas, data, +occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, K = rep(10, length(stateformulas)), - se = TRUE, threads=1){ + starts, method="BFGS", se = TRUE, threads=1, ...){ stopifnot(all(names(data@ylist) == names(stateformulas))) S <- length(stateformulas) dep <- create_dep_matrix(stateformulas) + + if(missing(modelOccupancy)){ + modelOccupancy <- rep(0, S) + names(modelOccupancy) <- names(stateformulas) + } else { + stopifnot(identical(names(modelOccupancy), names(stateformulas))) + if(is.logical(modelOccupancy)){ + modelOccupancy <- modelOccupancy*1 + } + not_allowed_occ <- colSums(dep) + sp_mismatch <- modelOccupancy & not_allowed_occ + if(any(sp_mismatch)){ + stop(paste("Species", + paste(names(modelOccupancy)[sp_mismatch], collapse=", "), + "cannot use occupancy models."), call.=FALSE) + } + } data <- as(data, "unmarkedFrameOccuRNMulti") gd <- getDesign(data, detformulas, stateformulas) @@ -22,17 +39,32 @@ occuRNMulti <- function(detformulas, stateformulas, data, Kmin <- sapply(gd$ylist, function(x) apply(x, 1, max)) nP <- max(gd$det_ind) - state_rng <- min(gd$state_ind):max(gd$state_ind) + + has_occ <- sum(modelOccupancy) > 0 + if(has_occ){ + sp_occ <- names(modelOccupancy)[modelOccupancy == 1] + sp_abun <- names(modelOccupancy)[modelOccupancy == 0] + occ_ind <- gd$state_ind[grepl(paste(paste0("^",sp_occ), collapse="|"), + rownames(gd$state_ind)),,drop=FALSE] + occ_rng <- min(occ_ind):max(occ_ind) + abun_ind <- gd$state_ind[grepl(paste(paste0("^",sp_abun), collapse="|"), + rownames(gd$state_ind)),,drop=FALSE] + abun_rng <- min(abun_ind):max(abun_ind) + } else { + state_rng <- min(gd$state_ind):max(gd$state_ind) + } + det_rng <- min(gd$det_ind):max(gd$det_ind) - starts <- rep(0, nP) + + if(missing(starts)) starts <- rep(0, nP) names(starts) <- gd$par_names nll_C <- function(pars){ - nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, + nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, modelOccupancy, gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads) } - fm <- optim(starts, nll_C, method="BFGS", hessian=TRUE) + fm <- optim(starts, nll_C, method=method, hessian=se, ...) covMat <- invertHessian(fm, nP, se) fmAIC <- 2 * fm$value + 2 * nP @@ -42,19 +74,33 @@ occuRNMulti <- function(detformulas, stateformulas, data, ests <- fm$par names(ests) <- gd$par_names - state <- unmarkedEstimate(name = "Abundance", short.name = "lam", - estimates = ests[state_rng], - covMat = as.matrix(covMat[state_rng, state_rng]), - invlink = "exp", - invlinkGrad = "exp") - det <- unmarkedEstimate(name = "Detection", short.name = "p", estimates = ests[det_rng], covMat = as.matrix(covMat[det_rng, det_rng]), invlink = "logistic", invlinkGrad = "logistic.grad") - estimateList <- unmarkedEstimateList(list(state=state, det=det)) + if(has_occ){ + lam <- unmarkedEstimate(name = "Abundance", short.name = "lam", + estimates = ests[abun_rng], + covMat = as.matrix(covMat[abun_rng, abun_rng]), + invlink = "exp", + invlinkGrad = "exp") + psi <- unmarkedEstimate(name = "Occupancy", short.name = "psi", + estimates = ests[occ_rng], + covMat = as.matrix(covMat[occ_rng, occ_rng]), + invlink = "logistic", + invlinkGrad = "logistic.grad") + + estimateList <- unmarkedEstimateList(list(lam=lam, psi=psi, det=det)) + } else { + lam <- unmarkedEstimate(name = "Abundance", short.name = "lam", + estimates = ests[state_rng], + covMat = as.matrix(covMat[state_rng, state_rng]), + invlink = "exp", + invlinkGrad = "exp") + estimateList <- unmarkedEstimateList(list(lam=lam, det=det)) + } umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(), detformulas = detformulas, stateformulas = stateformulas, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5d2eae36..ec850b27 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -454,8 +454,8 @@ BEGIN_RCPP END_RCPP } // nll_occuRNMulti -double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads); -RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, arma::ivec modOcc, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads); +RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP modOccSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -463,6 +463,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::mat >::type state_ind(state_indSEXP); Rcpp::traits::input_parameter< arma::mat >::type det_ind(det_indSEXP); Rcpp::traits::input_parameter< int >::type S(SSEXP); + Rcpp::traits::input_parameter< arma::ivec >::type modOcc(modOccSEXP); Rcpp::traits::input_parameter< Rcpp::List >::type ylist(ylistSEXP); Rcpp::traits::input_parameter< Rcpp::List >::type Xlist(XlistSEXP); Rcpp::traits::input_parameter< Rcpp::List >::type Vlist(VlistSEXP); @@ -470,7 +471,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::ivec >::type K(KSEXP); Rcpp::traits::input_parameter< arma::imat >::type Kmin(KminSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads)); + rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads)); return rcpp_result_gen; END_RCPP } @@ -589,7 +590,7 @@ static const R_CallMethodDef CallEntries[] = { {"_unmarked_nll_occuMulti", (DL_FUNC) &_unmarked_nll_occuMulti, 15}, {"_unmarked_nll_occuPEN", (DL_FUNC) &_unmarked_nll_occuPEN, 11}, {"_unmarked_nll_occuRN", (DL_FUNC) &_unmarked_nll_occuRN, 10}, - {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 11}, + {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 12}, {"_unmarked_nll_occuTTD", (DL_FUNC) &_unmarked_nll_occuTTD, 17}, {"_unmarked_nll_pcount", (DL_FUNC) &_unmarked_nll_pcount, 11}, {"_unmarked_nll_pcountOpen", (DL_FUNC) &_unmarked_nll_pcountOpen, 35}, diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp index d4059262..5402a9a3 100644 --- a/src/nll_occuRNMulti.cpp +++ b/src/nll_occuRNMulti.cpp @@ -8,10 +8,6 @@ using namespace Rcpp; using namespace arma; -//arma::mat inv_logit( arma::mat inp ){ -// return(1 / (1 + exp(-1 * inp))); -//} - //Calculate part of abundance likelihood to integrate over (abun f[g] * detect [g]) double calc_fg(int k, double lam, int J, arma::vec r, arma::rowvec y){ double f = R::dpois(k, lam, 0); @@ -29,7 +25,7 @@ double lik_subord_abun(int Kmin, int Kmax, double lam, int J, arma::vec r, arma::rowvec y){ double lik = 0.0; - + lam = exp(lam); for (int k = Kmin; k < (Kmax+1); k++){ lik += calc_fg(k, lam, J, r, y); } @@ -37,10 +33,20 @@ double lik_subord_abun(int Kmin, int Kmax, double lam, int J, arma::vec r, return lik; } +//Calculate occupancy likelihood for subordinate species +double lik_subord_occ(int Kmin, double psi, int J, arma::vec p, arma::rowvec y){ + double g = 1; + psi = inv_logit(psi); + for (int j = 0; j(Xlist[2]); lam3 = exp(X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1))); } @@ -134,7 +141,7 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i double lik_m = 0.0; - double lam2_m, lam3_m; + double par2_m, par3_m; // state parameters for species 2 and 3 for (int k1 = Kmin(m, 0); k1<(K(0)+1); k1++){ @@ -145,11 +152,12 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i if((S == 2) & dep(1,0)){ //sp1 --> sp2 lik_sp2 = 0; - - //ALLOW OCCUPANCY HERE - lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); - lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m)); - + par2_m = lp2(m) + k1 * lp2_sp1(m); + if(modOcc(1)){ + lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m)); + } else { + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m)); + } lik_m += fg1 * lik_sp2; //----------------------------------------------------------------------- @@ -157,17 +165,18 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i lik_sp2 = 0; - lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); + par2_m = exp(lp2(m) + k1 * lp2_sp1(m)); for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ - fg2 = calc_fg(k2, lam2_m, J, r2_sub, y2.row(m)); + fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m)); lik_sp3 = 0; - - //ALLOW OCCUPANCY HERE - lam3_m = exp(lp3(m) + k2 * lp3_sp2(m)); - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); - + par3_m = lp3(m) + k2 * lp3_sp2(m); + if(modOcc(2)){ + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + } else { + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + } lik_sp2 += fg2 * lik_sp3; } @@ -182,12 +191,13 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(m)); + par3_m = lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m); lik_sp3 = 0; - - //ALLOW OCCUPANCY HERE - lam3_m = exp(lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m)); - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); - + if(modOcc(1)){ + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + } else { + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + } lik_sp2 += fg2 * lik_sp3; } @@ -196,22 +206,27 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i //----------------------------------------------------------------------- } else if((S == 3) & (dep(1,0) & dep(2,0))){ // sp2 <-- sp1 --> sp3 + + par2_m = lp2(m) + k1 * lp2_sp1(m); lik_sp2 = 0; - //ALLOW OCCUPANCY HERE - lam2_m = exp(lp2(m) + k1 * lp2_sp1(m)); - lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m)); + if(modOcc(1)){ + lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m)); + } else { + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m)); + } + par3_m = lp3(m) + k1 * lp3_sp1(m); lik_sp3 = 0; - - //ALLOW OCCUPANCY HERE - lam3_m = exp(lp3(m) + k1 * lp3_sp1(m)); - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m)); + if(modOcc(2)){ + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + } else { + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + } lik_m += fg1 * lik_sp2 * lik_sp3; } } - nll -= log(lik_m); } return nll; From 3633e82add842fb9cd753147789890e05008a2ce Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 18 Dec 2023 15:46:49 -0500 Subject: [PATCH 03/21] Fix bug with sp1 --> sp3 <-- sp2 --- src/nll_occuRNMulti.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp index 5402a9a3..1bf54c9d 100644 --- a/src/nll_occuRNMulti.cpp +++ b/src/nll_occuRNMulti.cpp @@ -193,7 +193,7 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i par3_m = lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m); lik_sp3 = 0; - if(modOcc(1)){ + if(modOcc(2)){ lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); } else { lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); From 359725e221908abb90450d53ec7a01f8cd632aef Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 29 Mar 2024 21:34:15 -0400 Subject: [PATCH 04/21] Add sp1 --> sp2 --> sp3 & sp1 --> sp3 model cpp code --- src/nll_occuRNMulti.cpp | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp index 1bf54c9d..14a7434a 100644 --- a/src/nll_occuRNMulti.cpp +++ b/src/nll_occuRNMulti.cpp @@ -143,6 +143,8 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i double par2_m, par3_m; // state parameters for species 2 and 3 + //This can probably be hugely simplified by combining the cases + //in a smart way; kept them all separate for simplicity initially for (int k1 = Kmin(m, 0); k1<(K(0)+1); k1++){ //Dominant species @@ -224,6 +226,27 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i } lik_m += fg1 * lik_sp2 * lik_sp3; + } else if((S == 3) & (dep(1,0) & dep(2,0) & dep(2,1))){ // sp1 --> sp2 --> sp3 & sp1 --> sp3 + + lik_sp2 = 0; + + par2_m = exp(lp2(m) + k1 * lp2_sp1(m)); + + for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ + fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m)); + + lik_sp3 = 0; + par3_m = lp3(m) + k1 * lp2_sp1(m) + k2 * lp3_sp2(m); + if(modOcc(2)){ + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + } else { + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + } + lik_sp2 += fg2 * lik_sp3; + } + + lik_m += fg1 * lik_sp2; + } } From a6171dcb1c625d07b8fb2f4bc9127b2678033768 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 10 Jul 2024 14:39:23 -0400 Subject: [PATCH 05/21] Add predict method --- R/occuRNMulti.R | 140 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 139 insertions(+), 1 deletion(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index a8278c1a..ede71d7e 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -29,7 +29,7 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, if(any(sp_mismatch)){ stop(paste("Species", paste(names(modelOccupancy)[sp_mismatch], collapse=", "), - "cannot use occupancy models."), call.=FALSE) + "cannot use occupancy model."), call.=FALSE) } } @@ -216,4 +216,142 @@ create_dep_matrix <- function(stateformulas){ dep } +setMethod("predict", "unmarkedFitOccuRNMulti", + function(object, type, species, newdata, + level = 0.95, nsims = 100, ...){ + stopifnot(type %in% c("state", "det")) + if(!missing(species)){ + stopifnot(species %in% names(object@data@ylist)) + } + if(type == "state"){ + samps <- MASS::mvrnorm(nsims, mu=object@estimates@estimates$state@estimates, + Sigma = object@estimates@estimates$state@covMat) + + # Function to extract correct parts of state coefs vector + get_species_b <- function(b, object, sp){ + nm <- names(b) + allsp <- names(object@data@ylist) + ind <- which(allsp == sp) + out <- list() + other <- c(1:length(allsp))[-ind] + match1 <- grepl(paste0("[", allsp[ind], "]"), nm, fixed=TRUE) + out[[1]] <- b[match1] + + match2 <- grepl(paste0("[", sp, ":", allsp[other[1]],"]"), nm, fixed=TRUE) + if(sum(match2 > 0)){ + out[[allsp[other[1]]]] <- b[match2] + } + match3 <- grepl(paste0("[", sp, ":", allsp[other[2]],"]"), nm, fixed=TRUE) + if(sum(match3 > 0)){ + out[[allsp[other[2]]]] <- b[match3] + } + out + } + + depmat <- create_dep_matrix(object@stateformulas) + + top <- apply(depmat, 1, function(x) sum(x) == 0) + sp_top <- colnames(depmat)[top] + lev2 <- apply(depmat[!top,!top,drop=FALSE], 1, function(x) sum(x) == 0) + sp_2nd <- names(lev2)[lev2] + sp_3rd <- names(lev2)[!lev2] + + if(missing(newdata) || is.null(newdata)){ + + X <- lapply(object@stateformulas, function(x){ + if(is.list(x)){ + lapply(x, function(z) model.matrix(z, object@data@siteCovs)) + } else { + model.matrix(x, object@data@siteCovs) + } + }) + nr <- nrow(object@data@siteCovs) + } else { + X <- lapply(object@stateformulas, function(x){ + if(is.list(x)){ + lapply(x, function(z) unmarked:::make_mod_matrix(z, object@data@siteCovs, newdata)$X) + } else { + unmarked:::make_mod_matrix(x, object@data@siteCovs, newdata)$X + } + }) + nr <- nrow(newdata) + } + + coef_est <- matrix(object@estimates@estimates$state@estimates, nrow=1) + colnames(coef_est) <- names(object@estimates@estimates$state@estimates) + point_est <- calc_dependent_abundance(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd) + + message('Bootstrapping confidence intervals with ',nsims,' samples') + post <- calc_dependent_abundance(samps, X, object, nr, sp_top, sp_2nd, sp_3rd) + + # Summarize bootstrap + cis <- lapply(post, function(x){ + data.frame( + SE = apply(x, 1, sd, na.rm=TRUE), + upper = apply(x, 1, quantile, 0.025, na.rm=TRUE), + lower = apply(x, 1, quantile, 0.975, na.rm=TRUE) + ) + }) + + out <- mapply(function(x, y){ + cbind(Predicted = x, y) + }, point_est, cis, SIMPLIFY=FALSE) + + if(!missing(species)){ + out <- out[[species]] + } + return(out) + } else if(type == "det"){ + stop("det not supported yet") + } +}) + +calc_dependent_abundance <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd){ + nsims <- nrow(samps) + + all_species <- names(object@data@ylist) + + post <- lapply(1:length(all_species), function(x) matrix(NA, nr, nsims)) + names(post) <- all_species + + for (k in 1:nsims){ + + b <- samps[k,] + ests <- list() + + # Top-level species that don't depend on other species abundance + for (i in sp_top){ + beta <- get_species_b(b, object, i)[[1]] + ests[[i]] <- exp(X[[i]] %*% beta) + post[[i]][,k] <- ests[[i]] + } + # 2nd-level species that depend only on top-level species + for (i in sp_2nd){ + Xsub <- X[[i]] + beta <- get_species_b(b, object, i) + lp <- Xsub[[1]] %*% beta[[1]] + for (j in 2:length(Xsub)){ + other_sp <- names(Xsub)[j] + add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] + lp <- lp + add + } + ests[[i]] <- exp(lp) + post[[i]][,k] <- ests[[i]] + } + # 3rd-level species that depend on at least one 2nd-level species + for (i in sp_3rd){ + Xsub <- X[[i]] + beta <- get_species_b(b, object, i) + lp <- Xsub[[1]] %*% beta[[1]] + for (j in 2:length(Xsub)){ + other_sp <- names(Xsub)[j] + add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] + lp <- lp + add + } + ests[[i]] <- exp(lp) + post[[i]][,k] <- ests[[i]] + } + } + post +} From ff063b6be3582fb33d87e7875edba5a094e9f746 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 10 Jul 2024 14:57:28 -0400 Subject: [PATCH 06/21] Fix cpp problem --- src/nll_occuRNMulti.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp index 14a7434a..0de39b96 100644 --- a/src/nll_occuRNMulti.cpp +++ b/src/nll_occuRNMulti.cpp @@ -55,8 +55,8 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i double nll = 0; mat y1 = as(ylist[0]); - unsigned M = y1.n_rows; - unsigned J = y1.n_cols; + int M = y1.n_rows; + int J = y1.n_cols; //Parameters for dominant species mat X1 = as(Xlist[0]); From 378e84aa7b7c943e65e504008a49143127bcde8d Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 12 Jul 2024 16:25:47 -0400 Subject: [PATCH 07/21] Predict method works with occupancy species --- R/occuRNMulti.R | 78 ++++++++++++++++++++++++++++--------------------- 1 file changed, 45 insertions(+), 33 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index ede71d7e..4ff5f7e0 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -4,7 +4,8 @@ setClass("unmarkedFrameOccuRNMulti", setClass("unmarkedFitOccuRNMulti", representation( detformulas = "list", - stateformulas = "list"), + stateformulas = "list", + modelOccupancy = "numeric"), contains = "unmarkedFit") occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, @@ -104,6 +105,7 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(), detformulas = detformulas, stateformulas = stateformulas, + modelOccupancy = modelOccupancy, formula = ~1, data = data, #sitesRemoved = designMats$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, @@ -225,31 +227,9 @@ setMethod("predict", "unmarkedFitOccuRNMulti", stopifnot(species %in% names(object@data@ylist)) } if(type == "state"){ - samps <- MASS::mvrnorm(nsims, mu=object@estimates@estimates$state@estimates, - Sigma = object@estimates@estimates$state@covMat) - - # Function to extract correct parts of state coefs vector - get_species_b <- function(b, object, sp){ - nm <- names(b) - allsp <- names(object@data@ylist) - ind <- which(allsp == sp) - out <- list() - other <- c(1:length(allsp))[-ind] - match1 <- grepl(paste0("[", allsp[ind], "]"), nm, fixed=TRUE) - out[[1]] <- b[match1] - - match2 <- grepl(paste0("[", sp, ":", allsp[other[1]],"]"), nm, fixed=TRUE) - if(sum(match2 > 0)){ - out[[allsp[other[1]]]] <- b[match2] - } - match3 <- grepl(paste0("[", sp, ":", allsp[other[2]],"]"), nm, fixed=TRUE) - if(sum(match3 > 0)){ - out[[allsp[other[2]]]] <- b[match3] - } - out - } + samps <- MASS::mvrnorm(nsims, mu=coef(object), Sigma = vcov(object)) - depmat <- create_dep_matrix(object@stateformulas) + depmat <- unmarked:::create_dep_matrix(object@stateformulas) top <- apply(depmat, 1, function(x) sum(x) == 0) sp_top <- colnames(depmat)[top] @@ -277,13 +257,13 @@ setMethod("predict", "unmarkedFitOccuRNMulti", }) nr <- nrow(newdata) } - - coef_est <- matrix(object@estimates@estimates$state@estimates, nrow=1) - colnames(coef_est) <- names(object@estimates@estimates$state@estimates) - point_est <- calc_dependent_abundance(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd) + + coef_est <- matrix(coef(object), nrow=1) + colnames(coef_est) <- names(coef(object)) + point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink) message('Bootstrapping confidence intervals with ',nsims,' samples') - post <- calc_dependent_abundance(samps, X, object, nr, sp_top, sp_2nd, sp_3rd) + post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink) # Summarize bootstrap cis <- lapply(post, function(x){ @@ -307,7 +287,29 @@ setMethod("predict", "unmarkedFitOccuRNMulti", } }) -calc_dependent_abundance <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd){ +# Function to extract correct parts of state coefs vector +get_species_b <- function(b, object, sp){ + b <- b[!grepl("^p\\(", names(b))] # remove detection coefficients + nm <- names(b) + allsp <- names(object@data@ylist) + ind <- which(allsp == sp) + out <- list() + other <- c(1:length(allsp))[-ind] + match1 <- grepl(paste0("[", allsp[ind], "]"), nm, fixed=TRUE) + out[[1]] <- b[match1] + + match2 <- grepl(paste0("[", sp, ":", allsp[other[1]],"]"), nm, fixed=TRUE) + if(sum(match2 > 0)){ + out[[allsp[other[1]]]] <- b[match2] + } + match3 <- grepl(paste0("[", sp, ":", allsp[other[2]],"]"), nm, fixed=TRUE) + if(sum(match3 > 0)){ + out[[allsp[other[2]]]] <- b[match3] + } + out +} + +calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink){ nsims <- nrow(samps) all_species <- names(object@data@ylist) @@ -328,6 +330,11 @@ calc_dependent_abundance <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3r } # 2nd-level species that depend only on top-level species for (i in sp_2nd){ + if(object@modelOccupancy[i]){ + ilink <- plogis + } else { + ilink <- eval(str2lang(object@estimates@estimates$lam@invlink)) + } Xsub <- X[[i]] beta <- get_species_b(b, object, i) lp <- Xsub[[1]] %*% beta[[1]] @@ -336,11 +343,16 @@ calc_dependent_abundance <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3r add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] lp <- lp + add } - ests[[i]] <- exp(lp) + ests[[i]] <- ilink(lp) post[[i]][,k] <- ests[[i]] } # 3rd-level species that depend on at least one 2nd-level species for (i in sp_3rd){ + if(object@modelOccupancy[i]){ + ilink <- plogis + } else { + ilink <- eval(str2lang(object@estimates@estimates$lam@invlink)) + } Xsub <- X[[i]] beta <- get_species_b(b, object, i) lp <- Xsub[[1]] %*% beta[[1]] @@ -349,7 +361,7 @@ calc_dependent_abundance <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3r add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] lp <- lp + add } - ests[[i]] <- exp(lp) + ests[[i]] <- ilink(lp) post[[i]][,k] <- ests[[i]] } } From 2d9f18336bef53e3dd542ab9e60960d7deed1ad3 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 14 Jul 2024 20:17:18 -0400 Subject: [PATCH 08/21] Fix some bugs in input processing --- R/occuRNMulti.R | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 4ff5f7e0..30959541 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -135,7 +135,11 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", # Obs covs sc_long_ind <- rep(1:M, each=J) sc_long <- sc[sc_long_ind,] - oc <- cbind(sc_long, umf@obsCovs) + if(is.null(umf@obsCovs)){ + oc <- sc_long + } else { + oc <- cbind(sc_long, umf@obsCovs) + } vv <- lapply(detformulas, function(x){ model.matrix(x, oc) @@ -169,7 +173,7 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", det_n <- sapply(vv, ncol) det_end <- cumsum(det_n) + max(state_end) det_start <- det_end - det_n + 1 - det_ind <- cbind(det_end, det_start) + det_ind <- cbind(det_start, det_end) det_rng <- min(det_ind):max(det_ind) From 65d627c432631512eee8b09822e13bf0d624f287 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 15 Jul 2024 08:00:05 -0400 Subject: [PATCH 09/21] Easier specification of species to model occupancy for --- R/occuRNMulti.R | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 30959541..7b3c6abc 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -17,23 +17,20 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, S <- length(stateformulas) dep <- create_dep_matrix(stateformulas) - if(missing(modelOccupancy)){ - modelOccupancy <- rep(0, S) - names(modelOccupancy) <- names(stateformulas) - } else { - stopifnot(identical(names(modelOccupancy), names(stateformulas))) - if(is.logical(modelOccupancy)){ - modelOccupancy <- modelOccupancy*1 - } + modocc <- rep(0, S) + names(modocc) <- names(data@ylist) + if(!missing(modelOccupancy)){ + stopifnot(all(modelOccupancy %in% names(modocc))) + modocc[modelOccupancy] <- 1 not_allowed_occ <- colSums(dep) - sp_mismatch <- modelOccupancy & not_allowed_occ + sp_mismatch <- modocc & not_allowed_occ if(any(sp_mismatch)){ stop(paste("Species", - paste(names(modelOccupancy)[sp_mismatch], collapse=", "), + paste(names(modocc)[sp_mismatch], collapse=", "), "cannot use occupancy model."), call.=FALSE) } } - + data <- as(data, "unmarkedFrameOccuRNMulti") gd <- getDesign(data, detformulas, stateformulas) @@ -41,10 +38,10 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, nP <- max(gd$det_ind) - has_occ <- sum(modelOccupancy) > 0 + has_occ <- sum(modocc) > 0 if(has_occ){ - sp_occ <- names(modelOccupancy)[modelOccupancy == 1] - sp_abun <- names(modelOccupancy)[modelOccupancy == 0] + sp_occ <- names(modocc)[modocc == 1] + sp_abun <- names(modocc)[modocc == 0] occ_ind <- gd$state_ind[grepl(paste(paste0("^",sp_occ), collapse="|"), rownames(gd$state_ind)),,drop=FALSE] occ_rng <- min(occ_ind):max(occ_ind) @@ -61,7 +58,7 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, names(starts) <- gd$par_names nll_C <- function(pars){ - nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, modelOccupancy, + nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, modocc, gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads) } @@ -105,7 +102,7 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(), detformulas = detformulas, stateformulas = stateformulas, - modelOccupancy = modelOccupancy, + modelOccupancy = modocc, formula = ~1, data = data, #sitesRemoved = designMats$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, From a29d6981c0662355edc35e982e33e9fce20b1d29 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 15 Jul 2024 08:23:46 -0400 Subject: [PATCH 10/21] Make sure K length is correct --- R/occuRNMulti.R | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 7b3c6abc..3bc21e8a 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -34,6 +34,10 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, data <- as(data, "unmarkedFrameOccuRNMulti") gd <- getDesign(data, detformulas, stateformulas) + if(length(K) == 1){ + K <- rep(K, S) + } + stopifnot(length(K) == S) Kmin <- sapply(gd$ylist, function(x) apply(x, 1, max)) nP <- max(gd$det_ind) From d41d5f081b739aa37eba275850424511a74e4a46 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 15 Jul 2024 11:34:14 -0400 Subject: [PATCH 11/21] Handle missing values in y --- R/RcppExports.R | 4 ++-- R/occuRNMulti.R | 17 +++++++++++--- src/RcppExports.cpp | 14 +++++++----- src/nll_occuRNMulti.cpp | 49 +++++++++++++++++++++++------------------ 4 files changed, 52 insertions(+), 32 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 22fd89db..bb79441a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -73,8 +73,8 @@ nll_occuRN <- function(beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, thre .Call(`_unmarked_nll_occuRN`, beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, threads) } -nll_occuRNMulti <- function(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads) { - .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads) +nll_occuRNMulti <- function(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, miss, site_miss, threads) { + .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, miss, site_miss, threads) } nll_occuTTD <- function(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) { diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 3bc21e8a..1faa2dee 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -38,7 +38,12 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, K <- rep(K, S) } stopifnot(length(K) == S) - Kmin <- sapply(gd$ylist, function(x) apply(x, 1, max)) + Kmin <- sapply(gd$ylist, function(x){ + apply(x, 1, function(y){ + if(all(is.na(y))) return(0) + max(y, na.rm=TRUE) + }) + }) nP <- max(gd$det_ind) @@ -63,7 +68,8 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, nll_C <- function(pars){ nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, modocc, - gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads) + gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, + gd$miss, gd$site_miss, threads) } fm <- optim(starts, nll_C, method=method, hessian=se, ...) @@ -184,8 +190,13 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", det_par <- paste(det_prefix, det_par) par_names <- c(state_par, det_par) + # Missing values + miss <- is.na(umf@ylist[[1]]) * 1 + site_miss <- apply(umf@ylist[[1]], 1, function(x) all(is.na(x))) * 1 + list(par_names = par_names, state_ind = state_ind, det_ind = det_ind, - ylist = ylist, state_dm = mm, det_dm = vv) + ylist = ylist, state_dm = mm, det_dm = vv, + miss = miss, site_miss = site_miss) }) create_dep_matrix <- function(stateformulas){ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 455376b7..345a9a2e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -471,8 +471,8 @@ BEGIN_RCPP END_RCPP } // nll_occuRNMulti -double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, arma::ivec modOcc, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads); -RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP modOccSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, arma::ivec modOcc, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, arma::imat miss, arma::ivec site_miss, int threads); +RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP modOccSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP missSEXP, SEXP site_missSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -487,8 +487,10 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::imat >::type dep(depSEXP); Rcpp::traits::input_parameter< arma::ivec >::type K(KSEXP); Rcpp::traits::input_parameter< arma::imat >::type Kmin(KminSEXP); + Rcpp::traits::input_parameter< arma::imat >::type miss(missSEXP); + Rcpp::traits::input_parameter< arma::ivec >::type site_miss(site_missSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads)); + rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, miss, site_miss, threads)); return rcpp_result_gen; END_RCPP } @@ -586,8 +588,8 @@ BEGIN_RCPP END_RCPP } -RcppExport SEXP getDetVecs(void *, void *, void *, void *, void *); -RcppExport SEXP getSingleDetVec(void *, void *, void *); +RcppExport SEXP getDetVecs(SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP getSingleDetVec(SEXP, SEXP, SEXP); static const R_CallMethodDef CallEntries[] = { {"_unmarked_get_lik_trans", (DL_FUNC) &_unmarked_get_lik_trans, 2}, @@ -608,7 +610,7 @@ static const R_CallMethodDef CallEntries[] = { {"_unmarked_nll_occuMulti", (DL_FUNC) &_unmarked_nll_occuMulti, 15}, {"_unmarked_nll_occuPEN", (DL_FUNC) &_unmarked_nll_occuPEN, 11}, {"_unmarked_nll_occuRN", (DL_FUNC) &_unmarked_nll_occuRN, 10}, - {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 12}, + {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 14}, {"_unmarked_nll_occuTTD", (DL_FUNC) &_unmarked_nll_occuTTD, 17}, {"_unmarked_nll_pcount", (DL_FUNC) &_unmarked_nll_pcount, 11}, {"_unmarked_nll_pcountOpen", (DL_FUNC) &_unmarked_nll_pcountOpen, 35}, diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp index 0de39b96..f1029400 100644 --- a/src/nll_occuRNMulti.cpp +++ b/src/nll_occuRNMulti.cpp @@ -9,11 +9,13 @@ using namespace Rcpp; using namespace arma; //Calculate part of abundance likelihood to integrate over (abun f[g] * detect [g]) -double calc_fg(int k, double lam, int J, arma::vec r, arma::rowvec y){ +double calc_fg(int k, double lam, int J, arma::vec r, arma::rowvec y, + arma::irowvec miss){ double f = R::dpois(k, lam, 0); double p; double g = 0; for (int j = 0; j sp2 @@ -156,9 +163,9 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i lik_sp2 = 0; par2_m = lp2(m) + k1 * lp2_sp1(m); if(modOcc(1)){ - lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m)); + lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m), miss.row(m)); } else { - lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m)); + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m), miss.row(m)); } lik_m += fg1 * lik_sp2; @@ -170,14 +177,14 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i par2_m = exp(lp2(m) + k1 * lp2_sp1(m)); for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ - fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m)); + fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m), miss.row(m)); lik_sp3 = 0; par3_m = lp3(m) + k2 * lp3_sp2(m); if(modOcc(2)){ - lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } else { - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } lik_sp2 += fg2 * lik_sp3; } @@ -191,14 +198,14 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ - fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(m)); + fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(m), miss.row(m)); par3_m = lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m); lik_sp3 = 0; if(modOcc(2)){ - lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } else { - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } lik_sp2 += fg2 * lik_sp3; } @@ -212,17 +219,17 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i par2_m = lp2(m) + k1 * lp2_sp1(m); lik_sp2 = 0; if(modOcc(1)){ - lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m)); + lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m), miss.row(m)); } else { - lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m)); + lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m), miss.row(m)); } par3_m = lp3(m) + k1 * lp3_sp1(m); lik_sp3 = 0; if(modOcc(2)){ - lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } else { - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } lik_m += fg1 * lik_sp2 * lik_sp3; @@ -233,14 +240,14 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i par2_m = exp(lp2(m) + k1 * lp2_sp1(m)); for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){ - fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m)); + fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m), miss.row(m)); lik_sp3 = 0; par3_m = lp3(m) + k1 * lp2_sp1(m) + k2 * lp3_sp2(m); if(modOcc(2)){ - lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } else { - lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m)); + lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m), miss.row(m)); } lik_sp2 += fg2 * lik_sp3; } From ab39016f6018db0328846716d077da93ddd03a80 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 15 Jul 2024 11:51:25 -0400 Subject: [PATCH 12/21] Handle missing values in covariates --- R/occuRNMulti.R | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 1faa2dee..bdab5354 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -134,10 +134,20 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", mm <- lapply(stateformulas, function(s){ if(!is.list(s)){ - return(model.matrix(s, sc)) + mf <- model.frame(s, sc, na.action=stats::na.pass) + return(model.matrix(s, mf)) } - lapply(s, function(x) model.matrix(x, sc)) + lapply(s, function(x){ + mf <- model.frame(x, sc, na.action=stats::na.pass) + model.matrix(x, mf) + }) }) + all_mm <- lapply(mm, function(x){ + if(is.matrix(x)) return(x) + do.call(cbind, x) + }) + all_mm <- do.call(cbind, all_mm) + miss_site_covs <- apply(all_mm, 1, function(x) any(is.na(x))) # Obs covs sc_long_ind <- rep(1:M, each=J) @@ -149,8 +159,12 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", } vv <- lapply(detformulas, function(x){ - model.matrix(x, oc) + mf <- model.frame(x, oc, na.action=stats::na.pass) + model.matrix(x, mf) }) + all_vv <- do.call(cbind, vv) + miss_obs_covs <- apply(all_vv, 1, function(x) any(is.na(x))) + miss_obs_covs <- matrix(miss_obs_covs, M, J, byrow=TRUE) # Indices state_n <- unlist(lapply(mm, function(x){ @@ -191,8 +205,8 @@ setMethod("getDesign", "unmarkedFrameOccuRNMulti", par_names <- c(state_par, det_par) # Missing values - miss <- is.na(umf@ylist[[1]]) * 1 - site_miss <- apply(umf@ylist[[1]], 1, function(x) all(is.na(x))) * 1 + miss <- (is.na(umf@ylist[[1]]) | miss_obs_covs) * 1 + site_miss <- (apply(umf@ylist[[1]], 1, function(x) all(is.na(x))) | miss_site_covs) * 1 list(par_names = par_names, state_ind = state_ind, det_ind = det_ind, ylist = ylist, state_dm = mm, det_dm = vv, From 34e65448f4f08dcc8223a8e573b33335125c3c63 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 15 Jul 2024 13:10:02 -0400 Subject: [PATCH 13/21] For now, make sure detformulas are also in same order as y --- R/occuRNMulti.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index bdab5354..f03409a2 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -13,6 +13,7 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, starts, method="BFGS", se = TRUE, threads=1, ...){ stopifnot(all(names(data@ylist) == names(stateformulas))) + stopifnot(all(names(data@ylist) == names(detformulas))) S <- length(stateformulas) dep <- create_dep_matrix(stateformulas) From 66ceeae80b179804148973820721c77a1145dd71 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 17 Jul 2024 10:33:36 -0400 Subject: [PATCH 14/21] Handle out-of-order formula/data lists --- R/occuRNMulti.R | 87 +++++++++++++++++++++++++++++++++---------------- 1 file changed, 59 insertions(+), 28 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index f03409a2..c3ae6252 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -12,12 +12,21 @@ occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, K = rep(10, length(stateformulas)), starts, method="BFGS", se = TRUE, threads=1, ...){ - stopifnot(all(names(data@ylist) == names(stateformulas))) - stopifnot(all(names(data@ylist) == names(detformulas))) - S <- length(stateformulas) dep <- create_dep_matrix(stateformulas) + # Reorder everything to match + stopifnot(all(names(data@ylist) %in% names(stateformulas))) + stopifnot(all(names(data@ylist) %in% names(detformulas))) + species_order <- colnames(dep) + stateformulas <- stateformulas[species_order] + stateformulas <- lapply(stateformulas, function(x){ + if(is.call(x)) return(x) + x[order(match(names(x), c("", species_order)))] + }) + detformulas <- detformulas[species_order] + data@ylist <- data@ylist[species_order] + modocc <- rep(0, S) names(modocc) <- names(data@ylist) if(!missing(modelOccupancy)){ @@ -218,34 +227,56 @@ create_dep_matrix <- function(stateformulas){ S <- length(stateformulas) dep <- matrix(0, S, S) rownames(dep) <- colnames(dep) <- names(stateformulas) - allowed_species <- character(0) - - # Species 1 - stopifnot(is.call(stateformulas[[1]])) - allowed_species <- c(allowed_species, names(stateformulas)[1]) - -# Species 2 - if(is.list(stateformulas[[2]])){ - dep_sp <- names(stateformulas[[2]])[2:length(stateformulas[[2]])] - stopifnot(all(dep_sp %in% allowed_species)) - for (i in dep_sp){ - dep[names(stateformulas)[2], dep_sp] <- 1 - } + + # Dominant species + is_dom <- sapply(stateformulas, is.call) + if(sum(is_dom) == 0){ + stop("Must be at least one dominant species not depending on any others", + call.=FALSE) } - allowed_species <- c(allowed_species, names(stateformulas)[2]) - - # Species 3 - if(S > 2){ - if(is.list(stateformulas[[3]])){ - dep_sp <- names(stateformulas[[3]])[2:length(stateformulas[[3]])] - stopifnot(all(dep_sp %in% allowed_species)) - for (i in dep_sp){ - dep[names(stateformulas)[3], dep_sp] <- 1 - } - } - allowed_species <- c(allowed_species, names(stateformulas)[3]) + allowed_species <- sort(names(stateformulas)[is_dom]) + + # 2nd level + is_2nd <- sapply(stateformulas, function(x){ + if(is.call(x)) return(FALSE) + stopifnot(is.list(x)) + nm <- names(x)[names(x) != ""] + all(nm %in% allowed_species) + }) + species_2nd <- sort(names(stateformulas)[is_2nd]) + allowed_species <- c(allowed_species, species_2nd) + + for (i in species_2nd){ + dep_sp <- names(stateformulas[[i]]) + dep_sp <- dep_sp[dep_sp != ""] + dep[i, dep_sp] <- 1 } + # 3rd level + is_3rd <- sapply(1:length(stateformulas), function(i){ + if(names(stateformulas)[i] %in% allowed_species) return(FALSE) + x <- stateformulas[[i]] + stopifnot(is.list(x)) + nm <- names(x)[names(x) != ""] + all(nm %in% allowed_species) + }) + names(is_3rd) <- names(stateformulas) + + species_3rd <- sort(names(stateformulas)[is_3rd]) + allowed_species <- c(allowed_species, species_3rd) + + for (i in species_3rd){ + dep_sp <- names(stateformulas[[i]]) + dep_sp <- dep_sp[dep_sp != ""] + dep[i, dep_sp] <- 1 + } + + if(any(!names(stateformulas) %in% allowed_species)){ + stop("Invalid stateformula list", call.=FALSE) + } + + dep <- dep[allowed_species, allowed_species] + dep } From 492d65fe45ae389eafe54aa16d9fb1de7cb019a8 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 18 Jul 2024 09:17:26 -0400 Subject: [PATCH 15/21] Allow fixing dominant species abundance when predicting --- R/occuRNMulti.R | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index c3ae6252..349ddd7b 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -300,7 +300,7 @@ setMethod("predict", "unmarkedFitOccuRNMulti", sp_3rd <- names(lev2)[!lev2] if(missing(newdata) || is.null(newdata)){ - + newdata <- NULL X <- lapply(object@stateformulas, function(x){ if(is.list(x)){ lapply(x, function(z) model.matrix(z, object@data@siteCovs)) @@ -322,10 +322,10 @@ setMethod("predict", "unmarkedFitOccuRNMulti", coef_est <- matrix(coef(object), nrow=1) colnames(coef_est) <- names(coef(object)) - point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink) + point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) message('Bootstrapping confidence intervals with ',nsims,' samples') - post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink) + post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) # Summarize bootstrap cis <- lapply(post, function(x){ @@ -371,7 +371,8 @@ get_species_b <- function(b, object, sp){ out } -calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink){ +calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, + ilink, newdata){ nsims <- nrow(samps) all_species <- names(object@data@ylist) @@ -402,7 +403,15 @@ calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd lp <- Xsub[[1]] %*% beta[[1]] for (j in 2:length(Xsub)){ other_sp <- names(Xsub)[j] - add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] + + # Use fixed species abundances if they are in newdata + if(other_sp %in% names(newdata)){ + other_sp_abun <- newdata[[other_sp]] + } else { + other_sp_abun <- ests[[other_sp]] + } + + add <- Xsub[[other_sp]] %*% beta[[other_sp]] * other_sp_abun lp <- lp + add } ests[[i]] <- ilink(lp) @@ -420,7 +429,15 @@ calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd lp <- Xsub[[1]] %*% beta[[1]] for (j in 2:length(Xsub)){ other_sp <- names(Xsub)[j] - add <- Xsub[[other_sp]] %*% beta[[other_sp]] * ests[[other_sp]] + + # Use fixed species abundances if they are in newdata + if(other_sp %in% names(newdata)){ + other_sp_abun <- newdata[[other_sp]] + } else { + other_sp_abun <- ests[[other_sp]] + } + + add <- Xsub[[other_sp]] %*% beta[[other_sp]] * other_sp_abun lp <- lp + add } ests[[i]] <- ilink(lp) From c0c60504d701a29f1ff4734067b3e648c3ca7e98 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 25 Jul 2024 15:03:46 -0400 Subject: [PATCH 16/21] Swap lower/upper --- R/occuRNMulti.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 349ddd7b..4ded6319 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -331,8 +331,8 @@ setMethod("predict", "unmarkedFitOccuRNMulti", cis <- lapply(post, function(x){ data.frame( SE = apply(x, 1, sd, na.rm=TRUE), - upper = apply(x, 1, quantile, 0.025, na.rm=TRUE), - lower = apply(x, 1, quantile, 0.975, na.rm=TRUE) + lower = apply(x, 1, quantile, 0.025, na.rm=TRUE), + upper = apply(x, 1, quantile, 0.975, na.rm=TRUE) ) }) From 7b14ed70502c8aad0819bd9d3acfe23789e8f086 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 31 Jul 2024 10:11:33 -0400 Subject: [PATCH 17/21] Add detection prediction --- R/occuRNMulti.R | 122 +++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 110 insertions(+), 12 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 4ded6319..274cc69a 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -280,6 +280,15 @@ create_dep_matrix <- function(stateformulas){ dep } +# Compute linear combinations of estimates in unmarkedFit objects. +setMethod("linearComb", + signature(obj = "unmarkedFitOccuRNMulti", coefficients = "matrixOrVector"), + function(obj, coefficients, type, offset = NULL, re.form=NULL) +{ + stop("This method not supported for occuRNMulti fits, use predict instead.", + call.=FALSE) +}) + setMethod("predict", "unmarkedFitOccuRNMulti", function(object, type, species, newdata, level = 0.95, nsims = 100, ...){ @@ -289,7 +298,6 @@ setMethod("predict", "unmarkedFitOccuRNMulti", stopifnot(species %in% names(object@data@ylist)) } if(type == "state"){ - samps <- MASS::mvrnorm(nsims, mu=coef(object), Sigma = vcov(object)) depmat <- unmarked:::create_dep_matrix(object@stateformulas) @@ -324,17 +332,25 @@ setMethod("predict", "unmarkedFitOccuRNMulti", colnames(coef_est) <- names(coef(object)) point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) - message('Bootstrapping confidence intervals with ',nsims,' samples') - post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) - # Summarize bootstrap - cis <- lapply(post, function(x){ - data.frame( - SE = apply(x, 1, sd, na.rm=TRUE), - lower = apply(x, 1, quantile, 0.025, na.rm=TRUE), - upper = apply(x, 1, quantile, 0.975, na.rm=TRUE) - ) - }) + if(!is.null(level)){ + message('Bootstrapping confidence intervals with ',nsims,' samples') + samps <- MASS::mvrnorm(nsims, mu=coef(object), Sigma = vcov(object)) + post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) + + # Summarize bootstrap + cis <- lapply(post, function(x){ + data.frame( + SE = apply(x, 1, sd, na.rm=TRUE), + lower = apply(x, 1, quantile, (1-level)/2, na.rm=TRUE), + upper = apply(x, 1, quantile, 1 - (1-level)/2, na.rm=TRUE) + ) + }) + } else { + cis <- lapply(point_est, function(x){ + data.frame(SE = rep(NA, nrow(x)), lower=NA, upper=NA) + }) + } out <- mapply(function(x, y){ cbind(Predicted = x, y) @@ -345,7 +361,63 @@ setMethod("predict", "unmarkedFitOccuRNMulti", } return(out) } else if(type == "det"){ - stop("det not supported yet") + + det_species <- names(object@data@ylist) + if(!missing(species)){ + det_species <- species + } + + gd <- unmarked:::getDesign(object@data, object@detformulas, object@stateformulas) + od <- unmarked:::get_orig_data(object, "det") + chunk_size <- 70 + + if(missing(newdata)) newdata <- NULL + out <- lapply(det_species, function(i){ + form <- object@detformulas[[i]] + cf <- coef(object) + inds <- which(grepl(paste0("p([",i,"]"), names(cf), fixed=TRUE)) + new_est <- object@estimates@estimates$det + new_est@estimates <- cf[inds] + new_est@fixed <- 1:length(inds) + new_est@covMat <- vcov(object)[inds,inds,drop=FALSE] + new_est@covMatBS <- object@covMatBS[inds,inds,drop=FALSE] + + if(is.null(newdata)){ + X <- unmarked:::make_mod_matrix(form, od, newdata, re.form=NULL)$X + } else { + X <- gd$det_dm[[i]] + } + nr <- nrow(X) + sep <- rep(1:ceiling(nr/chunk_size), each=chunk_size, length.out=nr) + + x_chunk <- lapply(unique(sep), + function(i) as.matrix(X[sep==i,,drop=FALSE])) + + prmat <- lapply(x_chunk, function(x_i){ + has_na <- apply(x_i, 1, function(x_i) any(is.na(x_i))) + # Work around linearComb bug where there can't be NAs in inputs + x_i[has_na,] <- 0 + lc <- linearComb(new_est, x_i) + lc <- backTransform(lc) + out <- data.frame(Predicted=coef(lc), SE=NA, lower=NA, upper=NA) + if(!is.null(level)){ + se <- SE(lc) + ci <- confint(lc, level=level) + out$SE <- se + out$lower <- ci[,1] + out$upper <- ci[,2] + } + out[has_na,] <- NA + out + }) + prmat <- do.call(rbind, prmat) + rownames(prmat) <- NULL + prmat + }) + if(length(out) == 1) out <- out[[1]] + return(out) + } else { + stop("type must be state or det", call.=FALSE) } }) @@ -446,3 +518,29 @@ calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd } post } + + +setMethod("getP", "unmarkedFitOccuRNMulti", function(object) +{ + + ylist <- object@data@ylist + S <- length(ylist) + N <- nrow(ylist[[1]]) + + dm <- getDesign(object@data,object@detformulas,object@stateformulas, maxOrder=maxOrder) + pred <- predict(object,'det',se.fit=F) + dets <- do.call(cbind,lapply(pred,`[`,,1)) + #ugly mess + out <- list() + for (i in 1:S){ + pmat <- array(NA,dim(ylist[[1]])) + for (j in 1:N){ + ps <- dets[dm$yStart[j]:dm$yStop[j],i] + not_na <- !is.na(ylist[[i]][j,]) + pmat[j,not_na] <- ps + } + out[[i]] <- pmat + } + names(out) <- names(ylist) + out +}) From 7aeac0f8dc4be3ef35db2336333715499b87f6b1 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 31 Jul 2024 10:33:58 -0400 Subject: [PATCH 18/21] Add getP, fitted, residuals, SSE methods --- R/occuRNMulti.R | 76 +++++++++++++++++++++++++++++++++++++------------ 1 file changed, 58 insertions(+), 18 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 274cc69a..2987f8ba 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -293,7 +293,6 @@ setMethod("predict", "unmarkedFitOccuRNMulti", function(object, type, species, newdata, level = 0.95, nsims = 100, ...){ - stopifnot(type %in% c("state", "det")) if(!missing(species)){ stopifnot(species %in% names(object@data@ylist)) } @@ -382,7 +381,7 @@ setMethod("predict", "unmarkedFitOccuRNMulti", new_est@covMat <- vcov(object)[inds,inds,drop=FALSE] new_est@covMatBS <- object@covMatBS[inds,inds,drop=FALSE] - if(is.null(newdata)){ + if(!is.null(newdata)){ X <- unmarked:::make_mod_matrix(form, od, newdata, re.form=NULL)$X } else { X <- gd$det_dm[[i]] @@ -414,6 +413,7 @@ setMethod("predict", "unmarkedFitOccuRNMulti", rownames(prmat) <- NULL prmat }) + names(out) <- det_species if(length(out) == 1) out <- out[[1]] return(out) } else { @@ -525,22 +525,62 @@ setMethod("getP", "unmarkedFitOccuRNMulti", function(object) ylist <- object@data@ylist S <- length(ylist) - N <- nrow(ylist[[1]]) + M <- nrow(ylist[[1]]) + J <- ncol(ylist[[1]]) - dm <- getDesign(object@data,object@detformulas,object@stateformulas, maxOrder=maxOrder) - pred <- predict(object,'det',se.fit=F) - dets <- do.call(cbind,lapply(pred,`[`,,1)) - #ugly mess - out <- list() - for (i in 1:S){ - pmat <- array(NA,dim(ylist[[1]])) - for (j in 1:N){ - ps <- dets[dm$yStart[j]:dm$yStop[j],i] - not_na <- !is.na(ylist[[i]][j,]) - pmat[j,not_na] <- ps + pr <- predict(object, type = "det", level=NULL) + pr <- lapply(pr, function(x){ + matrix(x$Predicted, nrow=M, ncol=J, byrow=TRUE) + }) + pr +}) + +setMethod("fitted", "unmarkedFitOccuRNMulti", function(object, na.rm=FALSE){ + S <- length(object@data@ylist) + M <- nrow(object@data@ylist[[1]]) + J <- ncol(object@data@ylist[[1]]) + + # This will be lambda for RN models and psi for occupancy models + state <- predict(object, type="state", level=NULL) + # Reorganize into M x J matrices + state <- lapply(state, function(x){ + matrix(rep(x$Predicted, J), nrow=M, ncol=J) + }) + # This will be r for RN models and p for occupancy models + p <- getP(object) + + # Need to calculate this differently depending on if state is abundance or occupancy + fitted <- lapply(1:S, function(i){ + if(fit@modelOccupancy[i]){ + # If occupancy, fitted is psi * p + state[[i]] * p[[i]] + } else { + # If abundance + 1 - exp(-state[[i]]*p[[i]]) ## analytical integration. } - out[[i]] <- pmat - } - names(out) <- names(ylist) - out + }) + names(fitted) <- names(object@data@ylist) + fitted +}) + +setMethod("residuals", "unmarkedFitOccuRNMulti", function(object, ...){ + ylist <- object@data@ylist + S <- length(ylist) + ft <- fitted(object) + res <- lapply(1:S, function(i){ + ylist[[i]] - ft[[i]] + }) + names(res) <- names(ylist) + res +}) + +setMethod("SSE", "unmarkedFitOccuRNMulti", function(fit, ...){ + r <- do.call(rbind, residuals(fit)) + return(c(SSE = sum(r^2, na.rm=T))) +}) + +setMethod("nonparboot", "unmarkedFitOccuRNMulti", + function(object, B = 0, keepOldSamples = TRUE, ...) +{ + stop("Method not supported for this fit type", call.=FALSE) }) From 87daab300a858980ac0836141af0573726e9df07 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 31 Jul 2024 11:34:03 -0400 Subject: [PATCH 19/21] Add update, simulate, parboot methods --- R/occuRNMulti.R | 107 ++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 103 insertions(+), 4 deletions(-) diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 2987f8ba..3365e657 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -296,6 +296,12 @@ setMethod("predict", "unmarkedFitOccuRNMulti", if(!missing(species)){ stopifnot(species %in% names(object@data@ylist)) } + + se <- TRUE + if(is.null(hessian(object))){ + se = FALSE + } + if(type == "state"){ depmat <- unmarked:::create_dep_matrix(object@stateformulas) @@ -332,7 +338,7 @@ setMethod("predict", "unmarkedFitOccuRNMulti", point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) - if(!is.null(level)){ + if(se & !is.null(level)){ message('Bootstrapping confidence intervals with ',nsims,' samples') samps <- MASS::mvrnorm(nsims, mu=coef(object), Sigma = vcov(object)) post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) @@ -378,8 +384,13 @@ setMethod("predict", "unmarkedFitOccuRNMulti", new_est <- object@estimates@estimates$det new_est@estimates <- cf[inds] new_est@fixed <- 1:length(inds) - new_est@covMat <- vcov(object)[inds,inds,drop=FALSE] - new_est@covMatBS <- object@covMatBS[inds,inds,drop=FALSE] + if(se){ + new_est@covMat <- vcov(object)[inds,inds,drop=FALSE] + new_est@covMatBS <- object@covMatBS[inds,inds,drop=FALSE] + } else { + new_est@covMat <- matrix(NA, nrow=length(inds), ncol=length(inds)) + new_est@covMatBS <- matrix(NA, nrow=length(inds), ncol=length(inds)) + } if(!is.null(newdata)){ X <- unmarked:::make_mod_matrix(form, od, newdata, re.form=NULL)$X @@ -551,7 +562,7 @@ setMethod("fitted", "unmarkedFitOccuRNMulti", function(object, na.rm=FALSE){ # Need to calculate this differently depending on if state is abundance or occupancy fitted <- lapply(1:S, function(i){ - if(fit@modelOccupancy[i]){ + if(object@modelOccupancy[i]){ # If occupancy, fitted is psi * p state[[i]] * p[[i]] } else { @@ -584,3 +595,91 @@ setMethod("nonparboot", "unmarkedFitOccuRNMulti", { stop("Method not supported for this fit type", call.=FALSE) }) + +setMethod("simulate", "unmarkedFitOccuRNMulti", + function(object, nsim = 1, seed = NULL, na.rm=TRUE){ + + ylist <- object@data@ylist + M <- nrow(ylist[[1]]) + J <- ncol(ylist[[1]]) + S <- length(ylist) + occ <- object@modelOccupancy + + state <- predict(object, type="state", level=NULL) + state <- lapply(state, function(x) x$Predicted) + p <- getP(object) + + sims <- lapply(1:nsim, function(i){ + out <- lapply(1:S, function(s){ + y <- matrix(NA, M, J) + if(occ[s]){ + z <- rbinom(M, 1, state[[s]]) + for (j in 1:J){ + y[,j] <- rbinom(M, 1, p[[s]][,j]) + } + y + } else { + N <- rpois(M, state[[s]]) + for (j in 1:J){ + y[,j] <- rbinom(M, N, p[[s]][,j]) + } + y[y>0] <- 1 + y + } + }) + names(out) <- names(ylist) + out + }) + + sims +}) + +setMethod("replaceY", "unmarkedFrameOccuRNMulti", + function(object, newY, replNA=TRUE, ...){ + if(replNA){ + newY <- mapply(function(x, y){ is.na(x) <- is.na(y); x}, + newY , object@ylist, SIMPLIFY=FALSE) + } + object@ylist <- newY + object +}) + + +setMethod("update", "unmarkedFitOccuRNMulti", + function(object, detformulas, stateformulas, ..., evaluate = TRUE) +{ + + call <- object@call + if (is.null(call)) + stop("need an object with call slot") + if(!missing(detformulas)){ + call[["detformulas"]] <- detformulas + } else { + call[["detformulas"]] <- object@detformulas + } + if(!missing(stateformulas)){ + call[["stateformulas"]] <- stateformulas + } else { + call[["stateformulas"]] <- object@stateformulas + } + extras <- match.call(call=sys.call(-1), + expand.dots = FALSE)$... + if (length(extras) > 0) { + existing <- !is.na(match(names(extras), names(call))) + for (a in names(extras)[existing]) + call[[a]] <- extras[[a]] + if (any(!existing)) { + call <- c(as.list(call), extras[!existing]) + call <- as.call(call) + } + } + if (evaluate) + eval(call, parent.frame(2)) + else call +}) + +setMethod("ranef", "unmarkedFitOccuRNMulti", + function(object, B = 0, keepOldSamples = TRUE, ...) +{ + stop("Not currently supported for occuRNMulti", call.=FALSE) +}) From e6cbed5b450eab6253712811f452e34c99f8935b Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 31 Jul 2024 12:39:20 -0400 Subject: [PATCH 20/21] Add docs --- R/occuRNMulti.R | 6 +- man/SSE.Rd | 1 + man/fitted-methods.Rd | 1 + man/getP-methods.Rd | 1 + man/linearComb-methods.Rd | 3 +- man/nonparboot-methods.Rd | 2 +- man/occuRNMulti.Rd | 209 ++++++++++++++++++++++++++++++++++++++ man/predict-methods.Rd | 1 + man/ranef-methods.Rd | 1 + man/simulate-methods.Rd | 1 + man/unmarkedFit-class.Rd | 2 + 11 files changed, 223 insertions(+), 5 deletions(-) create mode 100644 man/occuRNMulti.Rd diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 3365e657..6029354c 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -335,13 +335,13 @@ setMethod("predict", "unmarkedFitOccuRNMulti", coef_est <- matrix(coef(object), nrow=1) colnames(coef_est) <- names(coef(object)) - point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) + point_est <- calc_dependent_response(coef_est, X, object, nr, sp_top, sp_2nd, sp_3rd, newdata) if(se & !is.null(level)){ message('Bootstrapping confidence intervals with ',nsims,' samples') samps <- MASS::mvrnorm(nsims, mu=coef(object), Sigma = vcov(object)) - post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, ilink, newdata) + post <- calc_dependent_response(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, newdata) # Summarize bootstrap cis <- lapply(post, function(x){ @@ -455,7 +455,7 @@ get_species_b <- function(b, object, sp){ } calc_dependent_response <- function(samps, X, object, nr, sp_top, sp_2nd, sp_3rd, - ilink, newdata){ + newdata){ nsims <- nrow(samps) all_species <- names(object@data@ylist) diff --git a/man/SSE.Rd b/man/SSE.Rd index e3204f2a..ce73cbec 100644 --- a/man/SSE.Rd +++ b/man/SSE.Rd @@ -5,6 +5,7 @@ \alias{SSE,unmarkedFitOccuMulti-method} \alias{SSE,unmarkedFitGDR-method} \alias{SSE,unmarkedFitIDS-method} +\alias{SSE,unmarkedFitOccuRNMulti-method} \title{Compute Sum of Squared Residuals for a Model Fit.} \description{ Compute the sum of squared residuals for an unmarked fit object. This diff --git a/man/fitted-methods.Rd b/man/fitted-methods.Rd index 9ead41cd..2a9439d3 100644 --- a/man/fitted-methods.Rd +++ b/man/fitted-methods.Rd @@ -18,6 +18,7 @@ \alias{fitted,unmarkedFitDailMadsen-method} \alias{fitted,unmarkedFitGOccu-method} \alias{fitted,unmarkedFitOccuCOP-method} +\alias{fitted,unmarkedFitOccuRNMulti-method} \title{Methods for Function fitted in Package `unmarked'} \description{Extracted fitted values from a fitted model. diff --git a/man/getP-methods.Rd b/man/getP-methods.Rd index 32028ab4..bd5c40c1 100644 --- a/man/getP-methods.Rd +++ b/man/getP-methods.Rd @@ -20,6 +20,7 @@ \alias{getP,unmarkedFitIDS-method} \alias{getP,unmarkedFitGOccu-method} \alias{getP,unmarkedFitOccuCOP-method} +\alias{getP,unmarkedFitOccuRNMulti-method} \title{Methods for Function getP in Package `unmarked'} \description{ diff --git a/man/linearComb-methods.Rd b/man/linearComb-methods.Rd index 83d3764b..9c103b92 100644 --- a/man/linearComb-methods.Rd +++ b/man/linearComb-methods.Rd @@ -4,6 +4,7 @@ \alias{linearComb-methods} \alias{linearComb,unmarkedEstimate,matrixOrVector-method} \alias{linearComb,unmarkedFit,matrixOrVector-method} +\alias{linearComb,unmarkedFitOccuRNMulti,matrixOrVector-method} \alias{show,unmarkedLinComb-method} \title{Methods for Function linearComb in Package `unmarked'} \description{ @@ -30,4 +31,4 @@ fm <- multinomPois(~ 1 ~ ufc + trba, ovenFrame) linearComb(fm, c(1, 0.5, 0.5), type = "state") linearComb(fm, matrix(c(1, 0.5, 0.5, 1, 0, 0, 1, 0, 0.5), 3, 3, byrow=TRUE), type="state") -} \ No newline at end of file +} diff --git a/man/nonparboot-methods.Rd b/man/nonparboot-methods.Rd index 53ff723f..796b87c0 100644 --- a/man/nonparboot-methods.Rd +++ b/man/nonparboot-methods.Rd @@ -20,7 +20,7 @@ \alias{nonparboot,unmarkedFitIDS-method} \alias{nonparboot,unmarkedFitDailMadsen-method} \alias{nonparboot,unmarkedFitOccuCOP-method} - +\alias{nonparboot,unmarkedFitOccuRNMulti-method} \title{ Nonparametric bootstrapping in unmarked } \description{ diff --git a/man/occuRNMulti.Rd b/man/occuRNMulti.Rd new file mode 100644 index 00000000..9fff3a6e --- /dev/null +++ b/man/occuRNMulti.Rd @@ -0,0 +1,209 @@ +\name{occuRNMulti} + +\alias{occuRNMulti} + +\title{Fit the Twining et al. (2024) Abundance-mediated Species Interaction Model} + +\usage{occuRNMulti(detformulas, stateformulas, data, modelOccupancy, + K = rep(10, length(stateformulas)), + starts, method="BFGS", se=TRUE, threads = 1, ...)} + +\arguments{ + \item{detformulas}{Named list of formulas for the detection models, one per species.} + \item{stateformulas}{A named list of formulas for the state model. The list should + have one element per species. Each list element is a single formula (for a dominant species) + or a list of formulas (for a subordinate species). In the latter case the first formula + corresponds to the part of the linear predictor unaffected by interactions, + and each subsequent formula corresponding to the linear predictor associated + with interactions with another species. See examples.} + \item{data}{An \code{\link{unmarkedFrameOccuMulti}} object.} + \item{modelOccupancy}{An (optional) character vector of species names, + corresponding to species for which you want to model occupancy as the + state parameter instead of abundance. This is only possible for + subordinate species that do not have an effect on any other species.} + \item{K}{Vector of upper summation indices (one per species) used to + numerically integrate out the latent abundance. This should be set high + enough so that it does not affect the parameter estimates. Computation + time will increase with K.} + \item{starts}{Vector of parameter starting values.} + \item{method}{Optimization method used by \code{\link{optim}}.} + \item{se}{Logical specifying whether or not to compute standard + errors.} + \item{threads}{Set the number of threads to use for optimization in C++, if + OpenMP is available on your system. Increasing the number of threads + may speed up optimization in some cases by running the likelihood + calculation in parallel. If \code{threads=1} (the default), OpenMP is disabled.} + \item{\dots}{Additional arguments to optim, such as lower and upper + bounds} + } + +\description{Fit the Twining et al. (2024) abundance-mediated species interaction model + using multi-species detection/non-detection data. The model uses the Royle-Nichols + model (Royle and Nichols 2003) to estimate abundance from detection/non-detection data. + The abundance of a dominant species can affect the abundance of subordinate + species, and this interaction can be mediated by covariates. A maximum + of 3 species are allowed.} + +\details{ + +See \code{\link{unmarkedFrame}} and \code{\link{unmarkedFrameOccuMulti}} for a +description of how to supply data to the \code{data} argument. + +Royle and Nichols (2003) described how species abundance at a site can be +obtained from detection/non-detection data by exploiting the expected positive correlation +between number of individuals of a species present at a site \eqn{i} (\eqn{N_i}) +and overall detection probability of at least one individual (\eqn{p_i}): + +\deqn{p_i = 1 - (1 - r_i)^{N_i}} + +where \eqn{r_i} is the probability of detecting an individual animal. Abundance +is then typically assumed to come from a Poisson distribution with parameter +\eqn{lambda_i} which can be a function of covariates. + +Twining et al. (2024) extend this model to multi-species data and additionally +model interaction(s) between a dominant species (\eqn{D}) and a subordinate +species (\eqn{S}) by including parameter(s) in the linear predictor for \eqn{\lambda_S} +that are interacted with \eqn{N_D}. For example, expected abundance of the subordinate +species can include an interaction that is dependent on some covariate \eqn{x}: + +\deqn{\mathrm{log}(\lambda_{S,i}) = \beta_0 + (\gamma_0 + \gamma_x \cdot x_i) \cdot N_{D,i}} + +where the \eqn{gamma} terms are the interaction parameters. Note that the +subordinate species \eqn{S} could not itself affect abundance of dominant species +\eqn{D}; thus the interactions are directional and not symmetric as with the +Rota et al. (2016) model. Subordinate species could have occupancy models instead of +abundance models (controlled with function argument \code{modelOccupancy}), +in which case the parameter of interest is \eqn{\psi_{S}}. +However, any species for which you want to model an effect on another species +(i.e., a dominant species), you cannot fit an occupancy model. + +In principle this model can work with any number of of species and interactions. +\code{occuRNMulti} handles only a maximum of 3 species. However, +it supports all possible interaction combinations for 2 or 3 species. +The interactions are specified using the \code{stateformulas} argument. +See the examples for details. + +Detection models for each species are independent and potentially include different +covariates. Detection models are provided to the \code{detformulas} argument, +which takes a named list of individual formulas with length equal to +the number of species. +} + +\value{unmarkedFitOccuRNMulti object describing the model fit.} + +\references{ + Twining 2024 +} + +\author{Ken Kellner \email{contact@kenkellner.com}} + +\seealso{\code{\link{unmarkedFrameOccuMulti}}, \code{\link{occuRN}}} + + +\examples{ + +\donttest{ + +# A model where a dominant species (coyote) affects an intermediate species +# (marten) which affects a subordinate species (fisher) +# e.g., sp1 --> sp2 --> sp3 +set.seed(1) + +M <- 200 +J <- 5 + +N1 <- N2 <- N3 <- 25 +while((max(N1) > 20) | (max(N2) > 20) | (max(N3) > 20)){ +# Covariates +sc <- data.frame(x1 = rnorm(M), x2 = rnorm(M)) +oc <- data.frame(x3 = rnorm(M*J)) + +# Coyote expected abundance +l1 <- 2 + +# True coyote abundance +N1 <- rpois(M, l1) + +# Lambda for marten +l2 <- exp(log(1.3) - 0.3*N1 + 0.5*sc$x1*N1) + +# True marten abundance +N2 <- rpois(M, l2) + +# Lambda for fisher +l3 <- exp(log(1) + -0.2 * N2) +N3 <- rpois(M, l3) +} + +# True detection prob for individuals of both species +rs <- 0.2 +ps1 <- 1 - (1-rs)^N1 +ps2 <- 1 - (1-rs)^N2 +ps3 <- 1 - (1-rs)^N3 + +# True parameter values +truth <- c(log(2), log(1.3), -0.3, 0.5, log(1), -0.2, rep(log(0.2/(1-0.2)), 3)) + +# Simulate observations +y1 <- y2 <- y3 <- matrix(NA, M, J) +for (m in 1:M){ + for (j in 1:J){ + y1[m,] <- rbinom(J, 1, ps1[m]) + y2[m,] <- rbinom(J, 1, ps2[m]) + y3[m,] <- rbinom(J, 1, ps3[m]) + } +} + +umf <- unmarkedFrameOccuMulti(y=list(coyote=y1, marten=y2, fisher=y3), + siteCovs=sc, obsCovs=oc) + +# Specify formulas for model as named list +# sp1 --> sp2 --> sp3 +# Coyote is unaffected by other species, so it includes just a single formula +# Marten is affected by coyote, so it has a list where the first formula +# is the part of the linear predictor not affected by coyote (here intercept-only), +# and the second (named) list element is the coyote effect, which is mediated +# by covariate x1. +# Finally, fisher has an intercept-only model for the independent part, +# and is affected by marten abundance, but without any covariates. +sf <- list(coyote = ~1, + marten = list(~1, coyote = ~x1), + fisher = list(~1, marten = ~1)) + +df <- list(coyote = ~1, marten = ~1, fisher=~1) + +K <- c(max(N1) + 3, max(N2) + 3, max(N3) + 3) + +fit <- occuRNMulti(df, sf, umf, K=K, threads=4) +summary(fit) + +cbind(unm=coef(fit), truth=truth) + +# Predicted abundance for each species at each site +# This plugs in expected abundance (lambda) of dominant species +# into the linear predictor, rather than using the latent abundance N +pr <- predict(fit, type='state') +lapply(pr, head) + +# Detection probability for each species +pr <- predict(fit, type = "state", level=NULL) + +# Example formula structures for other models +# (could add covariates to any formula) +# sp1 --> sp2 +sf <- list(coyote = ~1, marten = list(~1, coyote = ~1)) +# sp2 <-- sp1 --> sp3 +sf <- list(coyote = ~1, + marten = list(~1, coyote = ~1), + fisher = list(~1, coyote = ~1)) +# sp1 --> sp3 <-- sp2 +sf <- list(coyote = ~1, + marten = ~1, + fisher = list(~1, coyote = ~1, marten = ~1)) +# sp1 --> sp2 --> sp3 AND sp1 --> sp3 +sf <- list(coyote = ~1, + marten = list(~1, coyote = ~1), + fisher = list(~1, coyote = ~1, marten = ~1)) +} + +} diff --git a/man/predict-methods.Rd b/man/predict-methods.Rd index 79fa0196..bc02c48e 100644 --- a/man/predict-methods.Rd +++ b/man/predict-methods.Rd @@ -17,6 +17,7 @@ \alias{predict,unmarkedFitDSO-method} \alias{predict,unmarkedFitGDR-method} \alias{predict,unmarkedFitIDS-method} +\alias{predict,unmarkedFitOccuRNMulti-method} \alias{predict,unmarkedFitList-method} \alias{predict,unmarkedRanef-method} \title{ Methods for Function predict in Package `unmarked' } diff --git a/man/ranef-methods.Rd b/man/ranef-methods.Rd index 66452d79..3a5a0f09 100644 --- a/man/ranef-methods.Rd +++ b/man/ranef-methods.Rd @@ -23,6 +23,7 @@ \alias{ranef,unmarkedFitGOccu-method} \alias{ranef,unmarkedFitOccuCOP-method} \alias{ranef,unmarkedFitIDS-method} +\alias{ranef,unmarkedFitOccuRNMulti-method} \title{ Methods for Function \code{ranef} in Package \pkg{unmarked} } \description{ Estimate posterior distributions of the random variables (latent diff --git a/man/simulate-methods.Rd b/man/simulate-methods.Rd index 073686cb..4804ddbb 100644 --- a/man/simulate-methods.Rd +++ b/man/simulate-methods.Rd @@ -23,6 +23,7 @@ \alias{simulate,unmarkedFitDailMadsen-method} \alias{simulate,unmarkedFitGOccu-method} \alias{simulate,unmarkedFitOccuCOP-method} +\alias{simulate,unmarkedFitOccuRNMulti-method} \alias{simulate,character-method} \title{Methods for Function simulate in Package `unmarked'} diff --git a/man/unmarkedFit-class.Rd b/man/unmarkedFit-class.Rd index 8237311d..5422b153 100644 --- a/man/unmarkedFit-class.Rd +++ b/man/unmarkedFit-class.Rd @@ -30,6 +30,7 @@ \alias{residuals,unmarkedFitGDR-method} \alias{residuals,unmarkedFitIDS-method} \alias{residuals,unmarkedFitOccuCOP-method} +\alias{residuals,unmarkedFitOccuRNMulti-method} \alias{update,unmarkedFit-method} \alias{update,unmarkedFitColExt-method} \alias{update,unmarkedFitGMM-method} @@ -41,6 +42,7 @@ \alias{update,unmarkedFitIDS-method} \alias{update,unmarkedFitDailMadsen-method} \alias{update,unmarkedFitGOccu-method} +\alias{update,unmarkedFitOccuRNMulti-method} \alias{sampleSize} \alias{sampleSize,unmarkedFit-method} \alias{unmarkedFitOccu-class} From 1f2b78b72a80d06528a501d4b91f6b6adeecc3c0 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 31 Jul 2024 13:41:45 -0400 Subject: [PATCH 21/21] Tests --- R/occuRNMulti.R | 8 +- tests/testthat/test_occuRNMulti.R | 164 ++++++++++++++++++++++++++++++ 2 files changed, 170 insertions(+), 2 deletions(-) create mode 100644 tests/testthat/test_occuRNMulti.R diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R index 6029354c..a7a38a56 100644 --- a/R/occuRNMulti.R +++ b/R/occuRNMulti.R @@ -316,9 +316,13 @@ setMethod("predict", "unmarkedFitOccuRNMulti", newdata <- NULL X <- lapply(object@stateformulas, function(x){ if(is.list(x)){ - lapply(x, function(z) model.matrix(z, object@data@siteCovs)) + lapply(x, function(z){ + mf <- model.frame(z, object@data@siteCovs, na.action=stats::na.pass) + model.matrix(z, mf) + }) } else { - model.matrix(x, object@data@siteCovs) + mf <- model.frame(x, object@data@siteCovs, na.action=stats::na.pass) + model.matrix(x, mf) } }) nr <- nrow(object@data@siteCovs) diff --git a/tests/testthat/test_occuRNMulti.R b/tests/testthat/test_occuRNMulti.R new file mode 100644 index 00000000..c08c48b9 --- /dev/null +++ b/tests/testthat/test_occuRNMulti.R @@ -0,0 +1,164 @@ +context("occuRNMulti fitting function") +skip_on_cran() + +# Simulate data +set.seed(1) +M <- 300 +J <- 5 +N1 <- N2 <- N3 <- 25 +while((max(N1) > 20) | (max(N2) > 20) | (max(N3) > 20)){ +# Covariates +sc <- data.frame(x1 = rnorm(M), x2 = rnorm(M)) +oc <- data.frame(x3 = rnorm(M*J)) +# Coyote expected abundance +l1 <- 2 +# True coyote abundance +N1 <- rpois(M, l1) +# Lambda for marten +l2 <- exp(log(1.3) - 0.3*N1 + 0.5*sc$x1*N1) +# True marten abundance +N2 <- rpois(M, l2) +# Lambda for fisher +l3 <- exp(log(1) + -0.2 * N2) +N3 <- rpois(M, l3) +} +# True detection prob for individuals of both species +rs <- 0.2 +ps1 <- 1 - (1-rs)^N1 +ps2 <- 1 - (1-rs)^N2 +ps3 <- 1 - (1-rs)^N3 +# True parameter values +truth <- c(log(2), log(1.3), -0.3, 0.5, log(1), -0.2, rep(log(0.2/(1-0.2)), 3)) +# Simulate observations +y1 <- y2 <- y3 <- matrix(NA, M, J) +for (m in 1:M){ + for (j in 1:J){ + y1[m,] <- rbinom(J, 1, ps1[m]) + y2[m,] <- rbinom(J, 1, ps2[m]) + y3[m,] <- rbinom(J, 1, ps3[m]) + } +} +umf <- unmarkedFrameOccuMulti(y=list(coyote=y1, marten=y2, fisher=y3), + siteCovs=sc, obsCovs=oc) + +test_that("occuRNMulti can fit 2-species models",{ + umf2 <- umf + umf2@ylist <- umf@ylist[1:2] + sf <- list(coyote = ~1, + marten = list(~1, coyote = ~x1)) + df <- list(coyote=~1, marten=~1) + fit <- occuRNMulti(df, sf, umf2[1:50,]) + expect_is(fit, "unmarkedFitOccuRNMulti") + expect_equivalent(coef(fit), c(0.3888, 0.4018,-0.4229,0.6316,-1.0916,-1.8663), + tol=1e-4) + + pr <- suppressMessages(predict(fit, type='state')) + expect_equal(length(pr), 2) + expect_equal(nrow(pr[[1]]), 50) + + nd <- data.frame(x1=c(0,1)) + pr <- suppressMessages(predict(fit, type='state', newdata=nd)) + expect_equal(pr[[2]]$Predicted, c(0.8009, 2.0334), tol=1e-4) + + res <- residuals(fit) + expect_equal(length(res), 2) + expect_equal(res[[1]][1], -0.30977, tol=1e-4) + + s <- simulate(fit, nsim=2) + expect_equal(length(s), 2) + expect_equal(length(s[[1]]), 2) + + pb <- parboot(fit, nsim=2) + expect_is(pb, "parboot") + + # Unsupported methods + expect_error(ranef(fit)) + expect_error(nonparboot(fit)) + expect_error(linearComb(fit, type='state', c(0,1,1))) + + # Check with NAs + umf2@ylist[[1]][1,1] <- NA + umf2@ylist[[2]][1,1] <- NA + umf2@ylist[[1]][2,] <- NA + umf2@ylist[[2]][2,] <- NA + umf2@siteCovs$x1[3] <- NA # missing covariate value + + fit_na <- occuRNMulti(df, sf, umf2[1:50,]) + expect_equal(coef(fit), coef(fit_na), tol=0.2) + + pr <- suppressMessages(predict(fit_na, type='state')) + expect_true(is.na(pr[[2]]$Predicted[3])) + + df2 <- list(coyote=~x1, marten=~1) + fit_na2 <- occuRNMulti(df2, sf, umf2[1:20,]) + pr <- predict(fit_na2, type='det') + expect_true(all(is.na(pr$coyote$Predicted[11:15]))) + + # Occupancy model for subordinate + fit_occ <- occuRNMulti(df, sf, umf2, modelOccupancy='marten') + expect_equal(names(fit_occ), c("lam", "psi", "det")) + pr <- suppressMessages(predict(fit_occ, type='state')) + expect_true(max(pr$marten$Predicted, na.rm=TRUE) < 1) +}) + +test_that("occuRNMulti can fit 3-species models", { + + # `sp1 --> sp2 --> sp3` + sf <- list(coyote = ~1, + marten = list(~1, coyote = ~x1), + fisher = list(~1, marten = ~1)) + df <- list(coyote=~1, marten=~1, fisher=~1) + + fit <- occuRNMulti(df, sf, umf[1:10,]) + expect_equal(length(coef(fit)), 9) + pr <- suppressMessages(predict(fit, 'state')) + expect_equal(length(pr), 3) + + # `sp2 <-- sp1 --> sp3` + sf <- list(coyote = ~1, + marten = list(~1, coyote = ~x1), + fisher = list(~1, coyote = ~1)) + + fit <- occuRNMulti(df, sf, umf[1:10,]) + expect_equal(length(coef(fit)), 9) + expect_equal( + names(coef(fit, type='lam')), + c("lam([coyote] (Intercept))", "lam([fisher] (Intercept))", + "lam([fisher:coyote] (Intercept))", + "lam([marten] (Intercept))", "lam([marten:coyote] (Intercept))", + "lam([marten:coyote] x1)") + ) + pr <- suppressMessages(predict(fit, 'state')) + expect_equal(length(pr), 3) + + # `sp1 --> sp3 <-- sp2` + sf <- list(coyote = ~1, + marten = ~1, + fisher = list(~1, coyote = ~1, marten = ~1)) + fit <- occuRNMulti(df, sf, umf[1:10,]) + expect_equal(names(coef(fit, type='lam')), + c("lam([coyote] (Intercept))", "lam([marten] (Intercept))", + "lam([fisher] (Intercept))", "lam([fisher:coyote] (Intercept))", + "lam([fisher:marten] (Intercept))" + )) + pr <- suppressMessages(predict(fit, 'state')) + expect_equal(length(pr), 3) + + # `sp1 --> sp2 --> sp3` + # plus sp1 --> sp3 + sf <- list(coyote = ~1, + marten = list(~1, coyote = ~x1), + fisher = list(~1, coyote=~1, marten = ~1)) + df <- list(coyote=~1, marten=~1, fisher=~1) + + # This model fit doesn't work well due to tiny sample size + fit <- suppressWarnings(occuRNMulti(df, sf, umf[1:10,])) + expect_equal(names(coef(fit, 'lam')), + c("lam([coyote] (Intercept))", "lam([marten] (Intercept))", + "lam([marten:coyote] (Intercept))", "lam([marten:coyote] x1)", + "lam([fisher] (Intercept))", "lam([fisher:coyote] (Intercept))", + "lam([fisher:marten] (Intercept))")) + pr <- suppressMessages(predict(fit, 'state', level=NULL)) + expect_equal(length(pr), 3) + +})