diff --git a/DESCRIPTION b/DESCRIPTION index 536e76f..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 @@ -27,8 +28,10 @@ Collate: dHMM.R dOcc.R dNmixture.R + dNmixture_MNB.R + dNmixture_MP.R zzz.R -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.0 Suggests: rmarkdown, knitr, diff --git a/NAMESPACE b/NAMESPACE index 33a94ba..2fb3dcc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -30,6 +30,12 @@ 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(dNmixture_MP_v) +export(dNmixture_MP_s) + export(dOcc_s) export(dOcc_v) export(rCJS_ss) @@ -63,6 +69,11 @@ 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(rNmixture_MP_v) +export(rNmixture_MP_s) + export(nimNmixPois_logFac) export(rOcc_s) export(rOcc_v) diff --git a/R/dNmixture_MNB.R b/R/dNmixture_MNB.R new file mode 100644 index 0000000..a2d65b4 --- /dev/null +++ b/R/dNmixture_MNB.R @@ -0,0 +1,219 @@ +# +# 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 +#' dNmixture_MP_s dNmixture_MP_v rNmixture_MP_s rNmixture_MP_v +#' +#' @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 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. +#' 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. +#' @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 +#' # lambda: mean abundance +#' # theta: scale parameter on abundance distribution +#' # p: search-specific detection probabilities +#' +#' J <- 10 +#' lambda <- 5 +#' theta <- 2 +#' p <- runif(J, 0.4, 0.7) +#' +#' lambdat <- log(lambda) +#' thetat <- log(theta) +#' +#' # Generate data +#' +#' yv <- rNmixture_MNB_v(n = 1, lambda = lambda, p = p, theta = theta, J = J) +#' +#' # Write the model code +#' +#' nc <- nimbleCode({ +#' +#' lambdat ~ dnorm(0, 1/2) +#' lambda <- exp(lambdat) +#' +#' thetat ~ dnorm(0, 0.1) +#' theta <- exp(thetat) +#' +#' for (j in 1:J) { +#' p[j] ~ dunif(0, 1) +#' } +#' +#' x[1:J] ~ dNmixture_MNB_v(lambda = lambda, p = p[1:J], theta = theta, J = J) +#' +#' }) +#' +#' # Construct the model object +#' +#' nmix <- nimbleModel(nc, +#' constants = list(J = J), +#' data = list(x = yv), +#' inits = list(lambdat = lambdat, +#' thetat = thetat, +#' 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), + lambda = double(), + p = double(), + theta = double(), + J = double(), + log = integer(0, default = 0)) { + + p_vec <- rep(p, J) + + 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 + + + if (log) return(logProb) + else return(exp(logProb)) + returnType(double()) + +}) + +#' @rdname dNmixture_M +#' +#' @export +#' +rNmixture_MNB_s <- nimbleFunction( + run = function(n = integer(), + lambda = double(), + p = double(), + theta = double(), + J = double()) { + + + prob <- c(rep(p, J), 1 - (p*J)) + + ans <- numeric(J + 1) + n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) + 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), + 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)) + + 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(p)) + term4 <- -(x_tot + r) * log(r + lambda * 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(), + lambda = double(), + p = double(1), + theta = double(), + J = double()) { + + prob <- c(p, 1 - sum(p)) + + ans <- numeric(J + 1) + n <- rnbinom(n = 1, size = 1/theta, prob = 1/(1 + theta * lambda)) + if (n > 0) { + ans <- rmulti(n = 1, size = n, prob = prob) + } + + return(ans[1:J]) + returnType(double(1)) +}) diff --git a/R/dNmixture_MP.R b/R/dNmixture_MP.R new file mode 100644 index 0000000..7d81de5 --- /dev/null +++ b/R/dNmixture_MP.R @@ -0,0 +1,96 @@ + +#' @rdname dNmixture_M +#' +#' @export +#' +dNmixture_MP_s <- nimbleFunction( + 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)) + + p_vec <- rep(p, J) + + logProb <- x_tot * log(lambda) + sum(x * log(p_vec)) - + lambda * sum(p_vec) - 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(), + lambda = double(), + p = double(), + J = double()) { + + prob <- c(rep(p, J), 1 - (p * J)) + + ans <- numeric(J + 1) + n <- rpois(n = 1, lambda = lambda) + 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), + lambda = double(), + p = double(1), + J = double(), + log = integer(0, default = 0)) { + + + x_tot <- sum(x) + x_miss <- sum(x * seq(0, J - 1)) + + logProb <- x_tot * log(lambda) + sum(x * log(p)) - + lambda * sum(p) - 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(), + lambda = double(), + p = double(1), + J = double()) { + + prob <- c(p, 1 - sum(p)) + + ans <- numeric(J + 1) + n <- rpois(n = 1, lambda = lambda) + 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..ed91164 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,69 @@ ) + registerDistributions(list( + dNmixture_MNB_s = list( + BUGSdist = "dNmixture_MNB_s(lambda, p, theta, J)", + Rdist = "dNmixture_MNB_s(lambda, p, theta, J)", + discrete = TRUE, + types = c('value = double(1)', + 'lambda = double()', + 'p = double()', + 'theta = double()', + 'J = double()' + ), + mixedSizes = FALSE, + pqAvail = FALSE + )), verbose = FALSE + ) + + + registerDistributions(list( + dNmixture_MNB_v = list( + BUGSdist = "dNmixture_MNB_v(lambda, p, theta, J)", + Rdist = "dNmixture_MNB_v(lambda, p, theta, J)", + discrete = TRUE, + types = c('value = double(1)', + 'lambda = double()', + 'p = double(1)', + '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/man/dNmixture_M.Rd b/man/dNmixture_M.Rd new file mode 100644 index 0000000..877ef9b --- /dev/null +++ b/man/dNmixture_M.Rd @@ -0,0 +1,134 @@ +% Generated by roxygen2: do not edit by hand +% 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, lambda, p, theta, J, log = 0) + +rNmixture_MNB_s(n, lambda, p, theta, J) + +dNmixture_MNB_v(x, lambda, p, theta, J, log = 0) + +rNmixture_MNB_v(n, lambda, p, theta, J) + +dNmixture_MP_s(x, lambda, p, J, log = 0) + +rNmixture_MP_s(n, lambda, p, J) + +dNmixture_MP_v(x, lambda, p, J, log = 0) + +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{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. +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 +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{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 +# lambda: mean abundance +# theta: scale parameter on abundance distribution +# p: search-specific detection probabilities + + J <- 10 + lambda <- 5 + theta <- 2 + p <- runif(J, 0.4, 0.7) + + lambdat <- log(lambda) + thetat <- log(theta) + +# Generate data + + yv <- rNmixture_MNB_v(n = 1, lambda = lambda, p = p, theta = theta, J = J) + +# Write the model code + + nc <- nimbleCode({ + + lambdat ~ dnorm(0, 1/2) + lambda <- exp(lambdat) + + thetat ~ dnorm(0, 0.1) + theta <- exp(thetat) + + for (j in 1:J) { + p[j] ~ dunif(0, 1) + } + + x[1:J] ~ dNmixture_MNB_v(lambda = lambda, p = p[1:J], theta = theta, J = J) + + }) + +# Construct the model object + + nmix <- nimbleModel(nc, + constants = list(J = J), + data = list(x = yv), + inits = list(lambdat = lambdat, + thetat = thetat, + 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, Ben Goldstein +} diff --git a/tests/testthat/test-Nmixture_M.R b/tests/testthat/test-Nmixture_M.R new file mode 100644 index 0000000..33db71a --- /dev/null +++ b/tests/testthat/test-Nmixture_M.R @@ -0,0 +1,534 @@ +# 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) + lambda <- 8 + theta <- 0.5 + 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 + prob <- c(p, 1 - sum(p)) + + 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 = 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, + log = TRUE) + + l_vec[i] <- exp(lp_N + lp_x) + } + correctProbX <- sum(l_vec) + + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + + lProbX <- dNmixture_MNB_v(x = x, lambda = lambda, p = p, theta = theta, 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, lambda = lambda, p = p, theta = theta, J = J) + + expect_equal(CprobX, probX) + + 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(lambda = lambda, p = p[1:5], theta = theta, + J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + p = p, theta = theta), + 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(lambda = lambda, + p = p, theta = theta), + 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, 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, lambda = lambda, p = p, theta = theta, J = J) + } + expect_equal(xSim, CxSim) + + 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_equal(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_equal(CmxSim, mxSim) + }) + +# ----------------------------------------------------------------------------- +#### 2. Test dNmixture_MNB_s #### +test_that("dNmixture_MNB_s works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 2) + lambda <- 8 + theta <- 0.5 + p <- 0.05 + J <- 5 + + probX <- dNmixture_MNB_s(x = x, lambda = lambda, p = p, theta = theta, J = J) + + + # Calaculate likelihood using truncated infinite sum approach + K <- 1000 + prob <- c(rep(p, J), 1 - (p*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 = 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, + 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, 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, lambda = lambda, p = p, theta = theta, J = J) + expect_equal(CprobX, probX) + + 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(lambda = lambda, p = p, theta = theta, J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + p = p, theta = theta), + 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(lambda = lambda, + p = p, theta = theta), + 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, 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, lambda = lambda, p = p, theta = theta, J = J) + } + expect_equal(xSim, CxSim) + + 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_equal(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_equal(CmxSim, mxSim) + }) + + +# ----------------------------------------------------------------------------- +#### 3. Test dNmixture_MP_v #### +test_that("dNmixture_MP_v works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 0) + lambda <- 8 + 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 + prob <- c(p, 1 - sum(p)) + + 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 = lambda, 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, lambda = lambda, 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, lambda = lambda, p = p, J = J) + + expect_equal(CprobX, probX) + + 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(lambda = lambda, p = p[1:5], J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + p = p), + 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(lambda = lambda, + 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, 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, lambda = lambda, p = p, J = J) + } + expect_equal(xSim, CxSim) + + 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_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_equal(CmxSim, mxSim) + }) + + +# ----------------------------------------------------------------------------- +#### 4. Test dNmixture_MP_s #### + +test_that("dNmixture_MP_s works", + { + # Uncompiled calculation + x <- c(1, 0, 1, 3, 0) + lambda <- 8 + 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 <- c(rep(p, J), 1- p*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 = lambda, 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, lambda = lambda, 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, lambda = lambda, p = p, J = J) + + expect_equal(CprobX, probX) + + 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(lambda = lambda, p = p, J = J) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + p = p), + 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(lambda = lambda, + 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, 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, lambda = lambda, p = p, J = J) + } + expect_equal(xSim, CxSim) + + 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_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_equal(CmxSim, mxSim) + }) + diff --git a/vignettes/Multinomial_N_mixtures.Rmd b/vignettes/Multinomial_N_mixtures.Rmd new file mode 100644 index 0000000..3081068 --- /dev/null +++ b/vignettes/Multinomial_N_mixtures.Rmd @@ -0,0 +1,337 @@ +--- +title: "Multinomial-Negative Binomial N-Mixture Models" +author: "Juniper Simonis and Ben Goldstein" +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. +) +``` + + +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 + + +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 + +For more information on multinomial N-mixture models, we recommend ____ + + +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 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} +library(nimble) +library(nimbleEcology) +``` + +```{r} + + set.seed(1312) + + # Parameters and constants + + lambda <- 28.05 + p <- 0.11 + theta <- 11.4 + thetat <- log(theta) + b0 <- 3.077 + b1 <- 0.101 + g0 <- -1.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 + # lambda at the site level + # p at the site level for the scalar model + # at the visit level for the vector model + + lambda_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, 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, lambda = lambda_i[i], p = p_v_j[spots_in], + theta = theta, J = J_i[i])) + } +``` + +Next, we construct the scalar and vector model code: + + +```{r} + + + nc_s <- nimbleCode({ + + thetat ~ dnorm(0, 0.1) + + b0 ~ dnorm(0, 0.5) + b1 ~ dnorm(0, 0.1) + + g0 ~ dnorm(0, 0.5) + g1 ~ dnorm(0, 0.1) + + + theta <- exp(thetat) + + for (i in 1:R) { + + 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_v(lambda = lambda_i[i], + p = p_conv[i, 1:J_i[i]], + theta = theta, J = J_i[i]) + + } + + }) + + nc_v <- nimbleCode({ + + thetat ~ 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) + + theta <- exp(thetat) + + + for (i in 1:R) { + + lambdat_i[i] <- b0 + b1 * z[i] + lambda_i[i] <- exp(lambdat_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]) + + } + + }) + +``` + +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, + thetat = thetat)) + + + # 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, + thetat = thetat)) + + + # 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, + thetat = thetat)) + + + # 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, + thetat = thetat)) +``` + +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, "thetat"], type = "l", ylab = "thetat", 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, "thetat"], type = "l", ylab = "thetat", 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, "thetat"], type = "l", ylab = "thetat", 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, "thetat"], type = "l", ylab = "thetat", las = 1, xlab = "sample") +```