Multi-species dynamic occupancy model with R and JAGS

[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post is intended to provide a simple example of how to construct
and make inferences on a multi-species multi-year occupancy model using
R, JAGS, and the ‘rjags’ package. This is not intended to be a
standalone tutorial on dynamic community occupancy modeling. Useful
primary literature references include MacKenzie et al. (2002), Kery and
Royle (2007), Royle and Kery (2007), Russell et al. (2009), and Dorazio
et al. (2010). Royle and Dorazio’s Heirarchichal Modeling and
Inference in Ecology
also provides a clear explanation of simple one
species occupancy models, multispecies occupancy models, and dynamic
(multiyear) occupancy models, among other things. There’s also a wealth
of code provided
here by
Elise Zipkin, J. Andrew Royle, and others.

Before getting started, we can define two convenience functions:

1
2
3
4
5
6
7
<span class="line">logit <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    log<span class="p">(</span>x<span class="o">/</span><span class="p">(</span><span class="m">1</span> <span class="o">-</span> x<span class="p">))</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">antilogit <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    exp<span class="p">(</span>x<span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="m">1</span> <span class="o">+</span> exp<span class="p">(</span>x<span class="p">))</span>
</span><span class="line"><span class="p">}</span>
</span>

Then, initializing the number of sites, species, years, and repeat surveys (i.e. surveys within years, where the occupancy status of a site
is assumed to be constant),

1
2
3
4
<span class="line">nsite <span class="o"><-</span> <span class="m">150</span>
</span><span class="line">nspec <span class="o"><-</span> <span class="m">6</span>
</span><span class="line">nyear <span class="o"><-</span> <span class="m">4</span>
</span><span class="line">nrep <span class="o"><-</span> <span class="m">3</span>
</span>

we can begin to consider occupancy. We’re interested in making
inferences about the rates of colonization and population persistence
for each species in a community, while estimating and accounting for
imperfect detection.

Occupancy status at site , by species , in year is
represented by . For occupied sites ; for
unoccupied sites . However, is incompletely observed: it
is possible that a species is present at a site in some year
() but species was never
seen at at site in year across all repeat
surveys because of imperfect detection. These observations
are represented by . Here we assume that there are no
“false positive” observations. In other words, if
, then . If a site is
occupied, the probability that is represented as
a Bernoulli trial with probability of
detection , such that

The occupancy status of species at site in year
is modeled as a Markov Bernoulli trial. In other words whether a species
is present at a site in year is influenced by whether it was
present at year .

where for

and in year one

where the occupancy status in year 0,
, and
. and are
parameters that control the probabilities of colonization and
persistence. If a site was unoccupied by species in a previous
year , then the probability of colonization is
given by the antilogit of . If a site was previously
occupied , the probability of population
persistence is given by the anitlogit of . We
assume that the distributions of species specific parameters are
defined by community level hyperparameters such that
and
. We can generate
occupancy data as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
<span class="line"><span class="c1"># community level hyperparameters</span>
</span><span class="line">p_beta <span class="o">=</span> <span class="m">0.7</span>
</span><span class="line">mubeta <span class="o"><-</span> logit<span class="p">(</span>p_beta<span class="p">)</span>
</span><span class="line">sdbeta <span class="o"><-</span> <span class="m">2</span>
</span><span class="line">
</span><span class="line">p_rho <span class="o"><-</span> <span class="m">0.8</span>
</span><span class="line">murho <span class="o"><-</span> logit<span class="p">(</span>p_rho<span class="p">)</span>
</span><span class="line">sdrho <span class="o"><-</span> <span class="m">1</span>
</span><span class="line">
</span><span class="line"><span class="c1"># species specific random effects</span>
</span><span class="line">set.seed<span class="p">(</span><span class="m">1</span><span class="p">)</span>  <span class="c1"># for reproducibility</span>
</span><span class="line">beta <span class="o"><-</span> rnorm<span class="p">(</span>nspec<span class="p">,</span> mubeta<span class="p">,</span> sdbeta<span class="p">)</span>
</span><span class="line">set.seed<span class="p">(</span><span class="m">1008</span><span class="p">)</span>
</span><span class="line">rho <span class="o"><-</span> rnorm<span class="p">(</span>nspec<span class="p">,</span> murho<span class="p">,</span> sdrho<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial occupancy states</span>
</span><span class="line">set.seed<span class="p">(</span><span class="m">237</span><span class="p">)</span>
</span><span class="line">rho0 <span class="o"><-</span> runif<span class="p">(</span>nspec<span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">)</span>
</span><span class="line">z0 <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nspec<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    z0<span class="p">[,</span> i<span class="p">]</span> <span class="o"><-</span> rbinom<span class="p">(</span>nsite<span class="p">,</span> <span class="m">1</span><span class="p">,</span> rho0<span class="p">[</span>i<span class="p">])</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># subsequent occupancy</span>
</span><span class="line">z <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">,</span> nyear<span class="p">))</span>
</span><span class="line">lpsi <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">,</span> nyear<span class="p">))</span>
</span><span class="line">psi <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">,</span> nyear<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>j <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nsite<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nspec<span class="p">)</span> <span class="p">{</span>
</span><span class="line">        <span class="kr">for</span> <span class="p">(</span>t <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nyear<span class="p">)</span> <span class="p">{</span>
</span><span class="line">            <span class="kr">if</span> <span class="p">(</span>t <span class="o">==</span> <span class="m">1</span><span class="p">)</span> <span class="p">{</span>
</span><span class="line">                lpsi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> beta<span class="p">[</span>i<span class="p">]</span> <span class="o">+</span> rho<span class="p">[</span>i<span class="p">]</span> <span class="o">*</span> z0<span class="p">[</span>j<span class="p">,</span> i<span class="p">]</span>
</span><span class="line">                psi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> antilogit<span class="p">(</span>lpsi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">])</span>
</span><span class="line">                z<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> rbinom<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> psi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">])</span>
</span><span class="line">            <span class="p">}</span> <span class="kr">else</span> <span class="p">{</span>
</span><span class="line">                lpsi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> beta<span class="p">[</span>i<span class="p">]</span> <span class="o">+</span> rho<span class="p">[</span>i<span class="p">]</span> <span class="o">*</span> z<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t <span class="o">-</span> <span class="m">1</span><span class="p">]</span>
</span><span class="line">                psi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> antilogit<span class="p">(</span>lpsi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">])</span>
</span><span class="line">                z<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> rbinom<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> psi<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">])</span>
</span><span class="line">            <span class="p">}</span>
</span><span class="line">        <span class="p">}</span>
</span><span class="line">    <span class="p">}</span>
</span><span class="line"><span class="p">}</span>
</span>

For simplicity, we’ll assume that there are no differences in species detectability among sites, years, or repeat surveys, but that detectability varies among species. We’ll again use hyperparameters to specify a distribution of detection probabilities in our community, such that $logit(p_i) \sim Normal(\mu_p, \sigma_p^2)$.

1
2
3
4
5
6
<span class="line">p_p <span class="o"><-</span> <span class="m">0.7</span>
</span><span class="line">mup <span class="o"><-</span> logit<span class="p">(</span>p_p<span class="p">)</span>
</span><span class="line">sdp <span class="o"><-</span> <span class="m">1.5</span>
</span><span class="line">set.seed<span class="p">(</span><span class="m">222</span><span class="p">)</span>
</span><span class="line">lp <span class="o"><-</span> rnorm<span class="p">(</span>nspec<span class="p">,</span> mup<span class="p">,</span> sdp<span class="p">)</span>
</span><span class="line">p <span class="o"><-</span> antilogit<span class="p">(</span>lp<span class="p">)</span>
</span>

We can now generate our observations based on occupancy states and
detection probabilities. Although this could be vectorized for speed,
let’s stick with nested for loops in the interest of clarity.

1
2
3
4
5
6
7
8
9
10
<span class="line">x <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">,</span> nyear<span class="p">,</span> nrep<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>j <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nsite<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nspec<span class="p">)</span> <span class="p">{</span>
</span><span class="line">        <span class="kr">for</span> <span class="p">(</span>t <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nyear<span class="p">)</span> <span class="p">{</span>
</span><span class="line">            <span class="kr">for</span> <span class="p">(</span>k <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nrep<span class="p">)</span> <span class="p">{</span>
</span><span class="line">                x<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">,</span> k<span class="p">]</span> <span class="o"><-</span> rbinom<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> p<span class="p">[</span>i<span class="p">]</span> <span class="o">*</span> z<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">])</span>
</span><span class="line">            <span class="p">}</span>
</span><span class="line">        <span class="p">}</span>
</span><span class="line">    <span class="p">}</span>
</span><span class="line"><span class="p">}</span>
</span>

Now that we’ve collected some data, we can specify our model:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
<span class="line">cat<span class="p">(</span><span class="s">"</span>
</span><span class="line"><span class="s">model{</span>
</span><span class="line"><span class="s">  #### priors</span>
</span><span class="line"><span class="s">  # beta hyperparameters</span>
</span><span class="line"><span class="s">  p_beta ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">	mubeta <- log(p_beta / (1 - p_beta))</span>
</span><span class="line"><span class="s">  sigmabeta ~ dunif(0, 10)</span>
</span><span class="line"><span class="s">  taubeta <- (1 / (sigmabeta * sigmabeta))</span>
</span><span class="line">
</span><span class="line"><span class="s">  # rho hyperparameters</span>
</span><span class="line"><span class="s">  p_rho ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">	murho <- log(p_rho / (1 - p_rho))</span>
</span><span class="line"><span class="s">	sigmarho~dunif(0,10)</span>
</span><span class="line"><span class="s">	taurho<-1/(sigmarho*sigmarho)</span>
</span><span class="line">
</span><span class="line"><span class="s">  # p hyperparameters</span>
</span><span class="line"><span class="s">  p_p ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">  mup <- log(p_p / (1 - p_p))</span>
</span><span class="line"><span class="s">  sigmap ~ dunif(0,10)</span>
</span><span class="line"><span class="s">  taup <- (1 / (sigmap * sigmap))</span>
</span><span class="line">
</span><span class="line"><span class="s">  #### occupancy model</span>
</span><span class="line"><span class="s">  # species specific random effects</span>
</span><span class="line"><span class="s">  for (i in 1:(nspec)) {</span>
</span><span class="line"><span class="s">    rho0[i] ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    beta[i] ~ dnorm(mubeta, taubeta)</span>
</span><span class="line"><span class="s">    rho[i] ~ dnorm(murho, taurho)</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line"><span class="s">  </span>
</span><span class="line"><span class="s">  # occupancy states</span>
</span><span class="line"><span class="s">  for (j in 1:nsite) {</span>
</span><span class="line"><span class="s">    for (i in 1:nspec) {</span>
</span><span class="line"><span class="s">      z0[j, i] ~ dbern(rho0[i])</span>
</span><span class="line"><span class="s">      logit(psi[j, i, 1]) <- beta[i] + rho[i] * z0[j, i] </span>
</span><span class="line"><span class="s">      z[j, i, 1] ~ dbern(psi[j, i, 1]) </span>
</span><span class="line"><span class="s">      for (t in 2:nyear) {</span>
</span><span class="line"><span class="s">        logit(psi[j, i, t]) <- beta[i] + rho[i] * z[j, i, t-1]</span>
</span><span class="line"><span class="s">        z[j, i, t] ~ dbern(psi[j, i, t])</span>
</span><span class="line"><span class="s">      }</span>
</span><span class="line"><span class="s">    }</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line">
</span><span class="line"><span class="s">  #### detection model</span>
</span><span class="line"><span class="s">	for(i in 1:nspec){ </span>
</span><span class="line"><span class="s">		lp[i] ~ dnorm(mup, taup)</span>
</span><span class="line"><span class="s">		p[i] <- (exp(lp[i])) / (1 + exp(lp[i]))</span>
</span><span class="line"><span class="s">	}</span>
</span><span class="line"><span class="s">	</span>
</span><span class="line"><span class="s">  #### observation model</span>
</span><span class="line"><span class="s">  for (j in 1:nsite){</span>
</span><span class="line"><span class="s">    for (i in 1:nspec){</span>
</span><span class="line"><span class="s">      for (t in 1:nyear){</span>
</span><span class="line"><span class="s">        mu[j, i, t] <- z[j, i, t] * p[i] </span>
</span><span class="line"><span class="s">        for (k in 1:nrep){</span>
</span><span class="line"><span class="s">          x[j, i, t, k] ~ dbern(mu[j, i, t])</span>
</span><span class="line"><span class="s">        }</span>
</span><span class="line"><span class="s">      }</span>
</span><span class="line"><span class="s">    }</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line"><span class="s">}</span>
</span><span class="line"><span class="s">"</span><span class="p">,</span> fill<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> file<span class="o">=</span><span class="s">"com_occ.txt"</span><span class="p">)</span>
</span>

Next, bundle up the data.

1
<span class="line">data <span class="o"><-</span> list<span class="p">(</span>x <span class="o">=</span> x<span class="p">,</span> nrep <span class="o">=</span> nrep<span class="p">,</span> nsite <span class="o">=</span> nsite<span class="p">,</span> nspec <span class="o">=</span> nspec<span class="p">,</span> nyear <span class="o">=</span> nyear<span class="p">)</span>
</span>

Provide initial values.

1
2
3
4
5
6
7
8
9
10
11
12
13
<span class="line">zinit <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nsite<span class="p">,</span> nspec<span class="p">,</span> nyear<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>j <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nsite<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nspec<span class="p">)</span> <span class="p">{</span>
</span><span class="line">        <span class="kr">for</span> <span class="p">(</span>t <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nyear<span class="p">)</span> <span class="p">{</span>
</span><span class="line">            zinit<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> max<span class="p">(</span>x<span class="p">[</span>j<span class="p">,</span> i<span class="p">,</span> t<span class="p">,</span> <span class="p">])</span>
</span><span class="line">        <span class="p">}</span>
</span><span class="line">    <span class="p">}</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">inits <span class="o"><-</span> <span class="kr">function</span><span class="p">()</span> <span class="p">{</span>
</span><span class="line">    list<span class="p">(</span>p_beta <span class="o">=</span> runif<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">),</span> p_rho <span class="o">=</span> runif<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">),</span> sigmarho <span class="o">=</span> runif<span class="p">(</span><span class="m">1</span><span class="p">,</span>
</span><span class="line">        <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">),</span> sigmap <span class="o">=</span> runif<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">10</span><span class="p">),</span> sigmabeta <span class="o">=</span> runif<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">10</span><span class="p">),</span> z <span class="o">=</span> zinit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span>

As a side note, it is helpful in JAGS to provide initial values for the incompletely observed occupancy state $z$ that are consistent with observed presences, as provided in this example with zinit. In other words if $x(j, i, t, k)=1$, provide an intial value of $1$ for $z(j, i, t)$. Unlike WinBUGS and OpenBUGS, if you do not do this, you’ll often (but not always) encounter an error message such as:

1
2
3
<span class="line"><span class="c1"># Error in jags.model(file = 'com_occ.txt', data = data, n.chains = 3) :</span>
</span><span class="line"><span class="c1"># Error in node x[1,1,2,3] Observed node inconsistent with unobserved</span>
</span><span class="line"><span class="c1"># parents at initialization</span>
</span>

Now we’re ready to monitor and make inferences about some parameters of interest using JAGS.

1
2
3
4
<span class="line">params <span class="o"><-</span> c<span class="p">(</span><span class="s">"lp"</span><span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> <span class="s">"rho"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">require<span class="p">(</span>rjags<span class="p">)</span>
</span><span class="line">ocmod <span class="o"><-</span> jags.model<span class="p">(</span>file <span class="o">=</span> <span class="s">"com_occ.txt"</span><span class="p">,</span> inits <span class="o">=</span> inits<span class="p">,</span> data <span class="o">=</span> data<span class="p">,</span> n.chains <span class="o">=</span> <span class="m">3</span><span class="p">)</span>
</span>

1
2
3
4
<span class="line">nburn <span class="o"><-</span> <span class="m">2000</span>
</span><span class="line">update<span class="p">(</span>ocmod<span class="p">,</span> n.iter <span class="o">=</span> nburn<span class="p">)</span>
</span><span class="line">out <span class="o"><-</span> coda.samples<span class="p">(</span>ocmod<span class="p">,</span> n.iter <span class="o">=</span> <span class="m">17000</span><span class="p">,</span> variable.names <span class="o">=</span> params<span class="p">)</span>
</span><span class="line">summary<span class="p">(</span>out<span class="p">)</span>
</span>

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
<span class="line"><span class="c1">## </span>
</span><span class="line"><span class="c1">## Iterations = 3001:20000</span>
</span><span class="line"><span class="c1">## Thinning interval = 1 </span>
</span><span class="line"><span class="c1">## Number of chains = 3 </span>
</span><span class="line"><span class="c1">## Sample size per chain = 17000 </span>
</span><span class="line"><span class="c1">## </span>
</span><span class="line"><span class="c1">## 1. Empirical mean and standard deviation for each variable,</span>
</span><span class="line"><span class="c1">##    plus standard error of the mean:</span>
</span><span class="line"><span class="c1">## </span>
</span><span class="line"><span class="c1">##            Mean     SD Naive SE Time-series SE</span>
</span><span class="line"><span class="c1">## beta[1] -0.4539 0.1629 0.000721       0.002140</span>
</span><span class="line"><span class="c1">## beta[2]  0.9986 0.4274 0.001893       0.011038</span>
</span><span class="line"><span class="c1">## beta[3] -0.9036 0.1157 0.000513       0.001006</span>
</span><span class="line"><span class="c1">## beta[4]  4.7439 2.5492 0.011288       0.050756</span>
</span><span class="line"><span class="c1">## beta[5]  1.5512 0.3251 0.001440       0.007095</span>
</span><span class="line"><span class="c1">## beta[6] -0.8833 0.3021 0.001338       0.005858</span>
</span><span class="line"><span class="c1">## lp[1]    3.0660 0.1400 0.000620       0.000782</span>
</span><span class="line"><span class="c1">## lp[2]    0.8533 0.0580 0.000257       0.000409</span>
</span><span class="line"><span class="c1">## lp[3]    2.9045 0.2007 0.000889       0.001160</span>
</span><span class="line"><span class="c1">## lp[4]    0.3773 0.0491 0.000217       0.000315</span>
</span><span class="line"><span class="c1">## lp[5]    1.2224 0.0635 0.000281       0.000407</span>
</span><span class="line"><span class="c1">## lp[6]    0.4881 0.0615 0.000272       0.000551</span>
</span><span class="line"><span class="c1">## rho[1]   1.8640 0.2200 0.000974       0.002691</span>
</span><span class="line"><span class="c1">## rho[2]   2.6201 0.6009 0.002661       0.011576</span>
</span><span class="line"><span class="c1">## rho[3]  -0.0766 0.2229 0.000987       0.002010</span>
</span><span class="line"><span class="c1">## rho[4]   2.6567 2.2459 0.009945       0.035121</span>
</span><span class="line"><span class="c1">## rho[5]   1.1056 0.4091 0.001812       0.007679</span>
</span><span class="line"><span class="c1">## rho[6]   3.4425 0.4591 0.002033       0.009132</span>
</span><span class="line"><span class="c1">## </span>
</span><span class="line"><span class="c1">## 2. Quantiles for each variable:</span>
</span><span class="line"><span class="c1">## </span>
</span><span class="line"><span class="c1">##           2.5%    25%     50%    75%  97.5%</span>
</span><span class="line"><span class="c1">## beta[1] -0.777 -0.564 -0.4515 -0.343 -0.139</span>
</span><span class="line"><span class="c1">## beta[2]  0.168  0.707  0.9975  1.288  1.833</span>
</span><span class="line"><span class="c1">## beta[3] -1.132 -0.981 -0.9028 -0.826 -0.679</span>
</span><span class="line"><span class="c1">## beta[4]  1.117  3.206  4.3157  5.732 11.010</span>
</span><span class="line"><span class="c1">## beta[5]  0.896  1.338  1.5590  1.770  2.174</span>
</span><span class="line"><span class="c1">## beta[6] -1.509 -1.077 -0.8699 -0.675 -0.327</span>
</span><span class="line"><span class="c1">## lp[1]    2.803  2.970  3.0630  3.159  3.348</span>
</span><span class="line"><span class="c1">## lp[2]    0.740  0.814  0.8535  0.892  0.967</span>
</span><span class="line"><span class="c1">## lp[3]    2.528  2.766  2.8992  3.035  3.316</span>
</span><span class="line"><span class="c1">## lp[4]    0.282  0.344  0.3771  0.410  0.474</span>
</span><span class="line"><span class="c1">## lp[5]    1.100  1.179  1.2216  1.265  1.348</span>
</span><span class="line"><span class="c1">## lp[6]    0.369  0.446  0.4877  0.529  0.610</span>
</span><span class="line"><span class="c1">## rho[1]   1.439  1.714  1.8611  2.011  2.301</span>
</span><span class="line"><span class="c1">## rho[2]   1.504  2.205  2.5972  3.013  3.843</span>
</span><span class="line"><span class="c1">## rho[3]  -0.516 -0.226 -0.0747  0.074  0.359</span>
</span><span class="line"><span class="c1">## rho[4]  -0.928  1.270  2.3693  3.688  8.060</span>
</span><span class="line"><span class="c1">## rho[5]   0.319  0.831  1.1003  1.380  1.920</span>
</span><span class="line"><span class="c1">## rho[6]   2.601  3.123  3.4209  3.741  4.390</span>
</span>

1
<span class="line">plot<span class="p">(</span>out<span class="p">)</span>
</span>

plot of chunk unnamed-chunk-10
plot of chunk unnamed-chunk-10
plot of chunk unnamed-chunk-10
plot of chunk unnamed-chunk-10
plot of chunk unnamed-chunk-10

At this point, you’ll want to run through the usual MCMC diagnostics to
check for convergence and adjust the burn-in or number of iterations
accordingly. Once satisfied, we can check to see how well our model
performed based on our known parameter values.

1
2
3
<span class="line">require<span class="p">(</span>mcmcplots<span class="p">)</span>
</span><span class="line">caterplot<span class="p">(</span>out<span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> style <span class="o">=</span> <span class="s">"plain"</span><span class="p">)</span>
</span><span class="line">caterpoints<span class="p">(</span>beta<span class="p">)</span>
</span>

plot of chunk unnamed-chunk-11

1
2
<span class="line">caterplot<span class="p">(</span>out<span class="p">,</span> <span class="s">"lp"</span><span class="p">,</span> style <span class="o">=</span> <span class="s">"plain"</span><span class="p">)</span>
</span><span class="line">caterpoints<span class="p">(</span>lp<span class="p">)</span>
</span>

plot of chunk unnamed-chunk-11

1
2
<span class="line">caterplot<span class="p">(</span>out<span class="p">,</span> <span class="s">"rho"</span><span class="p">,</span> style <span class="o">=</span> <span class="s">"plain"</span><span class="p">)</span>
</span><span class="line">caterpoints<span class="p">(</span>rho<span class="p">)</span>
</span>

plot of chunk unnamed-chunk-11

To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)