diff --git a/DESCRIPTION b/DESCRIPTION index 57dba16c..4f109814 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -53,6 +53,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' 'IDS.R' 'plotEffects.R' 'mixedModelTools.R' diff --git a/NAMESPACE b/NAMESPACE index b1f0e70c..d150ab54 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -23,7 +23,7 @@ importFrom(Rcpp, evalCpp) 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, occuCOP, IDS) + multmixOpen, nmixTTD, gdistremoval, goccu, occuCOP, IDS, occuRNMulti) export(removalPiFun, doublePiFun) export(makeRemPiFun, makeCrPiFun, makeCrPiFunMb, makeCrPiFunMh) diff --git a/R/RcppExports.R b/R/RcppExports.R index 70b57ac2..bb79441a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -73,6 +73,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, 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) { .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..a7a38a56 --- /dev/null +++ b/R/occuRNMulti.R @@ -0,0 +1,689 @@ +setClass("unmarkedFrameOccuRNMulti", + contains="unmarkedFrameOccuMulti") + +setClass("unmarkedFitOccuRNMulti", + representation( + detformulas = "list", + stateformulas = "list", + modelOccupancy = "numeric"), + contains = "unmarkedFit") + +occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy, + K = rep(10, length(stateformulas)), + starts, method="BFGS", se = TRUE, threads=1, ...){ + + 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)){ + stopifnot(all(modelOccupancy %in% names(modocc))) + modocc[modelOccupancy] <- 1 + not_allowed_occ <- colSums(dep) + sp_mismatch <- modocc & not_allowed_occ + if(any(sp_mismatch)){ + stop(paste("Species", + paste(names(modocc)[sp_mismatch], collapse=", "), + "cannot use occupancy model."), call.=FALSE) + } + } + + 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, function(y){ + if(all(is.na(y))) return(0) + max(y, na.rm=TRUE) + }) + }) + + nP <- max(gd$det_ind) + + has_occ <- sum(modocc) > 0 + if(has_occ){ + 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) + 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) + + 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, modocc, + 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, ...) + covMat <- invertHessian(fm, nP, se) + + fmAIC <- 2 * fm$value + 2 * nP + #---------------------------------------------------------------------------- + + #Format output--------------------------------------------------------------- + ests <- fm$par + names(ests) <- gd$par_names + + det <- unmarkedEstimate(name = "Detection", short.name = "p", + estimates = ests[det_rng], + covMat = as.matrix(covMat[det_rng, det_rng]), + invlink = "logistic", + invlinkGrad = "logistic.grad") + + 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, + modelOccupancy = modocc, + 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)){ + mf <- model.frame(s, sc, na.action=stats::na.pass) + return(model.matrix(s, mf)) + } + 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) + sc_long <- sc[sc_long_ind,] + if(is.null(umf@obsCovs)){ + oc <- sc_long + } else { + oc <- cbind(sc_long, umf@obsCovs) + } + + vv <- lapply(detformulas, function(x){ + 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){ + 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_start, det_end) + + 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) + + # Missing values + 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, + miss = miss, site_miss = site_miss) +}) + +create_dep_matrix <- function(stateformulas){ + S <- length(stateformulas) + dep <- matrix(0, S, S) + rownames(dep) <- colnames(dep) <- names(stateformulas) + + # 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 <- 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 +} + +# 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, ...){ + + 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) + + 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)){ + newdata <- NULL + X <- lapply(object@stateformulas, function(x){ + if(is.list(x)){ + lapply(x, function(z){ + mf <- model.frame(z, object@data@siteCovs, na.action=stats::na.pass) + model.matrix(z, mf) + }) + } else { + mf <- model.frame(x, object@data@siteCovs, na.action=stats::na.pass) + model.matrix(x, mf) + } + }) + 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(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, 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, 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) + }, point_est, cis, SIMPLIFY=FALSE) + + if(!missing(species)){ + out <- out[[species]] + } + return(out) + } else if(type == "det"){ + + 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) + 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 + } 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 + }) + names(out) <- det_species + if(length(out) == 1) out <- out[[1]] + return(out) + } else { + stop("type must be state or det", call.=FALSE) + } +}) + +# 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, + newdata){ + 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){ + 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]] + for (j in 2:length(Xsub)){ + other_sp <- names(Xsub)[j] + + # 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) + 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]] + for (j in 2:length(Xsub)){ + other_sp <- names(Xsub)[j] + + # 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) + post[[i]][,k] <- ests[[i]] + } + } + post +} + + +setMethod("getP", "unmarkedFitOccuRNMulti", function(object) +{ + + ylist <- object@data@ylist + S <- length(ylist) + M <- nrow(ylist[[1]]) + J <- ncol(ylist[[1]]) + + 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(object@modelOccupancy[i]){ + # If occupancy, fitted is psi * p + state[[i]] * p[[i]] + } else { + # If abundance + 1 - exp(-state[[i]]*p[[i]]) ## analytical integration. + } + }) + 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) +}) + +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) +}) 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} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 455bc52c..345a9a2e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -470,6 +470,30 @@ 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, 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; + 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< 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); + 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, miss, site_miss, 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) { @@ -564,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}, @@ -586,6 +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, 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 new file mode 100644 index 00000000..f1029400 --- /dev/null +++ b/src/nll_occuRNMulti.cpp @@ -0,0 +1,263 @@ +#include +#include "utils.h" + +#ifdef _OPENMP + #include +#endif + +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, + arma::irowvec miss){ + double f = R::dpois(k, lam, 0); + double p; + double g = 0; + for (int j = 0; j(ylist[0]); + int M = y1.n_rows; + int 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 { + // Unused at the moment I think + 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; + 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), miss.row(m)); + } else { + 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; + + //----------------------------------------------------------------------- + } else if((S == 3) & (dep(1,0) & dep(2,1))){ // sp1 --> sp2 --> 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), 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), miss.row(m)); + } else { + 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; + } + + 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), 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), miss.row(m)); + } else { + 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; + } + + lik_m += fg1 * lik_sp2; + + //----------------------------------------------------------------------- + } else if((S == 3) & (dep(1,0) & dep(2,0))){ // sp2 <-- sp1 --> sp3 + + + 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), miss.row(m)); + } else { + 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), miss.row(m)); + } else { + 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; + } 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), 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), miss.row(m)); + } else { + 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; + } + + lik_m += fg1 * lik_sp2; + + } + } + + nll -= log(lik_m); + } + return nll; +} 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) + +})