From bb2e82914e35efda07cf12258ff2c94ff279a3c2 Mon Sep 17 00:00:00 2001 From: Juniper Simonis Date: Sat, 20 Aug 2022 17:16:52 -0700 Subject: [PATCH 01/14] adding MNB mixture distributions base density and random functions for scalar and vector p multinomial negative binomial mixture models with roxygen documentation tests for functions vignette showing a more indepth worked example following the original paper by Haines 2020 that articulated the closed form solution of the model --- R/dNmixture_MNB.R | 210 +++++++++++++++++++ R/zzz.R | 35 +++- man/dNmixture_M.Rd | 108 ++++++++++ tests/testthat/test-Nmixture_M.R | 227 ++++++++++++++++++++ vignettes/Multinomial_N_mixtures.Rmd | 303 +++++++++++++++++++++++++++ 5 files changed, 882 insertions(+), 1 deletion(-) create mode 100644 R/dNmixture_MNB.R create mode 100644 man/dNmixture_M.Rd create mode 100644 tests/testthat/test-Nmixture_M.R create mode 100644 vignettes/Multinomial_N_mixtures.Rmd diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R new file mode 100644 index 0000000..1b41644 --- /dev/null +++ b/R/dNmixture_MNB.R @@ -0,0 +1,210 @@ +# +# Multinomial Negative Binomial mixture models with scalar and vector p +# + + +# dNmixture_M +#' @title Multinomial N-mixture distribution for use in \code{nimble} models +#' +#' @description \code{dNmixture_MNB_s} and \code{dNmixture_MNB_v} provide Multinomial Negative Binomial (MNB) mixture distributions of abundance ("N-mixture") for use in \code{nimble} models. +#' +#' @name dNmixture_M +#' @aliases dNmixture_MNB_s dNmixture_MNB_v rNmixture_MNB_s rNmixture_MNB_v +#' +#' @author Juniper Simonis +#' +#' @param x vector of integer counts from a series of sampling occasions. +#' @param mu expected value of the Negative Binomial distribution of true abundance. +#' @param r shape parameter defining overdispersion. As \code{r} approaches 0, the Negative Binomial converges to a Poisson. +#' @param p detection probability (scalar for \code{dNmixture_MNB_s}, vector for \code{dNmixture_MNB_v}). +#' @param J integer number of searches. +#' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability. +#' @param n number of random draws, each returning a vector of length \code{len}. Currently only \code{n = 1} is supported, but the argument exists for standardization of "\code{r}" functions. +#' +#' @details These nimbleFunctions provide distributions that can be used directly in R or in \code{nimble} hierarchical models (via \code{\link[nimble]{nimbleCode}} and \code{\link[nimble]{nimbleModel}}). \cr +#' The distributions are implemented in closed-form following Haines (2020), which avoids infinite sum-based calculations. +#' +#' +#' @return For \code{dNmixture_s} and \code{dNmixture_v}: the probability (or likelihood) or log probability of observation vector \code{x}. \cr +#' For \code{rNmixture_s} and \code{rNmixture_v}: a simulated detection history, \code{x}. +#' +#' @seealso For binomial N-mixture models, see \code{\link{dNmixture}}. \cr +#' For occupancy models dealing with detection/nondetection data, see \code{\link{dOcc}}. \cr +#' For dynamic occupancy, see \code{\link{dDynOcc}}. +#' +#' @references +#' L. M. Haines. 2020. Multinomial N-mixture models for removal sampling. Biometrics 76:540-548. DOI 10.1111/biom.13147 +#' +#' @examples +#' +#' # Set up variables and parameters +#' # J: number of visits +#' # mu: mean abundance +#' # r: scale parameter on abundance distribution +#' # p: search-specific detection probabilities +#' +#' J <- 10 +#' mu <- 5 +#' r <- 2 +#' p <- runif(J, 0.4, 0.7) +#' +#' mut <- log(mu) +#' rt <- log(r) +#' +#' # Generate data +#' +#' yv <- rNmixture_MNB_v(n = 1, mu = mu, p = p, r = r, J = J) +#' +#' # Write the model code +#' +#' nc <- nimbleCode({ +#' +#' mut ~ dnorm(0, 1/2) +#' mu <- exp(mut) +#' +#' rt ~ dnorm(0, 0.1) +#' r <- exp(rt) +#' +#' for (j in 1:J) { +#' p[j] ~ dunif(0, 1) +#' } +#' +#' x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) +#' +#' }) +#' +#' # Construct the model object +#' +#' nmix <- nimbleModel(nc, +#' constants = list(J = J), +#' data = list(x = yv), +#' inits = list(mut = mut, +#' rt = rt, +#' p = p)) +#' +#' # Calculate log probability of data from the model +#' +#' nmix$calculate() +#' +#' +# nimbleOptions(checkNimbleFunction = FALSE) +NULL + +#' @rdname dNmixture_M +#' +#' @export +#' +dNmixture_MNB_s <- nimbleFunction( + run = function(x = double(1), + mu = double(), + p = double(), + r = double(), + J = double(), + log = integer(0, default = 0)) { + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + + term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + term2 <- r * log(r) + x_tot * log(mu) + term3 <- x_tot * log(p) + x_miss * log(1 - p) + term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) + logProb <- term1 + term2 + term3 + term4 + + if (log) return(logProb) + else return(exp(logProb)) + returnType(double()) + +}) + +#' @rdname dNmixture_M +#' +#' @export +#' +rNmixture_MNB_s <- nimbleFunction( + run = function(n = integer(), + mu = double(), + p = double(), + r = double(), + J = double()) { + + + prob <- numeric(J + 1) + for (i in 1:(J)) { + prob[i] <- pow(1 - p, i - 1) * p + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + ans <- numeric(J + 1) + n <- rnbinom(n = 1, size = r, prob = 1/(1 + (1/r) * mu)) + if (n > 0) { + ans <- rmulti(n = 1, size = n, prob = prob) + } + + return(ans[1:J]) + returnType(double(1)) +}) + + +#' @rdname dNmixture_M +#' +#' @export +#' +dNmixture_MNB_v <- nimbleFunction( + run = function(x = double(1), + mu = double(), + p = double(1), + r = double(), + J = double(), + log = integer(0, default = 0)) { + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + pp <- c(0, p) + prob <- numeric(J) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + ptot <- sum(prob) + + term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + term2 <- r * log(r) + x_tot * log(mu) + term3 <- sum(x * log(prob)) + term4 <- -(x_tot + r) * log(r + mu * ptot) + logProb <- term1 + term2 + term3 + term4 + + if (log) return(logProb) + else return(exp(logProb)) + returnType(double()) + +}) + + +#' @rdname dNmixture_M +#' +#' @export +#' +rNmixture_MNB_v <- nimbleFunction( + run = function(n = integer(), + mu = double(), + p = double(1), + r = double(), + J = double()) { + + + pp <- c(0, p) + prob <- numeric(J + 1) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + ans <- numeric(J + 1) + n <- rnbinom(n = 1, size = r, prob = 1/(1 + (1/r) * mu)) + if (n > 0) { + ans <- rmulti(n = 1, size = n, prob = prob) + } + + return(ans[1:J]) + returnType(double(1)) +}) diff --git a/R/zzz.R b/R/zzz.R index 20c1d0d..ba0080f 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -2,7 +2,7 @@ .onAttach <- function(libname, pkgname) { packageStartupMessage("Loading nimbleEcology. Registering the following distributions:\n ", - "dOcc", ", dDynOcc", ", dCJS", ", dHMM", ", dDHMM", ", dNmixture.\n") + "dOcc", ", dDynOcc", ", dCJS", ", dHMM", ", dDHMM", ", dNmixture", ", dNmixture_M.\n") # Register the distributions explicitly for two reasons: # 1. Avoid message to user about automatic registrations upon first use in a nimbleModel @@ -521,4 +521,37 @@ ) + registerDistributions(list( + dNmixture_MNB_s = list( + BUGSdist = "dNmixture_MNB_s(mu, p, r, J)", + Rdist = "dNmixture_MNB_s(mu, p, r, J)", + discrete = TRUE, + types = c('value = double(1)', + 'mu = double()', + 'p = double()', + 'r = double()', + 'J = double()' + ), + mixedSizes = FALSE, + pqAvail = FALSE + )), verbose = FALSE + ) + + + registerDistributions(list( + dNmixture_MNB_v = list( + BUGSdist = "dNmixture_MNB_v(mu, p, r, J)", + Rdist = "dNmixture_MNB_v(mu, p, r, J)", + discrete = TRUE, + types = c('value = double(1)', + 'mu = double()', + 'p = double(1)', + 'r = double()', + 'J = double()' + ), + mixedSizes = FALSE, + pqAvail = FALSE + )), verbose = FALSE + ) + })} diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd new file mode 100644 index 0000000..44809cb --- /dev/null +++ b/man/dNmixture_M.Rd @@ -0,0 +1,108 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dNmixture_MNB.R +\name{dNmixture_M} +\alias{dNmixture_M} +\alias{dNmixture_MNB_s} +\alias{dNmixture_MNB_v} +\alias{rNmixture_MNB_s} +\alias{rNmixture_MNB_v} +\title{Multinomial N-mixture distribution for use in \code{nimble} models} +\usage{ +dNmixture_MNB_s(x, mu, p, r, J, log = 0) + +rNmixture_MNB_s(n, mu, p, r, J) + +dNmixture_MNB_v(x, mu, p, r, J, log = 0) + +rNmixture_MNB_v(n, mu, p, r, J) +} +\arguments{ +\item{x}{vector of integer counts from a series of sampling occasions.} + +\item{mu}{expected value of the Negative Binomial distribution of true abundance.} + +\item{p}{detection probability (scalar for \code{dNmixture_MNB_s}, vector for \code{dNmixture_MNB_v}).} + +\item{r}{shape parameter defining overdispersion. As \code{r} approaches 0, the Negative Binomial converges to a Poisson.} + +\item{J}{integer number of searches.} + +\item{log}{\code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability.} + +\item{n}{number of random draws, each returning a vector of length \code{len}. Currently only \code{n = 1} is supported, but the argument exists for standardization of "\code{r}" functions.} +} +\value{ +For \code{dNmixture_s} and \code{dNmixture_v}: the probability (or likelihood) or log probability of observation vector \code{x}. \cr + For \code{rNmixture_s} and \code{rNmixture_v}: a simulated detection history, \code{x}. +} +\description{ +\code{dNmixture_MNB_s} and \code{dNmixture_MNB_v} provide Multinomial Negative Binomial (MNB) mixture distributions of abundance ("N-mixture") for use in \code{nimble} models. +} +\details{ +These nimbleFunctions provide distributions that can be used directly in R or in \code{nimble} hierarchical models (via \code{\link[nimble]{nimbleCode}} and \code{\link[nimble]{nimbleModel}}). \cr + The distributions are implemented in closed-form following Haines (2020), which avoids infinite sum-based calculations. +} +\examples{ + +# Set up variables and parameters +# J: number of visits +# mu: mean abundance +# r: scale parameter on abundance distribution +# p: search-specific detection probabilities + + J <- 10 + mu <- 5 + r <- 2 + p <- runif(J, 0.4, 0.7) + + mut <- log(mu) + rt <- log(r) + +# Generate data + + yv <- rNmixture_MNB_v(n = 1, mu = mu, p = p, r = r, J = J) + +# Write the model code + + nc <- nimbleCode({ + + mut ~ dnorm(0, 1/2) + mu <- exp(mut) + + rt ~ dnorm(0, 0.1) + r <- exp(rt) + + for (j in 1:J) { + p[j] ~ dunif(0, 1) + } + + x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) + + }) + +# Construct the model object + + nmix <- nimbleModel(nc, + constants = list(J = J), + data = list(x = yv), + inits = list(mut = mut, + rt = rt, + p = p)) + +# Calculate log probability of data from the model + + nmix$calculate() + + +} +\references{ +L. M. Haines. 2020. Multinomial N-mixture models for removal sampling. Biometrics 76:540-548. DOI 10.1111/biom.13147 +} +\seealso{ +For binomial N-mixture models, see \code{\link{dNmixture}}. \cr + For occupancy models dealing with detection/nondetection data, see \code{\link{dOcc}}. \cr + For dynamic occupancy, see \code{\link{dDynOcc}}. +} +\author{ +Juniper Simonis +} diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R new file mode 100644 index 0000000..e91f500 --- /dev/null +++ b/tests/testthat/test-Nmixture_M.R @@ -0,0 +1,227 @@ +# Test the N-mixture_M distribution nimbleFunction. + +# ----------------------------------------------------------------------------- +# 0. Load +# Set the context for testthat +context("Testing dNmixture_M-related functions.") + +# ----------------------------------------------------------------------------- +#### 1. Test dNmixture_MNB_v #### +test_that("dNmixture_MNB_v works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 0) + mu <- 8 + r <- 0.5 + p <- c(0.5, 0.3, 0.5, 0.4, 0.1) + J <- 5 + probX <- dNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J) + + # Manually calculate the correct answer + + x_tot <- sum(x) + x_miss <- sum(x * nimSeq(0, J - 1)) + pp <- nimC(0, p) + prob <- nimNumeric(J) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + ptot <- sum(prob) + term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + term2 <- r * log(r) + x_tot * log(mu) + term3 <- sum(x * log(prob)) + term4 <- -(x_tot + r) * log(r + mu * ptot) + logProb <- term1 + term2 + term3 + term4 + correctProbX <- exp(logProb) + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + + lProbX <- dNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + lCorrectProbX <- log(correctProbX) + + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + + CdNmixture_MNB_v <- compileNimble(dNmixture_MNB_v) + CprobX <- CdNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J) + + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + expect_equal(ClProbX, lProbX) + + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixture_MNB_v(mu = mu, p = p[1:5], r = r, + J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + m$calculate() + MlProbX <- m$getLogProb("x") + expect_equal(MlProbX, lProbX) + + # Compiled model + cm <- compileNimble(m) + cm$calculate() + CMlProbX <- cm$getLogProb("x") + expect_equal(CMlProbX, lProbX) + + # Test imputing value for all NAs + xNA <- c(NA, NA, NA, NA, NA) + mNA <- nimbleModel(nc, data = list(x = xNA), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + + + mNAConf <- configureMCMC(mNA) + mNAConf$addMonitors('x') + mNA_MCMC <- buildMCMC(mNAConf) + cmNA <- compileNimble(mNA, mNA_MCMC) + + set.seed(0) + cmNA$mNA_MCMC$run(10) + + # Did the imputed values come back? + expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[2]"]))) + + # Test simulation code + nSim <- 10 + xSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + xSim[i,] <- rNmixture_MNB_v(1, mu = mu, p = p, r = r, J = J) + } + + CrNmixture_MNB_v <- compileNimble(rNmixture_MNB_v) + CxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + CxSim[i,] <- CrNmixture_MNB_v(1, mu = mu, p = p, r = r, J = J) + } + expect_identical(xSim, CxSim) + + simNodes <- m$getDependencies(c('prob', 'mu'), self = FALSE) + mxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + m$simulate(simNodes, includeData = TRUE) + mxSim[i,] <- m$x + } + expect_identical(mxSim, xSim) + + CmxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + cm$simulate(simNodes, includeData = TRUE) + CmxSim[i,] <- cm$x + } + expect_identical(CmxSim, mxSim) + }) + +# ----------------------------------------------------------------------------- +#### 2. Test dNmixture_MNB_s #### +test_that("dNmixture_MNB_s works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 2) + mu <- 8 + p <- 0.4 + J <- 5 + + probX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J) + + # Uncompiled log probability + lProbX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + + # Compilation and compiled calculations + CdNmixture_MNB_s <- compileNimble(dNmixture_MNB_s) + CprobX <- CdNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + expect_equal(ClProbX, lProbX) + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixture_MNB_s(mu = mu, p = p, r = r, J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + m$calculate() + MlProbX <- m$getLogProb("x") + expect_equal(MlProbX, lProbX) + + # Compiled model + cm <- compileNimble(m) + cm$calculate() + CMlProbX <- cm$getLogProb("x") + expect_equal(CMlProbX, lProbX) + + # Test imputing value for all NAs + xNA <- c(NA, NA, NA, NA, NA) + mNA <- nimbleModel(nc, data = list(x = xNA), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + + + mNAConf <- configureMCMC(mNA) + mNAConf$addMonitors('x') + mNA_MCMC <- buildMCMC(mNAConf) + cmNA <- compileNimble(mNA, mNA_MCMC) + + set.seed(0) + cmNA$mNA_MCMC$run(10) + + # Did the imputed values come back? + expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[2]"]))) + + # Test simulation code + nSim <- 10 + xSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + xSim[i,] <- rNmixture_MNB_s(1, mu = mu, p = p, r = r, J = J) + } + + CrNmixture_MNB_s <- compileNimble(rNmixture_MNB_s) + CxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + CxSim[i,] <- CrNmixture_MNB_s(1, mu = mu, p = p, r = r, J = J) + } + expect_identical(xSim, CxSim) + + simNodes <- m$getDependencies(c('prob', 'mu'), self = FALSE) + mxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + m$simulate(simNodes, includeData = TRUE) + mxSim[i,] <- m$x + } + expect_identical(mxSim, xSim) + + CmxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + cm$simulate(simNodes, includeData = TRUE) + CmxSim[i,] <- cm$x + } + expect_identical(CmxSim, mxSim) + }) diff --git a/vignettes/Multinomial_N_mixtures.Rmd b/vignettes/Multinomial_N_mixtures.Rmd new file mode 100644 index 0000000..1908912 --- /dev/null +++ b/vignettes/Multinomial_N_mixtures.Rmd @@ -0,0 +1,303 @@ +--- +title: "Multinomial-Negative Binomial N-Mixture Models" +author: "Juniper Simonis" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +description: > + A worked example of multinomial mixture models. +vignette: > + %\VignetteIndexEntry{Multinomial-Negative Binomial N-Mixture Models} + %\VignetteEngine{knitr::rmarkdown} + \usepackage[utf8]{inputenc} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" + ,eval = TRUE ## uncomment this to build quickly without running code. +) +``` + + +This vignette details the use of two Multinomial - Negative Binomial mixture models in nimble: + - Nmixture_MNB_s is a multinomial with scalar p + - Nmixture_MNB_v is a multinomial with vector p + +This example follows that implemented in the foundational paper by Haines (2020), with additional complexities. + + +The distributions are parameterized in terms of mean (mu), scale (r) and visit-specific detection probability (p), where there are J visits to a given site. + +For extension of the basic model (with a single site), we replicate the distribution across R total sites (i in 1 ... R), each with variable numbers of visits (J_i) and J_tot total visits across all sites. + +We also include covariates on site abundance (z) and detection probability (w), the latter of which include site- and search-specific values (wR and wJ, respectively). +These covariates influence the mu and p parameters via generalized linear models (b and g parameters, respectively; log and logit links, respectively). + +To accommodate uneven sampling effort across sites, we use nested indexing approaches, and facilitate efficiency by pre-determining the index locations (spot1, spot2 corresponding to the first and last) for each site's values within the search-level vector. + + +```{r, results='hide', messages=FALSE,warnings=FALSE} +library(nimble) +library(nimbleEcology) +``` + +```{r} + + set.seed(1312) + + # Parameters and constants + + mu <- 28.05 + p <- 0.61 + r <- 11.4 + rt <- log(r) + b0 <- 3.077 + b1 <- 0.101 + g0 <- 0.5 + g1 <- 0.125 + g2 <- 0.125 + + R <- 20 + J_i <- rep(3:5, length.out = R) + J_tot <- sum(J_i) + + spot1 <- integer(R) + spot2 <- integer(R) + for (i in 1:R) { + spot1[i] <- (sum(J_i[1:i]) - J_i[i] + 1) + spot2[i] <- (sum(J_i[1:i])) + } + + + # Covariates + + z <- rnorm(R, 0, 1) + wR_draw <- rnorm(R, 0, 1) + wR <- rep(wR_draw, J_i) + wJ_draw <- rnorm(max(J_i), 0, 1) + wJ <- wJ_draw[sequence(J_i)] + + + # Derived parameters + # mu at the site level + # p at the site level for the scalar model + # at the visit level for the vector model + + mu_i <- exp(b0 + b1 * z) + p_s_i <- expit(g0 + g1 * wR_draw) + p_v_j <- expit(g0 + g1 * wR + g2 * wJ) + + + # Simulate observations from scalar and vector p + + ys <- NULL + for (i in 1:R) { + ys <- c(ys, rNmixture_MNB_s(n = 1, mu = mu_i[i], p = p_s_i[i], r = r, J = J_i[i])) + } + + yv <- NULL + for (i in 1:R) { + spots_in <- (sum(J_i[1:i]) - J_i[i] + 1):sum(J_i[1:i]) + yv <- c(yv, rNmixture_MNB_v(n = 1, mu = mu_i[i], p = p_v_j[spots_in], r = r, J = J_i[i])) + } +``` + +Next, we construct the scalar and vector model code: + + +```{r} + + + nc_s <- nimbleCode({ + + rt ~ dnorm(0, 0.1) + + b0 ~ dnorm(0, 0.5) + b1 ~ dnorm(0, 0.1) + + g0 ~ dnorm(0, 0.5) + g1 ~ dnorm(0, 0.1) + + + r <- exp(rt) + + for (i in 1:R) { + + mut_i[i] <- b0 + b1 * z[i] + mu_i[i] <- exp(mut_i[i]) + pt_i[i] <- g0 + g1 * wR[i] + p_i[i] <- expit(pt_i[i]) + + x[spot1[i]:spot2[i]] ~ dNmixture_MNB_s(mu = mu_i[i], p = p_i[i], r = r, J = J_i[i]) + + } + + }) + + nc_v <- nimbleCode({ + + rt ~ dnorm(0, 0.1) + + b0 ~ dnorm(0, 0.5) + b1 ~ dnorm(0, 0.1) + + g0 ~ dnorm(0, 0.5) + g1 ~ dnorm(0, 0.1) + g2 ~ dnorm(0, 0.1) + + r <- exp(rt) + + pt_j[1:J_tot] <- g0 + g1 * wR[1:J_tot] + g2 * wJ[1:J_tot] + p_j[1:J_tot] <- expit(pt_j[1:J_tot]) + + for (i in 1:R) { + + mut_i[i] <- b0 + b1 * z[i] + mu_i[i] <- exp(mut_i[i]) + + x[spot1[i]:spot2[i]] ~ dNmixture_MNB_v(mu = mu_i[i], p = p_j[spot1[i]:spot2[i]], r = r, J = J_i[i]) + + } + + }) + +``` + +This allows us to combine scalar and vector data and models cross-wise into four full model objects: + +```{r, results='hide', messages=FALSE,warnings=FALSE} + + # Scalar data, scalar model + + nmix_ss <- nimbleModel(nc_s, + constants = list(J_i = J_i, R = R, spot1 = spot1, spot2 = spot2, J_tot = J_tot), + data = list(x = ys, z = z, wR = wR), + inits = list(b0 = b0, + b1 = b1, + g0 = g0, + g1 = g1, + rt = rt)) + + + # Scalar data, vector model + + nmix_sv <- nimbleModel(nc_v, + constants = list(J_i = J_i, R = R, spot1 = spot1, spot2 = spot2, J_tot = J_tot), + data = list(x = ys, z = z, wR = wR, wJ = wJ), + inits = list(b0 = b0, + b1 = b1, + g0 = g0, + g1 = g1, + g2 = g2, + rt = rt)) + + + # Vector data, scalar model + + nmix_vs <- nimbleModel(nc_s, + constants = list(J_i = J_i, R = R, spot1 = spot1, spot2 = spot2, J_tot = J_tot), + data = list(x = yv, z = z, wR = wR), + inits = list(b0 = b0, + b1 = b1, + g0 = g0, + g1 = g1, + rt = rt)) + + + # Vector data, vector model + + nmix_vv <- nimbleModel(nc_v, + constants = list(J_i = J_i, R = R, spot1 = spot1, spot2 = spot2, J_tot = J_tot), + data = list(x = yv, z = z, wR = wR, wJ = wJ), + inits = list(b0 = b0, + b1 = b1, + g0 = g0, + g1 = g1, + g2 = g2, + rt = rt)) +``` + +Which we can then use in a variety of ways! + +```{r, results='hide', messages=FALSE,warnings=FALSE} + + # Calculate based on provided parameters + + nmix_ss$calculate() + nmix_sv$calculate() + nmix_vs$calculate() + nmix_vv$calculate() + + + # Compile and sample the models 10,000 times + + cnmix_ss <- compileNimble(nmix_ss) + nmix_ssMCMC <- buildMCMC(nmix_ss) + cnmix_ssMCMC <- compileNimble(nmix_ssMCMC, project = cnmix_ss) + samples_ss <- runMCMC(cnmix_ssMCMC, niter = 10000) + + cnmix_sv <- compileNimble(nmix_sv) + nmix_svMCMC <- buildMCMC(nmix_sv) + cnmix_svMCMC <- compileNimble(nmix_svMCMC, project = cnmix_sv) + samples_sv <- runMCMC(cnmix_svMCMC, niter = 10000) + + cnmix_vs <- compileNimble(nmix_vs) + nmix_vsMCMC <- buildMCMC(nmix_vs) + cnmix_vsMCMC <- compileNimble(nmix_vsMCMC, project = cnmix_vs) + samples_vs <- runMCMC(cnmix_vsMCMC, niter = 10000) + + cnmix_vv <- compileNimble(nmix_vv) + nmix_vvMCMC <- buildMCMC(nmix_vv) + cnmix_vvMCMC <- compileNimble(nmix_vvMCMC, project = cnmix_vv) + samples_vv <- runMCMC(cnmix_vvMCMC, niter = 10000) + +``` + +And plot the last 1,000 + +### Scalar data, scalar model + +```{r, echo = FALSE, results='hide', messages=FALSE,warnings=FALSE, fig.dim =c(6, 3), fig.align="center"} + par(mar = c(4, 4, 1, 1)) + plot(samples_ss[9001:10000, "b0"], type = "l", ylab = "b0", las = 1, xlab = "sample") + plot(samples_ss[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") + plot(samples_ss[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") + plot(samples_ss[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_ss[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") +``` + +### Scalar data, vector model + +```{r, echo = FALSE, results='hide', messages=FALSE,warnings=FALSE, fig.dim =c(6, 3), fig.align="center"} + par(mar = c(4, 4, 1, 1)) + plot(samples_sv[9001:10000, "b0"], type = "l", ylab = "b0", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "g2"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") +``` + +### Vector data, scalar model + +```{r, echo = FALSE, results='hide', messages=FALSE,warnings=FALSE, fig.dim =c(6, 3), fig.align="center"} + par(mar = c(4, 4, 1, 1)) + plot(samples_vs[9001:10000, "b0"], type = "l", ylab = "b0", las = 1, xlab = "sample") + plot(samples_vs[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") + plot(samples_vs[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") + plot(samples_vs[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_vs[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") +``` + +### Vector data, vector model + +```{r, echo = FALSE, results='hide', messages=FALSE,warnings=FALSE, fig.dim =c(6, 3), fig.align="center"} + par(mar = c(4, 4, 1, 1)) + plot(samples_vv[9001:10000, "b0"], type = "l", ylab = "b0", las = 1, xlab = "sample") + plot(samples_vv[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") + plot(samples_vv[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") + plot(samples_vv[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_vv[9001:10000, "g2"], type = "l", ylab = "g1", las = 1, xlab = "sample") + plot(samples_vv[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") +``` \ No newline at end of file From f1d8e573233c3213b25a07757065960fc94afdec Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Fri, 26 Aug 2022 10:29:00 -0700 Subject: [PATCH 02/14] add dNmixture_MNB and exported fns to DESCRIPTION, NAMESPACE for pkg install --- DESCRIPTION | 3 ++- NAMESPACE | 5 +++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 536e76f..1f04479 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,8 +27,9 @@ Collate: dHMM.R dOcc.R dNmixture.R + dNmixture_MNB.R zzz.R -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.0 Suggests: rmarkdown, knitr, diff --git a/NAMESPACE b/NAMESPACE index 33a94ba..b12a989 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -30,6 +30,9 @@ export(dNmixture_BBP_s) export(dNmixture_BBNB_oneObs) export(dNmixture_BBNB_v) export(dNmixture_BBNB_s) +export(dNmixture_BBNB_v) +export(dNmixture_MNB_v) +export(dNmixture_MNB_s) export(dOcc_s) export(dOcc_v) export(rCJS_ss) @@ -63,6 +66,8 @@ export(rNmixture_BBP_s) export(rNmixture_BBNB_oneObs) export(rNmixture_BBNB_v) export(rNmixture_BBNB_s) +export(rNmixture_MNB_s) +export(rNmixture_MNB_v) export(nimNmixPois_logFac) export(rOcc_s) export(rOcc_v) From 5efca5d04e2dcd75d98010fe3b84f3e76584844d Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Fri, 26 Aug 2022 15:28:40 -0700 Subject: [PATCH 03/14] fix bug in Nmixture_MNB_s test --- tests/testthat/test-Nmixture_M.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index e91f500..808dae5 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -136,6 +136,7 @@ test_that("dNmixture_MNB_s works", # Uncompiled calculation x <- c(1, 0, 1, 3, 2) mu <- 8 + r <- 0.5 p <- 0.4 J <- 5 From c8085947b1e53cc797268a90206e53a0605da837 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Fri, 26 Aug 2022 16:36:57 -0700 Subject: [PATCH 04/14] one more minor fix --- tests/testthat/test-Nmixture_M.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index 808dae5..7e0cbd5 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -111,7 +111,8 @@ test_that("dNmixture_MNB_v works", } expect_identical(xSim, CxSim) - simNodes <- m$getDependencies(c('prob', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for(i in 1:nSim) { @@ -209,7 +210,7 @@ test_that("dNmixture_MNB_s works", } expect_identical(xSim, CxSim) - simNodes <- m$getDependencies(c('prob', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for(i in 1:nSim) { From 75f2c1d15f6a08b57e127a5deeeda2c405821c4c Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Tue, 6 Sep 2022 17:29:48 -0700 Subject: [PATCH 05/14] use trunc inf sum for testing --- tests/testthat/test-Nmixture_M.R | 87 +++++++++++++++++++++++++++----- 1 file changed, 74 insertions(+), 13 deletions(-) diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index 7e0cbd5..911bebb 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -17,22 +17,53 @@ test_that("dNmixture_MNB_v works", J <- 5 probX <- dNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J) - # Manually calculate the correct answer - x_tot <- sum(x) - x_miss <- sum(x * nimSeq(0, J - 1)) - pp <- nimC(0, p) - prob <- nimNumeric(J) + # Calaculate likelihood using truncated infinite sum approach + K <- 1000 + pp <- c(0, p) + prob <- numeric(J + 1) for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] + prob[j] <- prod(1 - pp[1:j]) * p[j] } - ptot <- sum(prob) - term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - term2 <- r * log(r) + x_tot * log(mu) - term3 <- sum(x * log(prob)) - term4 <- -(x_tot + r) * log(r + mu * ptot) - logProb <- term1 + term2 + term3 + term4 - correctProbX <- exp(logProb) + prob[J + 1] <- 1 - sum(prob[1:J]) + + + min_N <- sum(x) + max_N <- K + + + # Manually calculate the correct answer via trunc inf sum + N_vec <- min_N:max_N + l_vec <- numeric(length(N_vec)) + for (i in 1:length(N_vec)) { + # log prob of N + lp_N <- dnbinom(x = N_vec[i], size = r, + prob = 1/(1 + (1/r) * mu), log = TRUE) + + # Log prob of data + lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, + log = TRUE) + + l_vec[i] <- exp(lp_N + lp_x) + } + correctProbX <- sum(l_vec) +# +# # Manually calculate the correct answer +# +# x_tot <- sum(x) +# x_miss <- sum(x * nimSeq(0, J - 1)) +# pp <- nimC(0, p) +# prob <- nimNumeric(J) +# for (j in 1:J) { +# prob[j] <- prod(1 - pp[1:j]) * p[j] +# } +# ptot <- sum(prob) +# term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) +# term2 <- r * log(r) + x_tot * log(mu) +# term3 <- sum(x * log(prob)) +# term4 <- -(x_tot + r) * log(r + mu * ptot) +# logProb <- term1 + term2 + term3 + term4 +# correctProbX <- exp(logProb) expect_equal(probX, correctProbX) @@ -143,6 +174,36 @@ test_that("dNmixture_MNB_s works", probX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J) + + # Calaculate likelihood using truncated infinite sum approach + K <- 1000 + prob <- numeric(J + 1) + for (i in 1:(J)) { + prob[i] <- pow(1 - p, i - 1) * p + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + min_N <- sum(x) + max_N <- K + + # Manually calculate the correct answer via trunc inf sum + N_vec <- min_N:max_N + l_vec <- numeric(length(N_vec)) + for (i in 1:length(N_vec)) { + # log prob of N + lp_N <- dnbinom(x = N_vec[i], size = r, + prob = 1/(1 + (1/r) * mu), log = TRUE) + + # Log prob of data + lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, + log = TRUE) + + l_vec[i] <- exp(lp_N + lp_x) + } + correctProbX <- sum(l_vec) + + expect_equal(correctProbX, probX) + # Uncompiled log probability lProbX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) From 6f6e30291f0398966e64e1c3020aaf47919d38c7 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Wed, 7 Sep 2022 11:30:30 -0700 Subject: [PATCH 06/14] add multinomial Poisson N-mixture --- DESCRIPTION | 1 + NAMESPACE | 6 + R/dNmixture_MNB.R | 36 ++-- R/dNmixture_MP.R | 126 ++++++++++++++ man/dNmixture_M.Rd | 24 ++- tests/testthat/test-Nmixture_M.R | 273 +++++++++++++++++++++++++++++++ 6 files changed, 446 insertions(+), 20 deletions(-) create mode 100644 R/dNmixture_MP.R diff --git a/DESCRIPTION b/DESCRIPTION index 1f04479..50cb959 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,6 +28,7 @@ Collate: dOcc.R dNmixture.R dNmixture_MNB.R + dNmixture_MP.R zzz.R RoxygenNote: 7.2.0 Suggests: diff --git a/NAMESPACE b/NAMESPACE index b12a989..d517272 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -33,6 +33,9 @@ export(dNmixture_BBNB_s) export(dNmixture_BBNB_v) export(dNmixture_MNB_v) export(dNmixture_MNB_s) +export(dNmixture_MP_v) +export(dNmixture_MP_s) + export(dOcc_s) export(dOcc_v) export(rCJS_ss) @@ -68,6 +71,9 @@ export(rNmixture_BBNB_v) export(rNmixture_BBNB_s) export(rNmixture_MNB_s) export(rNmixture_MNB_v) +export(dNmixture_MP_v) +export(dNmixture_MP_s) + export(nimNmixPois_logFac) export(rOcc_s) export(rOcc_v) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index 1b41644..0618002 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -6,17 +6,18 @@ # dNmixture_M #' @title Multinomial N-mixture distribution for use in \code{nimble} models #' -#' @description \code{dNmixture_MNB_s} and \code{dNmixture_MNB_v} provide Multinomial Negative Binomial (MNB) mixture distributions of abundance ("N-mixture") for use in \code{nimble} models. +#' @description \code{dNmixture_MNB_s} and \code{dNmixture_MNB_v} provide Multinomial Negative Binomial (MNB) mixture distributions of abundance ("N-mixture") for use in \code{nimble} models. #' #' @name dNmixture_M -#' @aliases dNmixture_MNB_s dNmixture_MNB_v rNmixture_MNB_s rNmixture_MNB_v +#' @aliases dNmixture_MNB_s dNmixture_MNB_v rNmixture_MNB_s rNmixture_MNB_v +#' dNmixture_MP_s dNmixture_MP_v rNmixture_MP_s rNmixture_MP_v #' #' @author Juniper Simonis #' #' @param x vector of integer counts from a series of sampling occasions. -#' @param mu expected value of the Negative Binomial distribution of true abundance. -#' @param r shape parameter defining overdispersion. As \code{r} approaches 0, the Negative Binomial converges to a Poisson. -#' @param p detection probability (scalar for \code{dNmixture_MNB_s}, vector for \code{dNmixture_MNB_v}). +#' @param mu expected value of the (negative binomial or Poisson) distribution of true abundance. +#' @param r shape parameter defining overdispersion. As \code{r} approaches 0, the negative binomial converges to a Poisson. +#' @param p detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}). #' @param J integer number of searches. #' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability. #' @param n number of random draws, each returning a vector of length \code{len}. Currently only \code{n = 1} is supported, but the argument exists for standardization of "\code{r}" functions. @@ -68,9 +69,9 @@ #' for (j in 1:J) { #' p[j] ~ dunif(0, 1) #' } -#' +#' #' x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) -#' +#' #' }) #' #' # Construct the model object @@ -113,7 +114,7 @@ dNmixture_MNB_s <- nimbleFunction( if (log) return(logProb) else return(exp(logProb)) - returnType(double()) + returnType(double()) }) @@ -161,11 +162,18 @@ dNmixture_MNB_v <- nimbleFunction( x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) pp <- c(0, p) - prob <- numeric(J) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - ptot <- sum(prob) + prob <- numeric(J) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + ptot <- sum(prob) + + # from _s: + # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + # term2 <- r * log(r) + x_tot * log(mu) + # term3 <- x_tot * log(p) + x_miss * log(1 - p) + # term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) + # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) term2 <- r * log(r) + x_tot * log(mu) @@ -175,7 +183,7 @@ dNmixture_MNB_v <- nimbleFunction( if (log) return(logProb) else return(exp(logProb)) - returnType(double()) + returnType(double()) }) diff --git a/R/dNmixture_MP.R b/R/dNmixture_MP.R new file mode 100644 index 0000000..ad8b3a3 --- /dev/null +++ b/R/dNmixture_MP.R @@ -0,0 +1,126 @@ + +#' @rdname dNmixture_M +#' +#' @export +#' +dNmixture_MP_s <- nimbleFunction( + run = function(x = double(1), + mu = double(), + p = double(), + J = double(), + log = integer(0, default = 0)) { + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + + + # Mathematica version of MNB: + ### loglik = Total[LogGamma[yrow+r]] - R*LogGamma[r] - ylogfact + + ### R*r*Log[r] + ytot*Log[mu] + + ### ytot*Log[p] + ysumj*Log[1-p] - + ### (ytot+R*r)*Log[r + mu*ptot]; + + # Mathematica version of MP: + ### loglik = ytot*Log[mu] + ytot*Log[p]+ ysumj*Log[1-p] - R*mu*ptot-ylogfact; + + # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + # term2 <- r * log(r) + x_tot * log(mu) + # term3 <- x_tot * log(p) + x_miss * log(1 - p) + # term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) + logProb <- x_tot * log(mu) + x_tot * log(p) + x_miss * log(1 - p) - + mu * (1 - (1 - p) ^ J) - sum(lfactorial(x)) + + + if (log) return(logProb) + else return(exp(logProb)) + returnType(double()) + + }) + +#' @rdname dNmixture_M +#' +#' @export +#' +rNmixture_MP_s <- nimbleFunction( + run = function(n = integer(), + mu = double(), + p = double(), + J = double()) { + + + prob <- numeric(J + 1) + for (i in 1:(J)) { + prob[i] <- pow(1 - p, i - 1) * p + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + ans <- numeric(J + 1) + n <- rpois(n = 1, lambda = mu) + if (n > 0) { + ans <- rmulti(n = 1, size = n, prob = prob) + } + + return(ans[1:J]) + returnType(double(1)) + }) + + +#' @rdname dNmixture_M +#' +#' @export +#' +dNmixture_MP_v <- nimbleFunction( + run = function(x = double(1), + mu = double(), + p = double(1), + J = double(), + log = integer(0, default = 0)) { + + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + + pp <- c(0, p) + prob <- numeric(J) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + ptot <- sum(prob) + + logProb <- x_tot * log(mu) + sum(x * log(prob)) - + mu * ptot - sum(lfactorial(x)) + + if (log) return(logProb) + else return(exp(logProb)) + returnType(double()) + + }) + + +#' @rdname dNmixture_M +#' +#' @export +#' +rNmixture_MP_v <- nimbleFunction( + run = function(n = integer(), + mu = double(), + p = double(1), + J = double()) { + + + pp <- c(0, p) + prob <- numeric(J + 1) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + ans <- numeric(J + 1) + n <- rpois(n = 1, lambda = mu) + if (n > 0) { + ans <- rmulti(n = 1, size = n, prob = prob) + } + + return(ans[1:J]) + returnType(double(1)) + }) diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd index 44809cb..a32466c 100644 --- a/man/dNmixture_M.Rd +++ b/man/dNmixture_M.Rd @@ -1,11 +1,15 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/dNmixture_MNB.R +% Please edit documentation in R/dNmixture_MNB.R, R/dNmixture_MP.R \name{dNmixture_M} \alias{dNmixture_M} \alias{dNmixture_MNB_s} \alias{dNmixture_MNB_v} \alias{rNmixture_MNB_s} \alias{rNmixture_MNB_v} +\alias{dNmixture_MP_s} +\alias{dNmixture_MP_v} +\alias{rNmixture_MP_s} +\alias{rNmixture_MP_v} \title{Multinomial N-mixture distribution for use in \code{nimble} models} \usage{ dNmixture_MNB_s(x, mu, p, r, J, log = 0) @@ -15,15 +19,23 @@ rNmixture_MNB_s(n, mu, p, r, J) dNmixture_MNB_v(x, mu, p, r, J, log = 0) rNmixture_MNB_v(n, mu, p, r, J) + +dNmixture_MP_s(x, mu, p, J, log = 0) + +rNmixture_MP_s(n, mu, p, J) + +dNmixture_MP_v(x, mu, p, J, log = 0) + +rNmixture_MP_v(n, mu, p, J) } \arguments{ \item{x}{vector of integer counts from a series of sampling occasions.} -\item{mu}{expected value of the Negative Binomial distribution of true abundance.} +\item{mu}{expected value of the (negative binomial or Poisson) distribution of true abundance.} -\item{p}{detection probability (scalar for \code{dNmixture_MNB_s}, vector for \code{dNmixture_MNB_v}).} +\item{p}{detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}).} -\item{r}{shape parameter defining overdispersion. As \code{r} approaches 0, the Negative Binomial converges to a Poisson.} +\item{r}{shape parameter defining overdispersion. As \code{r} approaches 0, the negative binomial converges to a Poisson.} \item{J}{integer number of searches.} @@ -75,9 +87,9 @@ These nimbleFunctions provide distributions that can be used directly in R or in for (j in 1:J) { p[j] ~ dunif(0, 1) } - + x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) - + }) # Construct the model object diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index 911bebb..3287f9d 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -288,3 +288,276 @@ test_that("dNmixture_MNB_s works", } expect_identical(CmxSim, mxSim) }) + + +# ----------------------------------------------------------------------------- +#### 1. Test dNmixture_MP_v #### +test_that("dNmixture_MP_v works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 0) + mu <- 8 + p <- c(0.5, 0.3, 0.5, 0.4, 0.1) + J <- 5 + probX <- dNmixture_MP_v(x = x, mu = mu, p = p, J = J) + + + # Calaculate likelihood using truncated infinite sum approach + K <- 1000 + pp <- c(0, p) + prob <- numeric(J + 1) + for (j in 1:J) { + prob[j] <- prod(1 - pp[1:j]) * p[j] + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + + min_N <- sum(x) + max_N <- K + + + # Manually calculate the correct answer via trunc inf sum + N_vec <- min_N:max_N + l_vec <- numeric(length(N_vec)) + for (i in 1:length(N_vec)) { + # log prob of N + lp_N <- dpois(x = N_vec[i], lambda = mu, log = TRUE) + + # Log prob of data + lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, + log = TRUE) + + l_vec[i] <- exp(lp_N + lp_x) + } + correctProbX <- sum(l_vec) + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + + lProbX <- dNmixture_MP_v(x = x, mu = mu, p = p, J = J, log = TRUE) + lCorrectProbX <- log(correctProbX) + + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + + CdNmixture_MP_v <- compileNimble(dNmixture_MP_v) + CprobX <- CdNmixture_MP_v(x = x, mu = mu, p = p, J = J) + + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_MP_v(x = x, mu = mu, p = p, J = J, log = TRUE) + expect_equal(ClProbX, lProbX) + + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixture_MP_v(mu = mu, p = p[1:5], J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + m$calculate() + MlProbX <- m$getLogProb("x") + expect_equal(MlProbX, lProbX) + + # Compiled model + cm <- compileNimble(m) + cm$calculate() + CMlProbX <- cm$getLogProb("x") + expect_equal(CMlProbX, lProbX) + + # Test imputing value for all NAs + xNA <- c(NA, NA, NA, NA, NA) + mNA <- nimbleModel(nc, data = list(x = xNA), + inits = list(mu = mu, + p = p), + constants = list(J = J)) + + + mNAConf <- configureMCMC(mNA) + mNAConf$addMonitors('x') + mNA_MCMC <- buildMCMC(mNAConf) + cmNA <- compileNimble(mNA, mNA_MCMC) + + set.seed(0) + cmNA$mNA_MCMC$run(10) + + # Did the imputed values come back? + expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[2]"]))) + + # Test simulation code + nSim <- 10 + xSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + xSim[i,] <- rNmixture_MP_v(1, mu = mu, p = p, J = J) + } + + CrNmixture_MP_v <- compileNimble(rNmixture_MP_v) + CxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + CxSim[i,] <- CrNmixture_MP_v(1, mu = mu, p = p, J = J) + } + expect_equal(xSim, CxSim) + + simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + + mxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + m$simulate(simNodes, includeData = TRUE) + mxSim[i,] <- m$x + } + expect_equal(mxSim, CxSim) + + CmxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + cm$simulate(simNodes, includeData = TRUE) + CmxSim[i,] <- cm$x + } + expect_identical(CmxSim, mxSim) + }) + + + +test_that("dNmixture_MP_s works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 0) + mu <- 8 + p <- c(0.4) + J <- 5 + probX <- dNmixture_MP_s(x = x, mu = mu, p = p, J = J) + + + # Calaculate likelihood using truncated infinite sum approach + K <- 1000 + prob <- numeric(J + 1) + for (i in 1:(J)) { + prob[i] <- pow(1 - p, i - 1) * p + } + prob[J + 1] <- 1 - sum(prob[1:J]) + + + min_N <- sum(x) + max_N <- K + + + # Manually calculate the correct answer via trunc inf sum + N_vec <- min_N:max_N + l_vec <- numeric(length(N_vec)) + for (i in 1:length(N_vec)) { + # log prob of N + lp_N <- dpois(x = N_vec[i], lambda = mu, log = TRUE) + + # Log prob of data + lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, + log = TRUE) + + l_vec[i] <- exp(lp_N + lp_x) + } + correctProbX <- sum(l_vec) + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + + lProbX <- dNmixture_MP_s(x = x, mu = mu, p = p, J = J, log = TRUE) + lCorrectProbX <- log(correctProbX) + + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + + CdNmixture_MP_s <- compileNimble(dNmixture_MP_s) + CprobX <- CdNmixture_MP_s(x = x, mu = mu, p = p, J = J) + + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_MP_s(x = x, mu = mu, p = p, J = J, log = TRUE) + expect_equal(ClProbX, lProbX) + + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixture_MP_s(mu = mu, p = p, J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(mu = mu, + p = p, r = r), + constants = list(J = J)) + m$calculate() + MlProbX <- m$getLogProb("x") + expect_equal(MlProbX, lProbX) + + # Compiled model + cm <- compileNimble(m) + cm$calculate() + CMlProbX <- cm$getLogProb("x") + expect_equal(CMlProbX, lProbX) + + # Test imputing value for all NAs + xNA <- c(NA, NA, NA, NA, NA) + mNA <- nimbleModel(nc, data = list(x = xNA), + inits = list(mu = mu, + p = p), + constants = list(J = J)) + + + mNAConf <- configureMCMC(mNA) + mNAConf$addMonitors('x') + mNA_MCMC <- buildMCMC(mNAConf) + cmNA <- compileNimble(mNA, mNA_MCMC) + + set.seed(0) + cmNA$mNA_MCMC$run(10) + + # Did the imputed values come back? + expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[2]"]))) + + # Test simulation code + nSim <- 10 + xSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + xSim[i,] <- rNmixture_MP_s(1, mu = mu, p = p, J = J) + } + + CrNmixture_MP_s <- compileNimble(rNmixture_MP_s) + CxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for (i in 1:nSim) { + CxSim[i,] <- CrNmixture_MP_s(1, mu = mu, p = p, J = J) + } + expect_equal(xSim, CxSim) + + simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + + mxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + m$simulate(simNodes, includeData = TRUE) + mxSim[i,] <- m$x + } + expect_equal(mxSim, CxSim) + + CmxSim <- array(NA, dim = c(nSim, J)) + set.seed(1) + for(i in 1:nSim) { + cm$simulate(simNodes, includeData = TRUE) + CmxSim[i,] <- cm$x + } + expect_identical(CmxSim, mxSim) + }) + From ebfa2f3295429f2e7f31937562211acaa7864669 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Wed, 7 Sep 2022 15:26:02 -0700 Subject: [PATCH 07/14] fix bugs in nmix mp test code --- NAMESPACE | 4 ++-- tests/testthat/test-Nmixture_M.R | 8 +++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index d517272..2fb3dcc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -71,8 +71,8 @@ export(rNmixture_BBNB_v) export(rNmixture_BBNB_s) export(rNmixture_MNB_s) export(rNmixture_MNB_v) -export(dNmixture_MP_v) -export(dNmixture_MP_s) +export(rNmixture_MP_v) +export(rNmixture_MP_s) export(nimNmixPois_logFac) export(rOcc_s) diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index 3287f9d..78fc109 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -291,7 +291,7 @@ test_that("dNmixture_MNB_s works", # ----------------------------------------------------------------------------- -#### 1. Test dNmixture_MP_v #### +#### 3. Test dNmixture_MP_v #### test_that("dNmixture_MP_v works", { # Uncompiled calculation @@ -360,7 +360,7 @@ test_that("dNmixture_MP_v works", m <- nimbleModel(code = nc, data = list(x = x), inits = list(mu = mu, - p = p, r = r), + p = p), constants = list(J = J)) m$calculate() MlProbX <- m$getLogProb("x") @@ -427,6 +427,8 @@ test_that("dNmixture_MP_v works", }) +# ----------------------------------------------------------------------------- +#### 4. Test dNmixture_MP_s #### test_that("dNmixture_MP_s works", { @@ -495,7 +497,7 @@ test_that("dNmixture_MP_s works", m <- nimbleModel(code = nc, data = list(x = x), inits = list(mu = mu, - p = p, r = r), + p = p), constants = list(J = J)) m$calculate() MlProbX <- m$getLogProb("x") From cdac84bb0d3927ceb71a4de58d7b578392cd48ec Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Fri, 9 Sep 2022 11:25:45 -0700 Subject: [PATCH 08/14] rename args and reparameterize NB in terms of lambda, theta to match dNmixture_*NB --- DESCRIPTION | 5 +- R/dNmixture_MNB.R | 94 ++++++++++++----------- R/dNmixture_MP.R | 28 +++---- R/zzz.R | 50 +++++++++--- tests/testthat/test-Nmixture_M.R | 128 +++++++++++++++---------------- 5 files changed, 171 insertions(+), 134 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 50cb959..893d01e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,16 +1,17 @@ Package: nimbleEcology Type: Package Title: Distributions for Ecological Models in 'nimble' -Version: 0.4.1 +Version: 0.5.0 Maintainer: Benjamin R. Goldstein Authors@R: c(person("Benjamin R.", "Goldstein", role = c("aut", "cre"), email = "ben.goldstein@berkeley.edu"), person("Daniel", "Turek", role = "aut"), person("Lauren", "Ponisio", role = "aut"), + person("Juniper", "Simonis", role = "ctb"), person("Perry", "de Valpine", role = "aut")) Date: 2021-11-1 Description: Common ecological distributions for 'nimble' models in the form of nimbleFunction objects. - Includes Cormack-Jolly-Seber, occupancy, dynamic occupancy, hidden Markov, dynamic hidden Markov, and N-mixture models. + Includes Cormack-Jolly-Seber, occupancy, dynamic occupancy, hidden Markov, dynamic hidden Markov, N-mixture, and multinomial N-mixture models. (Jolly (1965) , Seber (1965) , Turek et al. (2016) ). License: GPL-3 Copyright: Copyright (c) 2019, Perry de Valpine, Ben Goldstein, Daniel Turek, Lauren Ponisio diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index 0618002..07d1020 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -12,11 +12,11 @@ #' @aliases dNmixture_MNB_s dNmixture_MNB_v rNmixture_MNB_s rNmixture_MNB_v #' dNmixture_MP_s dNmixture_MP_v rNmixture_MP_s rNmixture_MP_v #' -#' @author Juniper Simonis +#' @author Juniper Simonis, Ben Goldstein #' #' @param x vector of integer counts from a series of sampling occasions. -#' @param mu expected value of the (negative binomial or Poisson) distribution of true abundance. -#' @param r shape parameter defining overdispersion. As \code{r} approaches 0, the negative binomial converges to a Poisson. +#' @param lambda expected value of the (negative binomial or Poisson) distribution of true abundance. +#' @param theta shape parameter defining overdispersion. #' @param p detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}). #' @param J integer number of searches. #' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability. @@ -40,37 +40,37 @@ #' #' # Set up variables and parameters #' # J: number of visits -#' # mu: mean abundance -#' # r: scale parameter on abundance distribution +#' # lambda: mean abundance +#' # theta: scale parameter on abundance distribution #' # p: search-specific detection probabilities #' #' J <- 10 -#' mu <- 5 -#' r <- 2 +#' lambda <- 5 +#' theta <- 2 #' p <- runif(J, 0.4, 0.7) #' -#' mut <- log(mu) -#' rt <- log(r) +#' lambdat <- log(mu) +#' thetat <- log(theta) #' #' # Generate data #' -#' yv <- rNmixture_MNB_v(n = 1, mu = mu, p = p, r = r, J = J) +#' yv <- rNmixture_MNB_v(n = 1, lambda = lambda, p = p, theta = theta, J = J) #' #' # Write the model code #' #' nc <- nimbleCode({ #' -#' mut ~ dnorm(0, 1/2) -#' mu <- exp(mut) +#' lambdat ~ dnorm(0, 1/2) +#' lambda <- exp(lambdat) #' -#' rt ~ dnorm(0, 0.1) -#' r <- exp(rt) +#' thetat ~ dnorm(0, 0.1) +#' theta <- exp(thetat) #' #' for (j in 1:J) { #' p[j] ~ dunif(0, 1) #' } #' -#' x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) +#' x[1:J] ~ dNmixture_MNB_v(lambda = lambda, p = p[1:J], theta = theta, J = J) #' #' }) #' @@ -79,8 +79,8 @@ #' nmix <- nimbleModel(nc, #' constants = list(J = J), #' data = list(x = yv), -#' inits = list(mut = mut, -#' rt = rt, +#' inits = list(lambdat = lambdat, +#' thetat = thetat, #' p = p)) #' #' # Calculate log probability of data from the model @@ -96,20 +96,22 @@ NULL #' @export #' dNmixture_MNB_s <- nimbleFunction( - run = function(x = double(1), - mu = double(), - p = double(), - r = double(), - J = double(), - log = integer(0, default = 0)) { + run = function(x = double(1), + lambda = double(), + p = double(), + theta = double(), + J = double(), + log = integer(0, default = 0)) { + + r <- 1 / theta x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - term2 <- r * log(r) + x_tot * log(mu) + term2 <- r * log(r) + x_tot * log(lambda) term3 <- x_tot * log(p) + x_miss * log(1 - p) - term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) + term4 <- -(x_tot + r) * log(r + lambda * (1 - (1 - p) ^ J)) logProb <- term1 + term2 + term3 + term4 if (log) return(logProb) @@ -123,11 +125,11 @@ dNmixture_MNB_s <- nimbleFunction( #' @export #' rNmixture_MNB_s <- nimbleFunction( - run = function(n = integer(), - mu = double(), - p = double(), - r = double(), - J = double()) { + run = function(n = integer(), + lambda = double(), + p = double(), + theta = double(), + J = double()) { prob <- numeric(J + 1) @@ -137,7 +139,7 @@ rNmixture_MNB_s <- nimbleFunction( prob[J + 1] <- 1 - sum(prob[1:J]) ans <- numeric(J + 1) - n <- rnbinom(n = 1, size = r, prob = 1/(1 + (1/r) * mu)) + n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) if (n > 0) { ans <- rmulti(n = 1, size = n, prob = prob) } @@ -152,12 +154,14 @@ rNmixture_MNB_s <- nimbleFunction( #' @export #' dNmixture_MNB_v <- nimbleFunction( - run = function(x = double(1), - mu = double(), - p = double(1), - r = double(), - J = double(), - log = integer(0, default = 0)) { + run = function(x = double(1), + lambda = double(), + p = double(1), + theta = double(), + J = double(), + log = integer(0, default = 0)) { + + r <- 1/theta x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) @@ -176,9 +180,9 @@ dNmixture_MNB_v <- nimbleFunction( # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - term2 <- r * log(r) + x_tot * log(mu) + term2 <- r * log(r) + x_tot * log(lambda) term3 <- sum(x * log(prob)) - term4 <- -(x_tot + r) * log(r + mu * ptot) + term4 <- -(x_tot + r) * log(r + lambda * ptot) logProb <- term1 + term2 + term3 + term4 if (log) return(logProb) @@ -193,11 +197,11 @@ dNmixture_MNB_v <- nimbleFunction( #' @export #' rNmixture_MNB_v <- nimbleFunction( - run = function(n = integer(), - mu = double(), - p = double(1), - r = double(), - J = double()) { + run = function(n = integer(), + lambda = double(), + p = double(1), + theta = double(), + J = double()) { pp <- c(0, p) @@ -208,7 +212,7 @@ rNmixture_MNB_v <- nimbleFunction( prob[J + 1] <- 1 - sum(prob[1:J]) ans <- numeric(J + 1) - n <- rnbinom(n = 1, size = r, prob = 1/(1 + (1/r) * mu)) + n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) if (n > 0) { ans <- rmulti(n = 1, size = n, prob = prob) } diff --git a/R/dNmixture_MP.R b/R/dNmixture_MP.R index ad8b3a3..bbae0af 100644 --- a/R/dNmixture_MP.R +++ b/R/dNmixture_MP.R @@ -4,11 +4,11 @@ #' @export #' dNmixture_MP_s <- nimbleFunction( - run = function(x = double(1), - mu = double(), - p = double(), - J = double(), - log = integer(0, default = 0)) { + run = function(x = double(1), + lambda = double(), + p = double(), + J = double(), + log = integer(0, default = 0)) { x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) @@ -27,8 +27,8 @@ dNmixture_MP_s <- nimbleFunction( # term2 <- r * log(r) + x_tot * log(mu) # term3 <- x_tot * log(p) + x_miss * log(1 - p) # term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) - logProb <- x_tot * log(mu) + x_tot * log(p) + x_miss * log(1 - p) - - mu * (1 - (1 - p) ^ J) - sum(lfactorial(x)) + logProb <- x_tot * log(lambda) + x_tot * log(p) + x_miss * log(1 - p) - + lambda * (1 - (1 - p) ^ J) - sum(lfactorial(x)) if (log) return(logProb) @@ -43,7 +43,7 @@ dNmixture_MP_s <- nimbleFunction( #' rNmixture_MP_s <- nimbleFunction( run = function(n = integer(), - mu = double(), + lambda = double(), p = double(), J = double()) { @@ -55,7 +55,7 @@ rNmixture_MP_s <- nimbleFunction( prob[J + 1] <- 1 - sum(prob[1:J]) ans <- numeric(J + 1) - n <- rpois(n = 1, lambda = mu) + n <- rpois(n = 1, lambda = lambda) if (n > 0) { ans <- rmulti(n = 1, size = n, prob = prob) } @@ -71,7 +71,7 @@ rNmixture_MP_s <- nimbleFunction( #' dNmixture_MP_v <- nimbleFunction( run = function(x = double(1), - mu = double(), + lambda = double(), p = double(1), J = double(), log = integer(0, default = 0)) { @@ -87,8 +87,8 @@ dNmixture_MP_v <- nimbleFunction( } ptot <- sum(prob) - logProb <- x_tot * log(mu) + sum(x * log(prob)) - - mu * ptot - sum(lfactorial(x)) + logProb <- x_tot * log(lambda) + sum(x * log(prob)) - + lambda * ptot - sum(lfactorial(x)) if (log) return(logProb) else return(exp(logProb)) @@ -103,7 +103,7 @@ dNmixture_MP_v <- nimbleFunction( #' rNmixture_MP_v <- nimbleFunction( run = function(n = integer(), - mu = double(), + lambda = double(), p = double(1), J = double()) { @@ -116,7 +116,7 @@ rNmixture_MP_v <- nimbleFunction( prob[J + 1] <- 1 - sum(prob[1:J]) ans <- numeric(J + 1) - n <- rpois(n = 1, lambda = mu) + n <- rpois(n = 1, lambda = lambda) if (n > 0) { ans <- rmulti(n = 1, size = n, prob = prob) } diff --git a/R/zzz.R b/R/zzz.R index ba0080f..ed91164 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -523,13 +523,13 @@ registerDistributions(list( dNmixture_MNB_s = list( - BUGSdist = "dNmixture_MNB_s(mu, p, r, J)", - Rdist = "dNmixture_MNB_s(mu, p, r, J)", + BUGSdist = "dNmixture_MNB_s(lambda, p, theta, J)", + Rdist = "dNmixture_MNB_s(lambda, p, theta, J)", discrete = TRUE, types = c('value = double(1)', - 'mu = double()', + 'lambda = double()', 'p = double()', - 'r = double()', + 'theta = double()', 'J = double()' ), mixedSizes = FALSE, @@ -540,18 +540,50 @@ registerDistributions(list( dNmixture_MNB_v = list( - BUGSdist = "dNmixture_MNB_v(mu, p, r, J)", - Rdist = "dNmixture_MNB_v(mu, p, r, J)", + BUGSdist = "dNmixture_MNB_v(lambda, p, theta, J)", + Rdist = "dNmixture_MNB_v(lambda, p, theta, J)", discrete = TRUE, types = c('value = double(1)', - 'mu = double()', + 'lambda = double()', 'p = double(1)', - 'r = double()', + 'theta = double()', 'J = double()' ), mixedSizes = FALSE, pqAvail = FALSE )), verbose = FALSE ) - + + + registerDistributions(list( + dNmixture_MP_s = list( + BUGSdist = "dNmixture_MP_s(lambda, p, J)", + Rdist = "dNmixture_MP_s(lambda, p, J)", + discrete = TRUE, + types = c('value = double(1)', + 'lambda = double()', + 'p = double()', + 'J = double()' + ), + mixedSizes = FALSE, + pqAvail = FALSE + )), verbose = FALSE + ) + + + registerDistributions(list( + dNmixture_MP_v = list( + BUGSdist = "dNmixture_MP_v(lambda, p, J)", + Rdist = "dNmixture_MP_v(lambda, p, J)", + discrete = TRUE, + types = c('value = double(1)', + 'lambda = double()', + 'p = double(1)', + 'J = double()' + ), + mixedSizes = FALSE, + pqAvail = FALSE + )), verbose = FALSE + ) + })} diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index 78fc109..ea179c6 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -11,11 +11,11 @@ test_that("dNmixture_MNB_v works", { # Uncompiled calculation x <- c(1, 0, 1, 3, 0) - mu <- 8 - r <- 0.5 + lambda <- 8 + theta <- 0.5 p <- c(0.5, 0.3, 0.5, 0.4, 0.1) J <- 5 - probX <- dNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J) + probX <- dNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, J = J) # Calaculate likelihood using truncated infinite sum approach @@ -37,8 +37,8 @@ test_that("dNmixture_MNB_v works", l_vec <- numeric(length(N_vec)) for (i in 1:length(N_vec)) { # log prob of N - lp_N <- dnbinom(x = N_vec[i], size = r, - prob = 1/(1 + (1/r) * mu), log = TRUE) + lp_N <- dnbinom(x = N_vec[i], size = 1/theta, + prob = 1/(1 + theta * lambda), log = TRUE) # Log prob of data lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, @@ -69,7 +69,7 @@ test_that("dNmixture_MNB_v works", # Uncompiled log probability - lProbX <- dNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + lProbX <- dNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, J = J, log = TRUE) lCorrectProbX <- log(correctProbX) expect_equal(lProbX, lCorrectProbX) @@ -77,25 +77,25 @@ test_that("dNmixture_MNB_v works", # Compilation and compiled calculations CdNmixture_MNB_v <- compileNimble(dNmixture_MNB_v) - CprobX <- CdNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J) + CprobX <- CdNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, J = J) expect_equal(CprobX, probX) - ClProbX <- CdNmixture_MNB_v(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + ClProbX <- CdNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, J = J, log = TRUE) expect_equal(ClProbX, lProbX) # Use in Nimble model nc <- nimbleCode({ - x[1:5] ~ dNmixture_MNB_v(mu = mu, p = p[1:5], r = r, + x[1:5] ~ dNmixture_MNB_v(lambda = lambda, p = p[1:5], theta = theta, J = J) }) m <- nimbleModel(code = nc, data = list(x = x), - inits = list(mu = mu, - p = p, r = r), + inits = list(lambda = lambda, + p = p, theta = theta), constants = list(J = J)) m$calculate() MlProbX <- m$getLogProb("x") @@ -110,8 +110,8 @@ test_that("dNmixture_MNB_v works", # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), - inits = list(mu = mu, - p = p, r = r), + inits = list(lambda = lambda, + p = p, theta = theta), constants = list(J = J)) @@ -131,18 +131,18 @@ test_that("dNmixture_MNB_v works", xSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - xSim[i,] <- rNmixture_MNB_v(1, mu = mu, p = p, r = r, J = J) + xSim[i,] <- rNmixture_MNB_v(1, lambda = lambda, p = p, theta = theta, J = J) } CrNmixture_MNB_v <- compileNimble(rNmixture_MNB_v) CxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - CxSim[i,] <- CrNmixture_MNB_v(1, mu = mu, p = p, r = r, J = J) + CxSim[i,] <- CrNmixture_MNB_v(1, lambda = lambda, p = p, theta = theta, J = J) } - expect_identical(xSim, CxSim) + expect_equal(xSim, CxSim) - simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'lambda'), self = FALSE) mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) @@ -150,7 +150,7 @@ test_that("dNmixture_MNB_v works", m$simulate(simNodes, includeData = TRUE) mxSim[i,] <- m$x } - expect_identical(mxSim, xSim) + expect_equal(mxSim, xSim) CmxSim <- array(NA, dim = c(nSim, J)) set.seed(1) @@ -158,7 +158,7 @@ test_that("dNmixture_MNB_v works", cm$simulate(simNodes, includeData = TRUE) CmxSim[i,] <- cm$x } - expect_identical(CmxSim, mxSim) + expect_equal(CmxSim, mxSim) }) # ----------------------------------------------------------------------------- @@ -167,12 +167,12 @@ test_that("dNmixture_MNB_s works", { # Uncompiled calculation x <- c(1, 0, 1, 3, 2) - mu <- 8 - r <- 0.5 + lambda <- 8 + theta <- 0.5 p <- 0.4 J <- 5 - probX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J) + probX <- dNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J) # Calaculate likelihood using truncated infinite sum approach @@ -191,8 +191,8 @@ test_that("dNmixture_MNB_s works", l_vec <- numeric(length(N_vec)) for (i in 1:length(N_vec)) { # log prob of N - lp_N <- dnbinom(x = N_vec[i], size = r, - prob = 1/(1 + (1/r) * mu), log = TRUE) + lp_N <- dnbinom(x = N_vec[i], size = 1/theta, + prob = 1/(1 + theta * lambda), log = TRUE) # Log prob of data lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, @@ -205,26 +205,26 @@ test_that("dNmixture_MNB_s works", expect_equal(correctProbX, probX) # Uncompiled log probability - lProbX <- dNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + lProbX <- dNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J, log = TRUE) # Compilation and compiled calculations CdNmixture_MNB_s <- compileNimble(dNmixture_MNB_s) - CprobX <- CdNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J) + CprobX <- CdNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J) expect_equal(CprobX, probX) - ClProbX <- CdNmixture_MNB_s(x = x, mu = mu, p = p, r = r, J = J, log = TRUE) + ClProbX <- CdNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J, log = TRUE) expect_equal(ClProbX, lProbX) # Use in Nimble model nc <- nimbleCode({ - x[1:5] ~ dNmixture_MNB_s(mu = mu, p = p, r = r, J = J) + x[1:5] ~ dNmixture_MNB_s(lambda = lambda, p = p, theta = theta, J = J) }) m <- nimbleModel(code = nc, data = list(x = x), - inits = list(mu = mu, - p = p, r = r), + inits = list(lambda = lambda, + p = p, theta = theta), constants = list(J = J)) m$calculate() MlProbX <- m$getLogProb("x") @@ -239,8 +239,8 @@ test_that("dNmixture_MNB_s works", # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), - inits = list(mu = mu, - p = p, r = r), + inits = list(lambda = lambda, + p = p, theta = theta), constants = list(J = J)) @@ -260,25 +260,25 @@ test_that("dNmixture_MNB_s works", xSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - xSim[i,] <- rNmixture_MNB_s(1, mu = mu, p = p, r = r, J = J) + xSim[i,] <- rNmixture_MNB_s(1, lambda = lambda, p = p, theta = theta, J = J) } CrNmixture_MNB_s <- compileNimble(rNmixture_MNB_s) CxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - CxSim[i,] <- CrNmixture_MNB_s(1, mu = mu, p = p, r = r, J = J) + CxSim[i,] <- CrNmixture_MNB_s(1, lambda = lambda, p = p, theta = theta, J = J) } - expect_identical(xSim, CxSim) + expect_equal(xSim, CxSim) - simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'lambda'), self = FALSE) mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for(i in 1:nSim) { m$simulate(simNodes, includeData = TRUE) mxSim[i,] <- m$x } - expect_identical(mxSim, xSim) + expect_equal(mxSim, xSim) CmxSim <- array(NA, dim = c(nSim, J)) set.seed(1) @@ -286,7 +286,7 @@ test_that("dNmixture_MNB_s works", cm$simulate(simNodes, includeData = TRUE) CmxSim[i,] <- cm$x } - expect_identical(CmxSim, mxSim) + expect_equal(CmxSim, mxSim) }) @@ -296,10 +296,10 @@ test_that("dNmixture_MP_v works", { # Uncompiled calculation x <- c(1, 0, 1, 3, 0) - mu <- 8 + lambda <- 8 p <- c(0.5, 0.3, 0.5, 0.4, 0.1) J <- 5 - probX <- dNmixture_MP_v(x = x, mu = mu, p = p, J = J) + probX <- dNmixture_MP_v(x = x, lambda = lambda, p = p, J = J) # Calaculate likelihood using truncated infinite sum approach @@ -321,7 +321,7 @@ test_that("dNmixture_MP_v works", l_vec <- numeric(length(N_vec)) for (i in 1:length(N_vec)) { # log prob of N - lp_N <- dpois(x = N_vec[i], lambda = mu, log = TRUE) + lp_N <- dpois(x = N_vec[i], lambda = lambda, log = TRUE) # Log prob of data lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, @@ -335,7 +335,7 @@ test_that("dNmixture_MP_v works", # Uncompiled log probability - lProbX <- dNmixture_MP_v(x = x, mu = mu, p = p, J = J, log = TRUE) + lProbX <- dNmixture_MP_v(x = x, lambda = lambda, p = p, J = J, log = TRUE) lCorrectProbX <- log(correctProbX) expect_equal(lProbX, lCorrectProbX) @@ -343,23 +343,23 @@ test_that("dNmixture_MP_v works", # Compilation and compiled calculations CdNmixture_MP_v <- compileNimble(dNmixture_MP_v) - CprobX <- CdNmixture_MP_v(x = x, mu = mu, p = p, J = J) + CprobX <- CdNmixture_MP_v(x = x, lambda = lambda, p = p, J = J) expect_equal(CprobX, probX) - ClProbX <- CdNmixture_MP_v(x = x, mu = mu, p = p, J = J, log = TRUE) + ClProbX <- CdNmixture_MP_v(x = x, lambda = lambda, p = p, J = J, log = TRUE) expect_equal(ClProbX, lProbX) # Use in Nimble model nc <- nimbleCode({ - x[1:5] ~ dNmixture_MP_v(mu = mu, p = p[1:5], J = J) + x[1:5] ~ dNmixture_MP_v(lambda = lambda, p = p[1:5], J = J) }) m <- nimbleModel(code = nc, data = list(x = x), - inits = list(mu = mu, + inits = list(lambda = lambda, p = p), constants = list(J = J)) m$calculate() @@ -375,7 +375,7 @@ test_that("dNmixture_MP_v works", # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), - inits = list(mu = mu, + inits = list(lambda = lambda, p = p), constants = list(J = J)) @@ -396,18 +396,18 @@ test_that("dNmixture_MP_v works", xSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - xSim[i,] <- rNmixture_MP_v(1, mu = mu, p = p, J = J) + xSim[i,] <- rNmixture_MP_v(1, lambda = lambda, p = p, J = J) } CrNmixture_MP_v <- compileNimble(rNmixture_MP_v) CxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - CxSim[i,] <- CrNmixture_MP_v(1, mu = mu, p = p, J = J) + CxSim[i,] <- CrNmixture_MP_v(1, lambda = lambda, p = p, J = J) } expect_equal(xSim, CxSim) - simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'lambda'), self = FALSE) mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) @@ -423,7 +423,7 @@ test_that("dNmixture_MP_v works", cm$simulate(simNodes, includeData = TRUE) CmxSim[i,] <- cm$x } - expect_identical(CmxSim, mxSim) + expect_equal(CmxSim, mxSim) }) @@ -434,10 +434,10 @@ test_that("dNmixture_MP_s works", { # Uncompiled calculation x <- c(1, 0, 1, 3, 0) - mu <- 8 + lambda <- 8 p <- c(0.4) J <- 5 - probX <- dNmixture_MP_s(x = x, mu = mu, p = p, J = J) + probX <- dNmixture_MP_s(x = x, lambda = lambda, p = p, J = J) # Calaculate likelihood using truncated infinite sum approach @@ -458,7 +458,7 @@ test_that("dNmixture_MP_s works", l_vec <- numeric(length(N_vec)) for (i in 1:length(N_vec)) { # log prob of N - lp_N <- dpois(x = N_vec[i], lambda = mu, log = TRUE) + lp_N <- dpois(x = N_vec[i], lambda = lambda, log = TRUE) # Log prob of data lp_x <- dmultinom(c(x, N_vec[i] - sum(x)), size = N_vec[i], prob = prob, @@ -472,7 +472,7 @@ test_that("dNmixture_MP_s works", # Uncompiled log probability - lProbX <- dNmixture_MP_s(x = x, mu = mu, p = p, J = J, log = TRUE) + lProbX <- dNmixture_MP_s(x = x, lambda = lambda, p = p, J = J, log = TRUE) lCorrectProbX <- log(correctProbX) expect_equal(lProbX, lCorrectProbX) @@ -480,23 +480,23 @@ test_that("dNmixture_MP_s works", # Compilation and compiled calculations CdNmixture_MP_s <- compileNimble(dNmixture_MP_s) - CprobX <- CdNmixture_MP_s(x = x, mu = mu, p = p, J = J) + CprobX <- CdNmixture_MP_s(x = x, lambda = lambda, p = p, J = J) expect_equal(CprobX, probX) - ClProbX <- CdNmixture_MP_s(x = x, mu = mu, p = p, J = J, log = TRUE) + ClProbX <- CdNmixture_MP_s(x = x, lambda = lambda, p = p, J = J, log = TRUE) expect_equal(ClProbX, lProbX) # Use in Nimble model nc <- nimbleCode({ - x[1:5] ~ dNmixture_MP_s(mu = mu, p = p, J = J) + x[1:5] ~ dNmixture_MP_s(lambda = lambda, p = p, J = J) }) m <- nimbleModel(code = nc, data = list(x = x), - inits = list(mu = mu, + inits = list(lambda = lambda, p = p), constants = list(J = J)) m$calculate() @@ -512,7 +512,7 @@ test_that("dNmixture_MP_s works", # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), - inits = list(mu = mu, + inits = list(lambda = lambda, p = p), constants = list(J = J)) @@ -533,18 +533,18 @@ test_that("dNmixture_MP_s works", xSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - xSim[i,] <- rNmixture_MP_s(1, mu = mu, p = p, J = J) + xSim[i,] <- rNmixture_MP_s(1, lambda = lambda, p = p, J = J) } CrNmixture_MP_s <- compileNimble(rNmixture_MP_s) CxSim <- array(NA, dim = c(nSim, J)) set.seed(1) for (i in 1:nSim) { - CxSim[i,] <- CrNmixture_MP_s(1, mu = mu, p = p, J = J) + CxSim[i,] <- CrNmixture_MP_s(1, lambda = lambda, p = p, J = J) } expect_equal(xSim, CxSim) - simNodes <- m$getDependencies(c('p', 'mu'), self = FALSE) + simNodes <- m$getDependencies(c('p', 'lambda'), self = FALSE) mxSim <- array(NA, dim = c(nSim, J)) set.seed(1) @@ -560,6 +560,6 @@ test_that("dNmixture_MP_s works", cm$simulate(simNodes, includeData = TRUE) CmxSim[i,] <- cm$x } - expect_identical(CmxSim, mxSim) + expect_equal(CmxSim, mxSim) }) From f52acdf525b9bc5ae3733d3507b878ce60fc71e5 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Sat, 10 Sep 2022 14:44:53 -0700 Subject: [PATCH 09/14] reparameterize multinomial Nmix in terms of lambda, theta, user-defined p --- R/dNmixture_MNB.R | 52 ++++++++++------------------ R/dNmixture_MP.R | 49 +++++++------------------- tests/testthat/test-Nmixture_M.R | 59 ++++++++------------------------ 3 files changed, 45 insertions(+), 115 deletions(-) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index 07d1020..b4f07a7 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -103,16 +103,21 @@ dNmixture_MNB_s <- nimbleFunction( J = double(), log = integer(0, default = 0)) { - r <- 1 / theta + p_vec <- rep(p, J) - x_tot <- sum(x) - x_miss <- sum(x * seq(0, J - 1)) + r <- 1/theta + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + + ptot <- sum(p_vec) + + term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) + term2 <- r * log(r) + x_tot * log(lambda) + term3 <- sum(x * log(p_vec)) + term4 <- -(x_tot + r) * log(r + lambda * ptot) + logProb <- term1 + term2 + term3 + term4 - term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - term2 <- r * log(r) + x_tot * log(lambda) - term3 <- x_tot * log(p) + x_miss * log(1 - p) - term4 <- -(x_tot + r) * log(r + lambda * (1 - (1 - p) ^ J)) - logProb <- term1 + term2 + term3 + term4 if (log) return(logProb) else return(exp(logProb)) @@ -132,11 +137,7 @@ rNmixture_MNB_s <- nimbleFunction( J = double()) { - prob <- numeric(J + 1) - for (i in 1:(J)) { - prob[i] <- pow(1 - p, i - 1) * p - } - prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(rep(p, J), 1 - (p*J)) ans <- numeric(J + 1) n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) @@ -165,23 +166,12 @@ dNmixture_MNB_v <- nimbleFunction( x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) - pp <- c(0, p) - prob <- numeric(J) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - ptot <- sum(prob) - - # from _s: - # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - # term2 <- r * log(r) + x_tot * log(mu) - # term3 <- x_tot * log(p) + x_miss * log(1 - p) - # term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) - # + + ptot <- sum(p) term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) term2 <- r * log(r) + x_tot * log(lambda) - term3 <- sum(x * log(prob)) + term3 <- sum(x * log(p)) term4 <- -(x_tot + r) * log(r + lambda * ptot) logProb <- term1 + term2 + term3 + term4 @@ -203,13 +193,7 @@ rNmixture_MNB_v <- nimbleFunction( theta = double(), J = double()) { - - pp <- c(0, p) - prob <- numeric(J + 1) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(p, 1 - sum(p)) ans <- numeric(J + 1) n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) diff --git a/R/dNmixture_MP.R b/R/dNmixture_MP.R index bbae0af..eff1786 100644 --- a/R/dNmixture_MP.R +++ b/R/dNmixture_MP.R @@ -13,22 +13,10 @@ dNmixture_MP_s <- nimbleFunction( x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) + p_vec <- rep(p, J) - # Mathematica version of MNB: - ### loglik = Total[LogGamma[yrow+r]] - R*LogGamma[r] - ylogfact + - ### R*r*Log[r] + ytot*Log[mu] + - ### ytot*Log[p] + ysumj*Log[1-p] - - ### (ytot+R*r)*Log[r + mu*ptot]; - - # Mathematica version of MP: - ### loglik = ytot*Log[mu] + ytot*Log[p]+ ysumj*Log[1-p] - R*mu*ptot-ylogfact; - - # term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) - # term2 <- r * log(r) + x_tot * log(mu) - # term3 <- x_tot * log(p) + x_miss * log(1 - p) - # term4 <- -(x_tot + r) * log(r + mu * (1 - (1 - p) ^ J)) - logProb <- x_tot * log(lambda) + x_tot * log(p) + x_miss * log(1 - p) - - lambda * (1 - (1 - p) ^ J) - sum(lfactorial(x)) + logProb <- x_tot * log(lambda) + sum(x * log(p_vec)) - + lambda * sum(p_vec) - sum(lfactorial(x)) if (log) return(logProb) @@ -47,12 +35,7 @@ rNmixture_MP_s <- nimbleFunction( p = double(), J = double()) { - - prob <- numeric(J + 1) - for (i in 1:(J)) { - prob[i] <- pow(1 - p, i - 1) * p - } - prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(rep(p, J), 1 - (p * J)) ans <- numeric(J + 1) n <- rpois(n = 1, lambda = lambda) @@ -80,15 +63,8 @@ dNmixture_MP_v <- nimbleFunction( x_tot <- sum(x) x_miss <- sum(x * seq(0, J - 1)) - pp <- c(0, p) - prob <- numeric(J) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - ptot <- sum(prob) - - logProb <- x_tot * log(lambda) + sum(x * log(prob)) - - lambda * ptot - sum(lfactorial(x)) + logProb <- x_tot * log(lambda) + sum(x * log(p)) - + lambda * sum(p) - sum(lfactorial(x)) if (log) return(logProb) else return(exp(logProb)) @@ -108,12 +84,13 @@ rNmixture_MP_v <- nimbleFunction( J = double()) { - pp <- c(0, p) - prob <- numeric(J + 1) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - prob[J + 1] <- 1 - sum(prob[1:J]) + # pp <- c(0, p) + # prob <- numeric(J + 1) + # for (j in 1:J) { + # prob[j] <- prod(1 - pp[1:j]) * p[j] + # } + # prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(p, 1 - sum(p)) ans <- numeric(J + 1) n <- rpois(n = 1, lambda = lambda) diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R index ea179c6..33db71a 100644 --- a/tests/testthat/test-Nmixture_M.R +++ b/tests/testthat/test-Nmixture_M.R @@ -13,20 +13,14 @@ test_that("dNmixture_MNB_v works", x <- c(1, 0, 1, 3, 0) lambda <- 8 theta <- 0.5 - p <- c(0.5, 0.3, 0.5, 0.4, 0.1) + p <- c(0.2, 0.1, 0.05, 0.2, 0.22) J <- 5 probX <- dNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, J = J) # Calaculate likelihood using truncated infinite sum approach K <- 1000 - pp <- c(0, p) - prob <- numeric(J + 1) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - prob[J + 1] <- 1 - sum(prob[1:J]) - + prob <- c(p, 1 - sum(p)) min_N <- sum(x) max_N <- K @@ -47,23 +41,7 @@ test_that("dNmixture_MNB_v works", l_vec[i] <- exp(lp_N + lp_x) } correctProbX <- sum(l_vec) -# -# # Manually calculate the correct answer -# -# x_tot <- sum(x) -# x_miss <- sum(x * nimSeq(0, J - 1)) -# pp <- nimC(0, p) -# prob <- nimNumeric(J) -# for (j in 1:J) { -# prob[j] <- prod(1 - pp[1:j]) * p[j] -# } -# ptot <- sum(prob) -# term1 <- lgamma(r + x_tot) - lgamma(r) - sum(lfactorial(x)) -# term2 <- r * log(r) + x_tot * log(mu) -# term3 <- sum(x * log(prob)) -# term4 <- -(x_tot + r) * log(r + mu * ptot) -# logProb <- term1 + term2 + term3 + term4 -# correctProbX <- exp(logProb) + expect_equal(probX, correctProbX) @@ -169,7 +147,7 @@ test_that("dNmixture_MNB_s works", x <- c(1, 0, 1, 3, 2) lambda <- 8 theta <- 0.5 - p <- 0.4 + p <- 0.05 J <- 5 probX <- dNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J) @@ -177,11 +155,7 @@ test_that("dNmixture_MNB_s works", # Calaculate likelihood using truncated infinite sum approach K <- 1000 - prob <- numeric(J + 1) - for (i in 1:(J)) { - prob[i] <- pow(1 - p, i - 1) * p - } - prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(rep(p, J), 1 - (p*J)) min_N <- sum(x) max_N <- K @@ -297,20 +271,14 @@ test_that("dNmixture_MP_v works", # Uncompiled calculation x <- c(1, 0, 1, 3, 0) lambda <- 8 - p <- c(0.5, 0.3, 0.5, 0.4, 0.1) + p <- c(0.4, 0.2, 0.1, 0.1, 0.05) J <- 5 probX <- dNmixture_MP_v(x = x, lambda = lambda, p = p, J = J) # Calaculate likelihood using truncated infinite sum approach K <- 1000 - pp <- c(0, p) - prob <- numeric(J + 1) - for (j in 1:J) { - prob[j] <- prod(1 - pp[1:j]) * p[j] - } - prob[J + 1] <- 1 - sum(prob[1:J]) - + prob <- c(p, 1 - sum(p)) min_N <- sum(x) max_N <- K @@ -435,18 +403,19 @@ test_that("dNmixture_MP_s works", # Uncompiled calculation x <- c(1, 0, 1, 3, 0) lambda <- 8 - p <- c(0.4) + p <- c(0.1) J <- 5 probX <- dNmixture_MP_s(x = x, lambda = lambda, p = p, J = J) # Calaculate likelihood using truncated infinite sum approach K <- 1000 - prob <- numeric(J + 1) - for (i in 1:(J)) { - prob[i] <- pow(1 - p, i - 1) * p - } - prob[J + 1] <- 1 - sum(prob[1:J]) + # prob <- numeric(J + 1) + # for (i in 1:(J)) { + # prob[i] <- pow(1 - p, i - 1) * p + # } + # prob[J + 1] <- 1 - sum(prob[1:J]) + prob <- c(rep(p, J), 1- p*J) min_N <- sum(x) From 4b49708842173fb574ad84cd45721648cd9c33a9 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Sat, 10 Sep 2022 14:54:45 -0700 Subject: [PATCH 10/14] update doc --- R/dNmixture_MNB.R | 22 ++++++++++++++------ R/dNmixture_MP.R | 7 ------- man/dNmixture_M.Rd | 52 ++++++++++++++++++++++++---------------------- 3 files changed, 43 insertions(+), 38 deletions(-) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index b4f07a7..f53714a 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -15,17 +15,27 @@ #' @author Juniper Simonis, Ben Goldstein #' #' @param x vector of integer counts from a series of sampling occasions. -#' @param lambda expected value of the (negative binomial or Poisson) distribution of true abundance. -#' @param theta shape parameter defining overdispersion. -#' @param p detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}). +#' @param lambda expected value of the (negative binomial or Poisson) +#' distribution of true abundance. +#' @param theta overdispersion parameter required for negative binomial +#' (*NB) models. theta is parameterized such that variance of +#' the negative binomial variable x is \code{lambda^2 * theta + lambda} +#' @param p detection probability (scalar for \code{dNmixture_M\*_s}, +#' vector for \code{dNmixture_M\*_v}). If a vector, p[i] is the probability +#' that an individual is observed in observation category i and no other +#' categories. \code{p} must sum to less than or equal to 1, although this +#' is not checked by the function for computational efficiency. The +#' probability of an individual going unobserved is \code{1 - sum(p)}. #' @param J integer number of searches. -#' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability. -#' @param n number of random draws, each returning a vector of length \code{len}. Currently only \code{n = 1} is supported, but the argument exists for standardization of "\code{r}" functions. +#' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or +#' \code{0} to return probability. +#' @param n number of random draws, each returning a vector of length +#' \code{len}. Currently only \code{n = 1} is supported, but the argument +#' exists for standardization of "\code{r}" functions. #' #' @details These nimbleFunctions provide distributions that can be used directly in R or in \code{nimble} hierarchical models (via \code{\link[nimble]{nimbleCode}} and \code{\link[nimble]{nimbleModel}}). \cr #' The distributions are implemented in closed-form following Haines (2020), which avoids infinite sum-based calculations. #' -#' #' @return For \code{dNmixture_s} and \code{dNmixture_v}: the probability (or likelihood) or log probability of observation vector \code{x}. \cr #' For \code{rNmixture_s} and \code{rNmixture_v}: a simulated detection history, \code{x}. #' diff --git a/R/dNmixture_MP.R b/R/dNmixture_MP.R index eff1786..7d81de5 100644 --- a/R/dNmixture_MP.R +++ b/R/dNmixture_MP.R @@ -83,13 +83,6 @@ rNmixture_MP_v <- nimbleFunction( p = double(1), J = double()) { - - # pp <- c(0, p) - # prob <- numeric(J + 1) - # for (j in 1:J) { - # prob[j] <- prod(1 - pp[1:j]) * p[j] - # } - # prob[J + 1] <- 1 - sum(prob[1:J]) prob <- c(p, 1 - sum(p)) ans <- numeric(J + 1) diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd index a32466c..8617b84 100644 --- a/man/dNmixture_M.Rd +++ b/man/dNmixture_M.Rd @@ -12,30 +12,32 @@ \alias{rNmixture_MP_v} \title{Multinomial N-mixture distribution for use in \code{nimble} models} \usage{ -dNmixture_MNB_s(x, mu, p, r, J, log = 0) +dNmixture_MNB_s(x, lambda, p, theta, J, log = 0) -rNmixture_MNB_s(n, mu, p, r, J) +rNmixture_MNB_s(n, lambda, p, theta, J) -dNmixture_MNB_v(x, mu, p, r, J, log = 0) +dNmixture_MNB_v(x, lambda, p, theta, J, log = 0) -rNmixture_MNB_v(n, mu, p, r, J) +rNmixture_MNB_v(n, lambda, p, theta, J) -dNmixture_MP_s(x, mu, p, J, log = 0) +dNmixture_MP_s(x, lambda, p, J, log = 0) -rNmixture_MP_s(n, mu, p, J) +rNmixture_MP_s(n, lambda, p, J) -dNmixture_MP_v(x, mu, p, J, log = 0) +dNmixture_MP_v(x, lambda, p, J, log = 0) -rNmixture_MP_v(n, mu, p, J) +rNmixture_MP_v(n, lambda, p, J) } \arguments{ \item{x}{vector of integer counts from a series of sampling occasions.} -\item{mu}{expected value of the (negative binomial or Poisson) distribution of true abundance.} +\item{lambda}{expected value of the (negative binomial or Poisson) distribution of true abundance.} \item{p}{detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}).} -\item{r}{shape parameter defining overdispersion. As \code{r} approaches 0, the negative binomial converges to a Poisson.} +\item{theta}{overdispersion parameter required for negative binomial +(*NB) models. theta is parameterized such that variance of +the negative binomial variable x is \code{lambda^2 * theta + lambda}} \item{J}{integer number of searches.} @@ -58,37 +60,37 @@ These nimbleFunctions provide distributions that can be used directly in R or in # Set up variables and parameters # J: number of visits -# mu: mean abundance -# r: scale parameter on abundance distribution +# lambda: mean abundance +# theta: scale parameter on abundance distribution # p: search-specific detection probabilities J <- 10 - mu <- 5 - r <- 2 + lambda <- 5 + theta <- 2 p <- runif(J, 0.4, 0.7) - mut <- log(mu) - rt <- log(r) + lambdat <- log(mu) + thetat <- log(theta) # Generate data - yv <- rNmixture_MNB_v(n = 1, mu = mu, p = p, r = r, J = J) + yv <- rNmixture_MNB_v(n = 1, lambda = lambda, p = p, theta = theta, J = J) # Write the model code nc <- nimbleCode({ - mut ~ dnorm(0, 1/2) - mu <- exp(mut) + lambdat ~ dnorm(0, 1/2) + lambda <- exp(lambdat) - rt ~ dnorm(0, 0.1) - r <- exp(rt) + thetat ~ dnorm(0, 0.1) + theta <- exp(thetat) for (j in 1:J) { p[j] ~ dunif(0, 1) } - x[1:J] ~ dNmixture_MNB_v(mu = mu, p = p[1:J], r = r, J = J) + x[1:J] ~ dNmixture_MNB_v(lambda = lambda, p = p[1:J], theta = theta, J = J) }) @@ -97,8 +99,8 @@ These nimbleFunctions provide distributions that can be used directly in R or in nmix <- nimbleModel(nc, constants = list(J = J), data = list(x = yv), - inits = list(mut = mut, - rt = rt, + inits = list(lambdat = lambdat, + thetat = thetat, p = p)) # Calculate log probability of data from the model @@ -116,5 +118,5 @@ For binomial N-mixture models, see \code{\link{dNmixture}}. \cr For dynamic occupancy, see \code{\link{dDynOcc}}. } \author{ -Juniper Simonis +Juniper Simonis, Ben Goldstein } From 4e58136bc0e11c8646adb06282df9dbe36d794e8 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Sat, 10 Sep 2022 14:54:58 -0700 Subject: [PATCH 11/14] update doc --- man/dNmixture_M.Rd | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd index 8617b84..cd46c48 100644 --- a/man/dNmixture_M.Rd +++ b/man/dNmixture_M.Rd @@ -31,9 +31,15 @@ rNmixture_MP_v(n, lambda, p, J) \arguments{ \item{x}{vector of integer counts from a series of sampling occasions.} -\item{lambda}{expected value of the (negative binomial or Poisson) distribution of true abundance.} +\item{lambda}{expected value of the (negative binomial or Poisson) +distribution of true abundance.} -\item{p}{detection probability (scalar for \code{dNmixture_M\*_s}, vector for \code{dNmixture_M\*_v}).} +\item{p}{detection probability (scalar for \code{dNmixture_M\*_s}, +vector for \code{dNmixture_M\*_v}). If a vector, p[i] is the probability +that an individual is observed in observation category i and no other +categories. \code{p} must sum to less than or equal to 1, although this +is not checked by the function for computational efficiency. The +probability of an individual going unobserved is \code{1 - sum(p)}.} \item{theta}{overdispersion parameter required for negative binomial (*NB) models. theta is parameterized such that variance of @@ -41,9 +47,12 @@ the negative binomial variable x is \code{lambda^2 * theta + lambda}} \item{J}{integer number of searches.} -\item{log}{\code{TRUE} or \code{1} to return log probability. \code{FALSE} or \code{0} to return probability.} +\item{log}{\code{TRUE} or \code{1} to return log probability. \code{FALSE} or +\code{0} to return probability.} -\item{n}{number of random draws, each returning a vector of length \code{len}. Currently only \code{n = 1} is supported, but the argument exists for standardization of "\code{r}" functions.} +\item{n}{number of random draws, each returning a vector of length +\code{len}. Currently only \code{n = 1} is supported, but the argument +exists for standardization of "\code{r}" functions.} } \value{ For \code{dNmixture_s} and \code{dNmixture_v}: the probability (or likelihood) or log probability of observation vector \code{x}. \cr From 0b0edf72cffa883eeaea691c5d50a23841f01603 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Tue, 18 Oct 2022 16:34:40 -0700 Subject: [PATCH 12/14] update vignette, doc --- R/dNmixture_MNB.R | 7 +- man/dNmixture_M.Rd | 7 +- vignettes/Multinomial_N_mixtures.Rmd | 122 +++++++++++++++++---------- 3 files changed, 88 insertions(+), 48 deletions(-) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index f53714a..f35d2af 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -24,8 +24,11 @@ #' vector for \code{dNmixture_M\*_v}). If a vector, p[i] is the probability #' that an individual is observed in observation category i and no other #' categories. \code{p} must sum to less than or equal to 1, although this -#' is not checked by the function for computational efficiency. The -#' probability of an individual going unobserved is \code{1 - sum(p)}. +#' is not checked by the function for computational efficiency. +#' If scalar, p is the probability that an individual is detected in each +#' observation type, and \code{p*length(x)} must be less than or equal to 1. +#' The probability of an individual going unobserved is \code{1 - sum(p)} +#' in the vector case and \code{1 - p*length(x)} in the scalar case. #' @param J integer number of searches. #' @param log \code{TRUE} or \code{1} to return log probability. \code{FALSE} or #' \code{0} to return probability. diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd index cd46c48..1bda8cb 100644 --- a/man/dNmixture_M.Rd +++ b/man/dNmixture_M.Rd @@ -38,8 +38,11 @@ distribution of true abundance.} vector for \code{dNmixture_M\*_v}). If a vector, p[i] is the probability that an individual is observed in observation category i and no other categories. \code{p} must sum to less than or equal to 1, although this -is not checked by the function for computational efficiency. The -probability of an individual going unobserved is \code{1 - sum(p)}.} +is not checked by the function for computational efficiency. +If scalar, p is the probability that an individual is detected in each +observation type, and \code{p*length(x)} must be less than or equal to 1. +The probability of an individual going unobserved is \code{1 - sum(p)} +in the vector case and \code{1 - p*length(x)} in the scalar case.} \item{theta}{overdispersion parameter required for negative binomial (*NB) models. theta is parameterized such that variance of diff --git a/vignettes/Multinomial_N_mixtures.Rmd b/vignettes/Multinomial_N_mixtures.Rmd index 1908912..3081068 100644 --- a/vignettes/Multinomial_N_mixtures.Rmd +++ b/vignettes/Multinomial_N_mixtures.Rmd @@ -1,6 +1,6 @@ --- title: "Multinomial-Negative Binomial N-Mixture Models" -author: "Juniper Simonis" +author: "Juniper Simonis and Ben Goldstein" date: "`r Sys.Date()`" output: rmarkdown::html_vignette description: > @@ -20,21 +20,35 @@ knitr::opts_chunk$set( ``` -This vignette details the use of two Multinomial - Negative Binomial mixture models in nimble: - - Nmixture_MNB_s is a multinomial with scalar p - - Nmixture_MNB_v is a multinomial with vector p +nimbleEcology now includes + - Nmixture_MNB_s is a multinomial N-mixture with negative binomial-distributed abundance and scalar (sample-invariant) detection + - Nmixture_MNB_v is a multinomial N-mixture with negative binomial-distributed abundance and vector (sample-variant) detection + - Nmixture_MP_s is a multinomial N-mixture with Poisson-distributed abundance and scalar (sample-invariant) detection + - Nmixture_MP_v is a multinomial N-mixture with Poisson-distributed abundance and vector (sample-variant) detection -This example follows that implemented in the foundational paper by Haines (2020), with additional complexities. +Multinomial N-mixture models describe a case where a population of individuals is each detected in one or more pools. Like the classical Examples of contexts where multinomial N-mixture models occur include removal sampling, multi-observer sampling, or -The distributions are parameterized in terms of mean (mu), scale (r) and visit-specific detection probability (p), where there are J visits to a given site. +For more information on multinomial N-mixture models, we recommend ____ -For extension of the basic model (with a single site), we replicate the distribution across R total sites (i in 1 ... R), each with variable numbers of visits (J_i) and J_tot total visits across all sites. -We also include covariates on site abundance (z) and detection probability (w), the latter of which include site- and search-specific values (wR and wJ, respectively). +This example follows that implemented in the [foundational paper by Haines](https://pubmed.ncbi.nlm.nih.gov/31513284/) (2020), with additional complexities. + +The distributions are parameterized in terms of mean expected abundance (lambda), scale (theta) and visit-specific detection probability (p), where there are J visits to a given site. + +For extension of the basic model (with a single site), we model counts across R total sites (i in 1 ... R), each with variable numbers of visits ($J_i$) and $J_{tot}$ total visits across all sites. + +We include covariates on site abundance (z) and detection probability (w), the latter of which include site- and search-specific values (wR and wJ, respectively). These covariates influence the mu and p parameters via generalized linear models (b and g parameters, respectively; log and logit links, respectively). -To accommodate uneven sampling effort across sites, we use nested indexing approaches, and facilitate efficiency by pre-determining the index locations (spot1, spot2 corresponding to the first and last) for each site's values within the search-level vector. +To accommodate uneven sampling effort across sites, we use nested indexing approaches and pre-determine the index locations (spot1, spot2 corresponding to the first and last) for each site's values within the search-level vector. + +We consider two cases: in one, sampling efficiency is constant across replicates. +In the other, sampling efficiency varies across replicates. Because the detection +probability for each replicate varies in both cases (i.e. the probability of +capture during sample 1 is different from sample 2 which requires that the +individual was NOT captured during sample 1) we use `dNmixture_MNB_v` in both +cases. ```{r, results='hide', messages=FALSE,warnings=FALSE} @@ -48,13 +62,13 @@ library(nimbleEcology) # Parameters and constants - mu <- 28.05 - p <- 0.61 - r <- 11.4 - rt <- log(r) + lambda <- 28.05 + p <- 0.11 + theta <- 11.4 + thetat <- log(theta) b0 <- 3.077 b1 <- 0.101 - g0 <- 0.5 + g0 <- -1.5 g1 <- 0.125 g2 <- 0.125 @@ -80,11 +94,11 @@ library(nimbleEcology) # Derived parameters - # mu at the site level + # lambda at the site level # p at the site level for the scalar model # at the visit level for the vector model - mu_i <- exp(b0 + b1 * z) + lambda_i <- exp(b0 + b1 * z) p_s_i <- expit(g0 + g1 * wR_draw) p_v_j <- expit(g0 + g1 * wR + g2 * wJ) @@ -93,13 +107,15 @@ library(nimbleEcology) ys <- NULL for (i in 1:R) { - ys <- c(ys, rNmixture_MNB_s(n = 1, mu = mu_i[i], p = p_s_i[i], r = r, J = J_i[i])) + ys <- c(ys, rNmixture_MNB_s(n = 1, lambda = lambda_i[i], p = p_s_i[i], + theta = theta, J = J_i[i])) } yv <- NULL for (i in 1:R) { spots_in <- (sum(J_i[1:i]) - J_i[i] + 1):sum(J_i[1:i]) - yv <- c(yv, rNmixture_MNB_v(n = 1, mu = mu_i[i], p = p_v_j[spots_in], r = r, J = J_i[i])) + yv <- c(yv, rNmixture_MNB_v(n = 1, lambda = lambda_i[i], p = p_v_j[spots_in], + theta = theta, J = J_i[i])) } ``` @@ -111,7 +127,7 @@ Next, we construct the scalar and vector model code: nc_s <- nimbleCode({ - rt ~ dnorm(0, 0.1) + thetat ~ dnorm(0, 0.1) b0 ~ dnorm(0, 0.5) b1 ~ dnorm(0, 0.1) @@ -120,24 +136,31 @@ Next, we construct the scalar and vector model code: g1 ~ dnorm(0, 0.1) - r <- exp(rt) - + theta <- exp(thetat) + for (i in 1:R) { - mut_i[i] <- b0 + b1 * z[i] - mu_i[i] <- exp(mut_i[i]) - pt_i[i] <- g0 + g1 * wR[i] - p_i[i] <- expit(pt_i[i]) + lambdat_i[i] <- b0 + b1 * z[i] + lambda_i[i] <- exp(lambdat_i[i]) + logit(p_i[i]) <- g0 + g1 * wR[spot1[i]] + + p_conv[i, 1] <- p_i[i] + for (j in 2:J_i[i]) { + p_conv[i, j] <- p_i[i] * (1 - p_conv[i, j - 1]) + } - x[spot1[i]:spot2[i]] ~ dNmixture_MNB_s(mu = mu_i[i], p = p_i[i], r = r, J = J_i[i]) + + x[spot1[i]:spot2[i]] ~ dNmixture_MNB_v(lambda = lambda_i[i], + p = p_conv[i, 1:J_i[i]], + theta = theta, J = J_i[i]) } - + }) nc_v <- nimbleCode({ - rt ~ dnorm(0, 0.1) + thetat ~ dnorm(0, 0.1) b0 ~ dnorm(0, 0.5) b1 ~ dnorm(0, 0.1) @@ -146,17 +169,28 @@ Next, we construct the scalar and vector model code: g1 ~ dnorm(0, 0.1) g2 ~ dnorm(0, 0.1) - r <- exp(rt) + theta <- exp(thetat) + - pt_j[1:J_tot] <- g0 + g1 * wR[1:J_tot] + g2 * wJ[1:J_tot] - p_j[1:J_tot] <- expit(pt_j[1:J_tot]) - for (i in 1:R) { - mut_i[i] <- b0 + b1 * z[i] - mu_i[i] <- exp(mut_i[i]) + lambdat_i[i] <- b0 + b1 * z[i] + lambda_i[i] <- exp(lambdat_i[i]) - x[spot1[i]:spot2[i]] ~ dNmixture_MNB_v(mu = mu_i[i], p = p_j[spot1[i]:spot2[i]], r = r, J = J_i[i]) + logit(p_j[i, 1]) <- g0 + g1 * wR[spot1[i]] + g2 * wJ[spot1[i]] + p_conv[i, 1] <- p_j[i, 1] + + for (j in 2:J_i[i]) { + logit(p_j[i, j]) <- g0 + g1 * wR[spot1[i] + j-1] + g2 * wJ[spot1[i] + j-1] + + p_conv[i, j] <- p_j[i, j] * (1 - p_conv[i, j - 1]) + } + + + x[spot1[i]:spot2[i]] ~ dNmixture_MNB_v(lambda = lambda_i[i], + p = p_conv[i, 1:J_i[i]], + theta = theta, + J = J_i[i]) } @@ -177,7 +211,7 @@ This allows us to combine scalar and vector data and models cross-wise into four b1 = b1, g0 = g0, g1 = g1, - rt = rt)) + thetat = thetat)) # Scalar data, vector model @@ -190,7 +224,7 @@ This allows us to combine scalar and vector data and models cross-wise into four g0 = g0, g1 = g1, g2 = g2, - rt = rt)) + thetat = thetat)) # Vector data, scalar model @@ -202,7 +236,7 @@ This allows us to combine scalar and vector data and models cross-wise into four b1 = b1, g0 = g0, g1 = g1, - rt = rt)) + thetat = thetat)) # Vector data, vector model @@ -215,7 +249,7 @@ This allows us to combine scalar and vector data and models cross-wise into four g0 = g0, g1 = g1, g2 = g2, - rt = rt)) + thetat = thetat)) ``` Which we can then use in a variety of ways! @@ -264,7 +298,7 @@ And plot the last 1,000 plot(samples_ss[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") plot(samples_ss[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") plot(samples_ss[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") - plot(samples_ss[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") + plot(samples_ss[9001:10000, "thetat"], type = "l", ylab = "thetat", las = 1, xlab = "sample") ``` ### Scalar data, vector model @@ -276,7 +310,7 @@ And plot the last 1,000 plot(samples_sv[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") plot(samples_sv[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") plot(samples_sv[9001:10000, "g2"], type = "l", ylab = "g1", las = 1, xlab = "sample") - plot(samples_sv[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") + plot(samples_sv[9001:10000, "thetat"], type = "l", ylab = "thetat", las = 1, xlab = "sample") ``` ### Vector data, scalar model @@ -287,7 +321,7 @@ And plot the last 1,000 plot(samples_vs[9001:10000, "b1"], type = "l", ylab = "b1", las = 1, xlab = "sample") plot(samples_vs[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") plot(samples_vs[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") - plot(samples_vs[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") + plot(samples_vs[9001:10000, "thetat"], type = "l", ylab = "thetat", las = 1, xlab = "sample") ``` ### Vector data, vector model @@ -299,5 +333,5 @@ And plot the last 1,000 plot(samples_vv[9001:10000, "g0"], type = "l", ylab = "g0", las = 1, xlab = "sample") plot(samples_vv[9001:10000, "g1"], type = "l", ylab = "g1", las = 1, xlab = "sample") plot(samples_vv[9001:10000, "g2"], type = "l", ylab = "g1", las = 1, xlab = "sample") - plot(samples_vv[9001:10000, "rt"], type = "l", ylab = "rt", las = 1, xlab = "sample") -``` \ No newline at end of file + plot(samples_vv[9001:10000, "thetat"], type = "l", ylab = "thetat", las = 1, xlab = "sample") +``` From 0f92489f1ebc129bddb87fbf3284b52469ee8384 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Tue, 18 Oct 2022 16:46:48 -0700 Subject: [PATCH 13/14] fix typo in example --- R/dNmixture_MNB.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R index f35d2af..a2d65b4 100644 --- a/R/dNmixture_MNB.R +++ b/R/dNmixture_MNB.R @@ -62,7 +62,7 @@ #' theta <- 2 #' p <- runif(J, 0.4, 0.7) #' -#' lambdat <- log(mu) +#' lambdat <- log(lambda) #' thetat <- log(theta) #' #' # Generate data From 0d2c75720e5b5790136c47e6ec1bf19588f209d8 Mon Sep 17 00:00:00 2001 From: Ben Goldstein Date: Tue, 18 Oct 2022 16:47:44 -0700 Subject: [PATCH 14/14] roxygenize --- man/dNmixture_M.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/dNmixture_M.Rd b/man/dNmixture_M.Rd index 1bda8cb..877ef9b 100644 --- a/man/dNmixture_M.Rd +++ b/man/dNmixture_M.Rd @@ -81,7 +81,7 @@ These nimbleFunctions provide distributions that can be used directly in R or in theta <- 2 p <- runif(J, 0.4, 0.7) - lambdat <- log(mu) + lambdat <- log(lambda) thetat <- log(theta) # Generate data