Description Usage Arguments Details Value Author(s) References See Also Examples
Population growth model to be used in model fitting
via pva
.
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)
|
obs.error |
Character, describing the observation error.
Can be |
fixed |
Named numeric vector or list with fixed parameter names and values. Can be used for providing alternative prior specifications, see Details and Examples. |
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.
An S4 class of 'pvamodel' (see pvamodel-class
)
Khurram Nadeem and Peter Solymos
Nadeem, K., Lele S. R., 2012. Likelihood based population viability analysis in the presence of observation error. Oikos 121, 1656–1664.
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
|
Loading required package: dcmle
Loading required package: dclone
Loading required package: coda
Loading required package: parallel
Loading required package: Matrix
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
Loading required package: stats4
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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.