Dynamic occupancy models in Stan

[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.

Occupancy modeling is possible in Stan as shown here, despite the lack of support for integer parameters.
In many Bayesian applications of occupancy modeling, the true occupancy states (0 or 1) are directly modeled, but this can be avoided by marginalizing out the true occupancy state.
The Stan manual (pg. 96) gives an example of this kind of marginalization for a discrete change-point model.

Here’s a Stan implementation of a dynamic (multi-year) occupancy model of the sort described by MacKenzie et al. (2003).

First, the model statement:

dyn_occ.stan
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
<span class="line"><span class="n">data</span> <span class="p">{</span>
</span><span class="line">   <span class="kt">int</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="o">></span> <span class="n">nsite</span><span class="p">;</span>
</span><span class="line">   <span class="kt">int</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="o">></span> <span class="n">nyear</span><span class="p">;</span>
</span><span class="line">   <span class="kt">int</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="o">></span> <span class="n">nrep</span><span class="p">;</span>
</span><span class="line">   <span class="kt">int</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">upper</span><span class="o">=</span><span class="mi">1</span><span class="o">></span> <span class="n">Y</span><span class="p">[</span><span class="n">nsite</span><span class="p">,</span><span class="n">nyear</span><span class="p">,</span><span class="n">nrep</span><span class="p">];</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line"><span class="n">parameters</span> <span class="p">{</span>
</span><span class="line">   <span class="n">real</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">upper</span><span class="o">=</span><span class="mi">1</span><span class="o">></span> <span class="n">p</span><span class="p">;</span>
</span><span class="line">   <span class="n">real</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">upper</span><span class="o">=</span><span class="mi">1</span><span class="o">></span> <span class="n">gamma</span><span class="p">;</span>
</span><span class="line">   <span class="n">real</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">upper</span><span class="o">=</span><span class="mi">1</span><span class="o">></span> <span class="n">phi</span><span class="p">;</span>
</span><span class="line">   <span class="n">real</span><span class="o"><</span><span class="n">lower</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">upper</span><span class="o">=</span><span class="mi">1</span><span class="o">></span> <span class="n">psi1</span><span class="p">;</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line"><span class="n">transformed</span> <span class="n">parameters</span> <span class="p">{</span>
</span><span class="line">   <span class="n">matrix</span><span class="p">[</span><span class="n">nsite</span><span class="p">,</span> <span class="n">nyear</span><span class="p">]</span> <span class="n">psi</span><span class="p">;</span>
</span><span class="line">   <span class="k">for</span> <span class="p">(</span><span class="n">r</span> <span class="n">in</span> <span class="mi">1</span><span class="o">:</span><span class="n">nsite</span><span class="p">){</span>
</span><span class="line">     <span class="k">for</span> <span class="p">(</span><span class="n">t</span> <span class="n">in</span> <span class="mi">1</span><span class="o">:</span><span class="n">nyear</span><span class="p">){</span>
</span><span class="line">       <span class="k">if</span> <span class="p">(</span><span class="n">t</span> <span class="o"><</span> <span class="mi">2</span><span class="p">){</span>
</span><span class="line">          <span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span> <span class="o"><-</span> <span class="n">psi1</span><span class="p">;</span>
</span><span class="line">       <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
</span><span class="line">          <span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span> <span class="o"><-</span> <span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">phi</span> <span class="o">+</span> <span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="o">-</span><span class="mi">1</span><span class="p">])</span> <span class="o">*</span> <span class="n">gamma</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 class="p">}</span>
</span><span class="line"><span class="n">model</span> <span class="p">{</span>
</span><span class="line">   <span class="c1">// priors</span>
</span><span class="line">  <span class="n">psi1</span> <span class="o">~</span> <span class="n">uniform</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">);</span>
</span><span class="line">  <span class="n">gamma</span> <span class="o">~</span> <span class="n">uniform</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">);</span>
</span><span class="line">  <span class="n">phi</span> <span class="o">~</span> <span class="n">uniform</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">);</span>
</span><span class="line">  <span class="n">p</span> <span class="o">~</span> <span class="n">uniform</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">);</span>
</span><span class="line">
</span><span class="line">   <span class="c1">// likelihood</span>
</span><span class="line">  <span class="k">for</span> <span class="p">(</span><span class="n">r</span> <span class="n">in</span> <span class="mi">1</span><span class="o">:</span><span class="n">nsite</span><span class="p">){</span>
</span><span class="line">    <span class="k">for</span> <span class="p">(</span><span class="n">t</span> <span class="n">in</span> <span class="mi">1</span><span class="o">:</span><span class="n">nyear</span><span class="p">){</span>
</span><span class="line">      <span class="k">if</span> <span class="p">(</span><span class="n">sum</span><span class="p">(</span><span class="n">Y</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">])</span> <span class="o">></span> <span class="mi">0</span><span class="p">){</span>
</span><span class="line">        <span class="n">increment_log_prob</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">])</span> <span class="o">+</span> <span class="n">bernoulli_log</span><span class="p">(</span><span class="n">Y</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">],</span> <span class="n">p</span><span class="p">));</span>
</span><span class="line">      <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
</span><span class="line">        <span class="n">increment_log_prob</span><span class="p">(</span><span class="n">log_sum_exp</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">])</span> <span class="o">+</span> <span class="n">bernoulli_log</span><span class="p">(</span><span class="n">Y</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">],</span><span class="n">p</span><span class="p">),</span>
</span><span class="line">                                      <span class="n">log</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">psi</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">t</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 class="p">}</span>
</span>

This model can be made faster by storing values for log(psi) and log(1 – psi), as done in Bob Carpenter’s single season example.

Fitting the model (in parallel):

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
<span class="line">library<span class="p">(</span>rstan<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>parallel<span class="p">)</span>
</span><span class="line">
</span><span class="line">simulate_data <span class="o"><-</span> <span class="kr">function</span><span class="p">(){</span>
</span><span class="line">  nsite <span class="o"><-</span> <span class="m">100</span><span class="p">;</span>
</span><span class="line">  nrep <span class="o"><-</span> <span class="m">2</span><span class="p">;</span>     <span class="c1"># repeat surveys</span>
</span><span class="line">  nyear <span class="o"><-</span> <span class="m">10</span><span class="p">;</span>  <span class="c1"># nyears</span>
</span><span class="line">  p <span class="o"><-</span> <span class="m">0.8</span><span class="p">;</span>
</span><span class="line">  gamma <span class="o"><-</span> <span class="m">.2</span>
</span><span class="line">  phi <span class="o"><-</span> <span class="m">.8</span>
</span><span class="line">  psi1 <span class="o"><-</span> <span class="m">.8</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> nyear<span class="p">))</span>
</span><span class="line">  psi<span class="p">[,</span> <span class="m">1</span><span class="p">]</span> <span class="o"><-</span> psi1
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>t <span class="kr">in</span> <span class="m">2</span><span class="o">:</span>nyear<span class="p">){</span>
</span><span class="line">    psi<span class="p">[,</span> t<span class="p">]</span> <span class="o"><-</span> psi<span class="p">[,</span> t <span class="o">-</span> <span class="m">1</span><span class="p">]</span> <span class="o">*</span> phi <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">-</span> psi<span class="p">[,</span> t<span class="m">-1</span><span class="p">])</span> <span class="o">*</span> gamma
</span><span class="line">  <span class="p">}</span>
</span><span class="line">
</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> nyear<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><span class="line">    Z<span class="p">[,</span> t<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> psi<span class="p">[,</span> t<span class="p">])</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line">
</span><span class="line">  Y <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nsite<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>r <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nsite<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><span class="line">      Y<span class="p">[</span>r<span class="p">,</span> t<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> rbinom<span class="p">(</span>nrep<span class="p">,</span> <span class="m">1</span><span class="p">,</span> Z<span class="p">[</span>r<span class="p">,</span> t<span class="p">]</span> <span class="o">*</span> p<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="kr">return</span><span class="p">(</span>list<span class="p">(</span>nsite<span class="o">=</span>nsite<span class="p">,</span> nrep<span class="o">=</span>nrep<span class="p">,</span> nyear<span class="o">=</span>nyear<span class="p">,</span>
</span><span class="line">              p<span class="o">=</span>p<span class="p">,</span> gamma<span class="o">=</span>gamma<span class="p">,</span> phi<span class="o">=</span>phi<span class="p">,</span>
</span><span class="line">              psi1<span class="o">=</span>psi1<span class="p">,</span> Y<span class="o">=</span>Y<span class="p">))</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1">## parallel run</span>
</span><span class="line">d <span class="o"><-</span> simulate_data<span class="p">()</span>
</span><span class="line"><span class="c1"># initialize model</span>
</span><span class="line">mod <span class="o"><-</span> stan<span class="p">(</span><span class="s">"dyn_occ.stan"</span><span class="p">,</span> data <span class="o">=</span> d<span class="p">,</span> chains <span class="o">=</span> <span class="m">0</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">estimate_params <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>model<span class="p">,</span> data<span class="p">,</span> iter<span class="o">=</span><span class="m">2000</span><span class="p">){</span>
</span><span class="line">  sflist <span class="o"><-</span> mclapply<span class="p">(</span><span class="m">1</span><span class="o">:</span><span class="m">4</span><span class="p">,</span> mc.cores <span class="o">=</span> <span class="m">4</span><span class="p">,</span>
</span><span class="line">             <span class="kr">function</span><span class="p">(</span>i<span class="p">)</span> stan<span class="p">(</span>fit <span class="o">=</span> model<span class="p">,</span> data <span class="o">=</span> data<span class="p">,</span>
</span><span class="line">                              chains <span class="o">=</span> <span class="m">1</span><span class="p">,</span> chain_id <span class="o">=</span> i<span class="p">,</span>
</span><span class="line">                              iter<span class="o">=</span>iter<span class="p">,</span>
</span><span class="line">                              refresh <span class="o">=</span> <span class="m">-1</span><span class="p">))</span>
</span><span class="line">  fit <span class="o"><-</span> sflist2stanfit<span class="p">(</span>sflist<span class="p">)</span>
</span><span class="line">  <span class="kr">return</span><span class="p">(</span>fit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">fit <span class="o"><-</span> estimate_params<span class="p">(</span>model<span class="o">=</span>mod<span class="p">,</span> data<span class="o">=</span>d<span class="p">)</span>
</span>

Does it work? Let’s whip up 1000 simulated datasets and their corresponding estimates for colonization and extinction rates.

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
<span class="line">one_run <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>initialized_mod<span class="p">,</span> params<span class="p">){</span>
</span><span class="line">  require<span class="p">(</span>modeest<span class="p">)</span>
</span><span class="line">  require<span class="p">(</span>reshape2<span class="p">)</span>
</span><span class="line">  source<span class="p">(</span><span class="s">"HDI.R"</span><span class="p">)</span>
</span><span class="line">  d <span class="o"><-</span> simulate_data<span class="p">()</span>
</span><span class="line">  fit <span class="o"><-</span> estimate_params<span class="p">(</span>model<span class="o">=</span>initialized_mod<span class="p">,</span> data<span class="o">=</span>d<span class="p">)</span>
</span><span class="line">  <span class="c1"># store HDI for certain params</span>
</span><span class="line">  post <span class="o"><-</span> extract<span class="p">(</span>fit<span class="p">,</span> params<span class="p">)</span>
</span><span class="line">  vals <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>length<span class="p">(</span>params<span class="p">),</span> <span class="m">3</span><span class="p">))</span>
</span><span class="line">  rownames<span class="p">(</span>vals<span class="p">)</span> <span class="o"><-</span> params
</span><span class="line">  colnames<span class="p">(</span>vals<span class="p">)</span> <span class="o"><-</span> c<span class="p">(</span><span class="s">"LCI"</span><span class="p">,</span> <span class="s">"mode"</span><span class="p">,</span> <span class="s">"UCI"</span><span class="p">)</span>
</span><span class="line">  mode <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> length<span class="p">(</span>params<span class="p">))</span>
</span><span class="line">  lci <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> length<span class="p">(</span>params<span class="p">))</span>
</span><span class="line">  uci <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> length<span class="p">(</span>params<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>length<span class="p">(</span>params<span class="p">)){</span>
</span><span class="line">    <span class="c1"># calculate posterior mode</span>
</span><span class="line">    mode<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> mlv<span class="p">(</span>post<span class="p">[[</span>i<span class="p">]],</span> method<span class="o">=</span><span class="s">"mfv"</span><span class="p">)</span><span class="o">$</span>M
</span><span class="line">    <span class="c1"># calculate HDI</span>
</span><span class="line">    interv <span class="o"><-</span> HDI<span class="p">(</span>post<span class="p">[[</span>i<span class="p">]])</span>
</span><span class="line">    lci<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> interv<span class="p">[</span><span class="m">1</span><span class="p">]</span>
</span><span class="line">    uci<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> interv<span class="p">[</span><span class="m">2</span><span class="p">]</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line">  vals <span class="o"><-</span> data.frame<span class="p">(</span>parameter <span class="o">=</span> params<span class="p">,</span> mode<span class="o">=</span>mode<span class="p">,</span>
</span><span class="line">                     lci<span class="o">=</span>lci<span class="p">,</span> uci<span class="o">=</span>uci<span class="p">)</span>
</span><span class="line">  <span class="kr">return</span><span class="p">(</span>vals<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">iterations <span class="o"><-</span> <span class="m">1000</span>
</span><span class="line">check <span class="o"><-</span> one_run<span class="p">(</span>mod<span class="p">,</span> c<span class="p">(</span><span class="s">"gamma"</span><span class="p">,</span> <span class="s">"phi"</span><span class="p">))</span>
</span><span class="line">check<span class="o">$</span>iteration <span class="o"><-</span> rep<span class="p">(</span><span class="m">1</span><span class="p">,</span> nrow<span class="p">(</span>check<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">2</span><span class="o">:</span>iterations<span class="p">){</span>
</span><span class="line">  new_d <span class="o"><-</span> one_run<span class="p">(</span>mod<span class="p">,</span> c<span class="p">(</span><span class="s">"gamma"</span><span class="p">,</span> <span class="s">"phi"</span><span class="p">))</span>
</span><span class="line">  new_d<span class="o">$</span>iteration <span class="o"><-</span> rep<span class="p">(</span>i<span class="p">,</span> nrow<span class="p">(</span>new_d<span class="p">))</span>
</span><span class="line">  check <span class="o"><-</span> rbind<span class="p">(</span>check<span class="p">,</span> new_d<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">true.vals <span class="o"><-</span> data.frame<span class="p">(</span>parameter <span class="o">=</span> c<span class="p">(</span><span class="s">"gamma"</span><span class="p">,</span> <span class="s">"phi"</span><span class="p">),</span>
</span><span class="line">                        value <span class="o">=</span> c<span class="p">(</span>d<span class="o">$</span>gamma<span class="p">,</span> d<span class="o">$</span>phi<span class="p">))</span>
</span><span class="line">post.vals <span class="o"><-</span> data.frame<span class="p">(</span>parameter <span class="o">=</span> c<span class="p">(</span><span class="s">"gamma"</span><span class="p">,</span> <span class="s">"phi"</span><span class="p">),</span>
</span><span class="line">                        value <span class="o">=</span> c<span class="p">(</span>mean<span class="p">(</span>check<span class="o">$</span>mode<span class="p">[</span>check<span class="o">$</span>parameter <span class="o">==</span> <span class="s">"gamma"</span><span class="p">]),</span>
</span><span class="line">                                  mean<span class="p">(</span>check<span class="o">$</span>mode<span class="p">[</span>check<span class="o">$</span>parameter <span class="o">==</span> <span class="s">"phi"</span><span class="p">])))</span>
</span><span class="line">
</span><span class="line">ggplot<span class="p">(</span>check<span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_point<span class="p">(</span>aes<span class="p">(</span>x<span class="o">=</span>mode<span class="p">,</span> y<span class="o">=</span>iteration<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  facet_wrap<span class="p">(</span><span class="o">~</span>parameter<span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_segment<span class="p">(</span>aes<span class="p">(</span>x<span class="o">=</span>lci<span class="p">,</span> xend<span class="o">=</span>uci<span class="p">,</span>
</span><span class="line">                   y<span class="o">=</span>iteration<span class="p">,</span> yend<span class="o">=</span>iteration<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  theme_bw<span class="p">()</span> <span class="o">+</span>
</span><span class="line">  geom_vline<span class="p">(</span>aes<span class="p">(</span>xintercept<span class="o">=</span>value<span class="p">),</span> true.vals<span class="p">,</span>
</span><span class="line">             color<span class="o">=</span><span class="s">"blue"</span><span class="p">,</span> linetype<span class="o">=</span><span class="s">"dashed"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_vline<span class="p">(</span>aes<span class="p">(</span>xintercept<span class="o">=</span>value<span class="p">),</span> post.vals<span class="p">,</span>
</span><span class="line">             color<span class="o">=</span><span class="s">"red"</span><span class="p">,</span> linetype<span class="o">=</span><span class="s">"dashed"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"value"</span><span class="p">)</span>
</span>

Here are the results for the probability of colonization $gamma$, and the probability of persistence $phi$.
The blue dashed line shows the true value, and the dashed red lines shows the mean of all 1000 posterior modes.
The black lines represent the HPDI for each interation, and the black points represent the posterior modes.
This example uses a uniform prior on both of these parameters – probably an overrepresentation of prior ignorance in most real systems.

Based on some initial exploration, this approach seems much much (much?) faster than explicitly modeling the latent occurrence states in JAGS, with better chain mixing and considerably less autocorrelation. Extension to multi-species models should be straightforward too – huzzah!

References

MacKenzie DI, Nichols JD, Hines JE, Knutson MG, Franklin AB. 2003. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84(8): 2200-2207. pdf

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)