# growthmodels: Growth models In PVAClone: Population Viability Analysis with Data Cloning

## Description

Population growth model to be used in model fitting via `pva`.

## Usage

 ```1 2 3 4 5``` ```gompertz(obs.error = "none", fixed) ricker(obs.error = "none", fixed) thetalogistic(obs.error = "none", fixed) thetalogistic_D(obs.error = "none", fixed) bevertonholt(obs.error = "none", fixed) ```

## Arguments

 `obs.error` Character, describing the observation error. Can be `"none"`, `"poisson"`, or `"normal"`. `fixed` Named numeric vector or list with fixed parameter names and values. Can be used for providing alternative prior specifications, see Details and Examples.

## Details

These functions can be called in `pva` to fit the following growth models to a given population time series assuming both with and without observation error. When assuming the presence of observation error, either the Normal or the Poisson observation error model must be assumed within the state-space model formulation (Nadeem and Lele, 2012). The growth models are defined as follows.

Gompertz (`gompertz`):

x[t] = a + x[t-1] + bx[t-1] + E[t]

where x[t] is log abundance at time t and E[t] ~ Normal(0, sigma^2).

Ricker (`ricker`):

x[t] = x[t-1] + a + be^(x[t-1]) + E[t]

where x[t] is log abundance at time t and E[t] ~ Normal(0, sigma^2).

Theta-Logistic (`thetalogistic`):

x[t] = x[t-1] + r[1-(e^(x[t-1])/K)^theta] + E[t]

where x[t] is log abundance at time t and E[t] ~ Normal(0, sigma^2).

Theta-Logistic with Demographic Variability (`thetalogistic_D`):

x[t] = x[t-1] + r[1-(e^(x[t-1])/K)^theta] + E[t]

where x[t] is log abundance at time t and E[t] ~ Normal(0, sigma^2+sigma.d^2), where sigma.d^2 is the demographic variability. If sigma.d^2 is missing or fixed at zero, Theta-Logistic model is fitted instead.

Generilzed Beverton-Holt (`bevertonholt`):

x[t] = x[t-1] + r-[1+(e^(x[t-1])/K)^theta] + E[t]

where x[t] is log abundance at time t and E[t] ~ Normal(0, sigma^2).

Observation error models are described in the help page of `pva`.

The argument `fixed` can be used to fit the model assuming a priori values of a subset of the parameters. For instance, fixing theta equal to one reduces Theta-Logistic and Gnerelized Beverton-Holt models to Logistic and Beverton-Holt models respectively. The number of parameters that should be fixed at most is p-1, where p is the dimension of the full model. See examples below and in `pva` and `model.select`.

The `fixed` argument can be used to provide alternative prior specification using the BUGS language. In this case, values in `fixed` must be numeric. Use a list when real fixed values (numeric) and priors (character) are provided at the same time (see Examples). Alternative priors can be useful for testing insensitivity to priors, which is a diagnostic sign of data cloning convergence.

## Value

An S4 class of 'pvamodel' (see `pvamodel-class`)

## References

Nadeem, K., Lele S. R., 2012. Likelihood based population viability analysis in the presence of observation error. Oikos 121, 1656–1664.

`pvamodel-class`, `pva`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10``` ```gompertz() gompertz("poisson") ricker("normal") ricker("normal", fixed=c(a=5, sigma=0.5)) thetalogistic("none", fixed=c(theta=1)) bevertonholt("normal", fixed=c(theta=1)) ## alternative priors ricker("normal", fixed=c(a="a ~ dnorm(2, 1)"))@model ricker("normal", fixed=list(a="a ~ dnorm(2, 1)", sigma=0.5))@model ```

### Example output

```Loading required package: dcmle
dclone 2.2-0 	 2018-02-26
dcmle 0.3-1 	 2016-03-11

Attaching package: 'dcmle'

The following objects are masked from 'package:coda':

chanames, crosscorr.plot, gelman.diag, gelman.plot, geweke.diag,
heidel.diag, raftery.diag, varnames

PVAClone 0.1-6 	 2016-03-11
check out ?PVA for an overview
An object of class "pvamodel"
Slot "growth.model":
[1] "gompertz"

Slot "obs.error":
[1] "none"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
x[1,i] ~ dnorm(m, prc)
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + (1+b)*x[j-1,i]
}
}
c <- 1 + b
z <- 0.5 * log((1+c) / (1-c))
m <- a / (1-c)
prc <- (1-(c*c)) / (sigma*sigma)
a ~ dnorm(0, 0.01)
b ~ dunif(-0.999, -0.0001)
sigma <- exp(lnsigma)
lnsigma ~ dnorm(0, 1)
prcx <- 1/sigma^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 3

Slot "support":
Min Max
a              -Inf Inf
b     -2.000000e+00   0
sigma  2.220446e-16 Inf

Slot "params":
[1] "a"       "z"       "lnsigma"

Slot "varnames":
[1] "a"     "b"     "sigma"

Slot "fixed":
NULL

Slot "fancy":
[1] "Gompertz" NA

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("b" %in% vn)
mcmc[[i]][, "b"] <- atanh(mcmc[[i]][, "b"] + 1)
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("b" %in% vn)
vn[vn == "b"] <- "z"
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x69818c8>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("z" %in% vn)
mcmc[[i]][, "z"] <- tanh(mcmc[[i]][, "z"]) - 1
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("z" %in% vn)
vn[vn == "z"] <- "b"
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x69818c8>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(logx)
m <- which(is.na(data))
if (length(m) > 0) {
if (alt_obserror) {
stop("not yet implemented")
}
else {
y <- data
y[m] <- logx[m]
rval <- sum(dnorm(y[-1], mean = mle["a"] + mle["b"] *
y[-T], sd = mle["sigma"], log = TRUE))
}
}
else {
logd1 <- dnorm(data[-1], mean = mle["a"] + mle["b"] *
data[-T], sd = mle["sigma"], log = TRUE)
logd2 <- if (alt_obserror) {
dnorm(logx[-1], mean = mle["a"] + mle["b"] * logx[-T],
sd = mle["sigma"], log = TRUE)
}
else {
0
}
rval <- sum(logd1) + sum(logd2)
}
rval
}
<environment: 0x69818c8>

Slot "neffective":
function (obs)
sum(!is.na(obs))
<environment: 0x69818c8>

An object of class "pvamodel"
Slot "growth.model":
[1] "gompertz"

Slot "obs.error":
[1] "poisson"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
x[1,i] ~ dnorm(m, prc)
N[1,i] <- min(exp(x[1,i]), 10000)
O[1,i] ~ dpois(N[1,i])
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + (1+b)*x[j-1,i]
N[j,i] <- min(exp(x[j,i]), 10000)
O[j,i] ~ dpois(N[j,i])
}
}
c <- 1 + b
z <- 0.5 * log((1+c) / (1-c))
m <- a / (1-c)
prc <- (1-(c*c)) / (sigma*sigma)
a ~ dnorm(0, 0.01)
b ~ dunif(-0.999, -0.0001)
sigma <- exp(lnsigma)
lnsigma ~ dnorm(0, 1)
prcx <- 1/sigma^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 3

Slot "support":
Min Max
a              -Inf Inf
b     -2.000000e+00   0
sigma  2.220446e-16 Inf

Slot "params":
[1] "a"       "z"       "lnsigma"

Slot "varnames":
[1] "a"     "b"     "sigma"

Slot "fixed":
NULL

Slot "fancy":
[1] "Gompertz" "Poisson"

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("b" %in% vn)
mcmc[[i]][, "b"] <- atanh(mcmc[[i]][, "b"] + 1)
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("b" %in% vn)
vn[vn == "b"] <- "z"
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<bytecode: 0x87b1fe8>
<environment: 0x689c9c8>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("z" %in% vn)
mcmc[[i]][, "z"] <- tanh(mcmc[[i]][, "z"]) - 1
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("z" %in% vn)
vn[vn == "z"] <- "b"
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<bytecode: 0x8179cf0>
<environment: 0x689c9c8>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1], mean = mle["a"] + mle["b"] *
logx[-T], sd = mle["sigma"], log = TRUE)
logd2 <- dpois(data[-1], exp(logx[-1]), log = TRUE)
rval <- sum(logd1) + sum(logd2, na.rm = TRUE)
}
else {
stop("not yet implemented")
}
rval
}
<bytecode: 0x6886bb8>
<environment: 0x689c9c8>

Slot "neffective":
function (obs)
sum(!is.na(obs))
<bytecode: 0x2bb5f68>
<environment: 0x689c9c8>

An object of class "pvamodel"
Slot "growth.model":
[1] "ricker"

Slot "obs.error":
[1] "normal"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
N[1,i] <- exp(y[1,i])
x[1,i] <- y[1,i]
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]
N[j,i] <- min(exp(x[j,i]), 10000)
y[j,i] ~ dnorm(x[j,i], prcy)
}
}
sigma <- exp(lnsigma)
tau <- exp(lntau)
lnsigma ~ dnorm(0, 1)
lntau ~ dnorm(0, 1)
a ~ dnorm(0, 0.01)
b ~ dnorm(0, 10)
prcx <- 1/sigma^2
prcy <- 1/tau^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 4

Slot "support":
Min Max
a             -Inf Inf
b             -Inf Inf
sigma 2.220446e-16 Inf
tau   2.220446e-16 Inf

Slot "params":
[1] "a"       "b"       "lnsigma" "lntau"

Slot "varnames":
[1] "a"     "b"     "sigma" "tau"

Slot "fixed":
NULL

Slot "fancy":
[1] "Ricker" "Normal"

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7f14ab8>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7f14ab8>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1], mean = logx[-T] + mle["a"] +
mle["b"] * exp(logx[-T]), sd = mle["sigma"], log = TRUE)
logd2 <- dnorm(data[-1], mean = logx[-1], sd = mle["tau"],
log = TRUE)
rval <- sum(logd1) + sum(logd2, na.rm = TRUE)
}
else {
stop("not yet implemented")
}
rval
}
<environment: 0x7f14ab8>

Slot "neffective":
function (obs)
sum(!is.na(obs)) - 1
<environment: 0x7f14ab8>

An object of class "pvamodel"
Slot "growth.model":
[1] "ricker"

Slot "obs.error":
[1] "normal"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
N[1,i] <- exp(y[1,i])
x[1,i] <- y[1,i]
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]
N[j,i] <- min(exp(x[j,i]), 10000)
y[j,i] ~ dnorm(x[j,i], prcy)
}
}
sigma <- 0.5
tau <- exp(lntau)
lnsigma <- log(sigma)
lntau ~ dnorm(0, 1)
a <- 5
b ~ dnorm(0, 10)
prcx <- 1/sigma^2
prcy <- 1/tau^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 4

Slot "support":
Min Max
a             -Inf Inf
b             -Inf Inf
sigma 2.220446e-16 Inf
tau   2.220446e-16 Inf

Slot "params":
[1] "b"     "lntau"

Slot "varnames":
[1] "a"     "b"     "sigma" "tau"

Slot "fixed":
a sigma
5.0   0.5

Slot "fancy":
[1] "Ricker" "Normal"

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<bytecode: 0x851f998>
<environment: 0x7f6eca8>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<bytecode: 0x80a1070>
<environment: 0x7f6eca8>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1], mean = logx[-T] + mle["a"] +
mle["b"] * exp(logx[-T]), sd = mle["sigma"], log = TRUE)
logd2 <- dnorm(data[-1], mean = logx[-1], sd = mle["tau"],
log = TRUE)
rval <- sum(logd1) + sum(logd2, na.rm = TRUE)
}
else {
stop("not yet implemented")
}
rval
}
<bytecode: 0x5191b80>
<environment: 0x7f6eca8>

Slot "neffective":
function (obs)
sum(!is.na(obs)) - 1
<bytecode: 0x2b88d70>
<environment: 0x7f6eca8>

An object of class "pvamodel"
Slot "growth.model":
[1] "thetalogistic"

Slot "obs.error":
[1] "none"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i],prcx)
mu[j,i] <- x[j-1,i] + r*( 1- (exp(x[j-1,i])/K)^theta )
}
}
r ~ dnorm(0,1)
K ~ dexp(0.005)
theta <- 1
sigma <- exp(lnsigma)
lnsigma ~ dnorm(0, 1)
prcx <- 1/sigma^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 4

Slot "support":
Min Max
r             -Inf Inf
K     2.220446e-16 Inf
theta         -Inf Inf
sigma 2.220446e-16 Inf

Slot "params":
[1] "r"       "K"       "lnsigma"

Slot "varnames":
[1] "r"     "K"     "theta" "sigma"

Slot "fixed":
theta
1

Slot "fancy":
[1] "Theta Logistic" NA

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("b" %in% vn)
mcmc[[i]][, "b"] <- atanh(mcmc[[i]][, "b"] + 1)
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("b" %in% vn)
vn[vn == "b"] <- "z"
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7c83500>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("z" %in% vn)
mcmc[[i]][, "z"] <- tanh(mcmc[[i]][, "z"]) - 1
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("z" %in% vn)
vn[vn == "z"] <- "b"
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7c83500>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(logx)
m <- which(is.na(data))
if (length(m) > 0) {
if (alt_obserror) {
stop("not yet implemented")
}
else {
y <- data
y[m] <- logx[m]
rval <- sum(dnorm(y[-1], mean = y[-T] + mle["r"] *
(1 - (exp(y[-T])/mle["K"])^mle["theta"]), sd = mle["sigma"],
log = TRUE))
}
}
else {
logd1 <- dnorm(data[-1], mean = data[-T] + mle["r"] *
(1 - (exp(data[-T])/mle["K"])^mle["theta"]), sd = mle["sigma"],
log = TRUE)
logd2 <- if (alt_obserror) {
dnorm(logx[-1], mean = logx[-T] + mle["r"] * (1 -
(exp(logx[-T])/mle["K"])^mle["theta"]), sd = mle["sigma"],
log = TRUE)
}
else {
0
}
rval <- sum(logd1) + sum(logd2)
}
rval
}
<environment: 0x7c83500>

Slot "neffective":
function (obs)
sum(!is.na(obs))
<environment: 0x7c83500>

An object of class "pvamodel"
Slot "growth.model":
[1] "bevertonholt"

Slot "obs.error":
[1] "normal"

Slot "model":
Object of class "custommodel":

model {
for (i in 1:kk) {
N[1,i] <- exp(y[1,i])
x[1,i] <- y[1,i]
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i],prcx)
mu[j,i] <- x[j-1,i] + r - log(1 + (N[j-1,i]/K)^theta)
N[j,i] <- max(exp(x[j,i]) , 1)
y[j,i]~dnorm(x[j,i],prcy)
}
}
r ~ dnorm(0,4)
K ~ dexp(0.01)
theta <- 1
sigma <- exp(lnsigma)
lnsigma ~ dnorm(0, 1)
prcx <- 1/sigma^2
tau <- exp(lntau)
lntau ~ dnorm(0, 1)
prcy <- 1/tau^2
}

Slot "genmodel":
function ()
NULL

Slot "p":
[1] 5

Slot "support":
Min Max
r             -Inf Inf
K     2.220446e-16 Inf
theta         -Inf Inf
sigma 2.220446e-16 Inf
tau   2.220446e-16 Inf

Slot "params":
[1] "r"       "K"       "lnsigma" "lntau"

Slot "varnames":
[1] "r"     "K"     "theta" "sigma" "tau"

Slot "fixed":
theta
1

Slot "fancy":
[1] "Generalized Beverton-Holt" "Normal"

Slot "transf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("b" %in% vn)
mcmc[[i]][, "b"] <- atanh(mcmc[[i]][, "b"] + 1)
if ("sigma" %in% vn)
mcmc[[i]][, "sigma"] <- log(mcmc[[i]][, "sigma"])
if ("tau" %in% vn)
mcmc[[i]][, "tau"] <- log(mcmc[[i]][, "tau"])
}
if ("b" %in% vn)
vn[vn == "b"] <- "z"
if ("sigma" %in% vn)
vn[vn == "sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn == "tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7cfcf38>

Slot "backtransf":
function (mcmc, obs.error)
{
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("z" %in% vn)
mcmc[[i]][, "z"] <- tanh(mcmc[[i]][, "z"]) - 1
if ("lnsigma" %in% vn)
mcmc[[i]][, "lnsigma"] <- exp(mcmc[[i]][, "lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][, "lntau"] <- exp(mcmc[[i]][, "lntau"])
}
if ("z" %in% vn)
vn[vn == "z"] <- "b"
if ("lnsigma" %in% vn)
vn[vn == "lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn == "lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
<environment: 0x7cfcf38>

Slot "logdensity":
function (logx, mle, data, null_obserror = FALSE, alt_obserror = FALSE)
{
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1], mean = logx[-T] + mle["r"] +
log(1 + (exp(logx[-T])/mle["K"])^mle["theta"]), sd = mle["sigma"],
log = TRUE)
logd2 <- dnorm(data[-1], mean = logx[-1], sd = mle["tau"],
log = TRUE)
rval <- sum(logd1) + sum(logd2, na.rm = TRUE)
}
else {
stop("not yet implemented")
}
rval
}
<environment: 0x7cfcf38>

Slot "neffective":
function (obs)
sum(!is.na(obs))
<environment: 0x7cfcf38>

Object of class "custommodel":

model {
for (i in 1:kk) {
N[1,i] <- exp(y[1,i])
x[1,i] <- y[1,i]
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]
N[j,i] <- min(exp(x[j,i]), 10000)
y[j,i] ~ dnorm(x[j,i], prcy)
}
}
sigma <- exp(lnsigma)
tau <- exp(lntau)
lnsigma ~ dnorm(0, 1)
lntau ~ dnorm(0, 1)
a ~ dnorm(2, 1)
b ~ dnorm(0, 10)
prcx <- 1/sigma^2
prcy <- 1/tau^2
}
Object of class "custommodel":

model {
for (i in 1:kk) {
N[1,i] <- exp(y[1,i])
x[1,i] <- y[1,i]
for (j in 2:T) {
x[j,i] ~ dnorm(mu[j,i], prcx)
mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]
N[j,i] <- min(exp(x[j,i]), 10000)
y[j,i] ~ dnorm(x[j,i], prcy)
}
}
sigma <- 0.5
tau <- exp(lntau)
lnsigma <- log(sigma)
lntau ~ dnorm(0, 1)
a ~ dnorm(2, 1)
b ~ dnorm(0, 10)
prcx <- 1/sigma^2
prcy <- 1/tau^2
}
```

PVAClone documentation built on May 2, 2019, 5:49 a.m.