Occupancy model fit & AUC

[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 models are used to understand species distributions while
accounting for imperfect detection. In this post, I’ll demonstrate a
method to evaluate the performance of occupancy models based on the
area under a receiver operating characteristic curve (AUC), as published
last year by Elise Zipkin and colleagues in
Ecological Applications.

Suppose we are to fit a multi-year occupancy model for one species. We
will evaluate the fit based on how well the model predicts occupancy in
the final year of the project. Start by simulating some data
(for details on the structure of these simulated data, refer to
this post and
references therein):

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
<span class="line"><span class="c1"># define convenience functions</span>
</span><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><span class="line">
</span><span class="line"><span class="c1"># initialize parameters</span>
</span><span class="line">nsite <span class="o"><-</span> <span class="m">100</span>
</span><span class="line">nyear <span class="o"><-</span> <span class="m">5</span>
</span><span class="line">nrep <span class="o"><-</span> <span class="m">3</span>
</span><span class="line">
</span><span class="line"><span class="c1"># yearly Pr(colonization)</span>
</span><span class="line">p_beta <span class="o">=</span> <span class="m">0.5</span>
</span><span class="line">beta <span class="o"><-</span> logit<span class="p">(</span>p_beta<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># yearly Pr(persistence)</span>
</span><span class="line">p_rho <span class="o">=</span> <span class="m">0.6</span>
</span><span class="line">rho <span class="o"><-</span> logit<span class="p">(</span>p_rho<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># assume there are 2 site level covariates: covariate and treatment</span>
</span><span class="line"><span class="c1"># response to covariate</span>
</span><span class="line">p_alpha1 <span class="o"><-</span> <span class="m">.6</span>
</span><span class="line">alpha1 <span class="o"><-</span> logit<span class="p">(</span>p_alpha1<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># response to treatment </span>
</span><span class="line">p_alpha2 <span class="o"><-</span> <span class="m">.5</span>
</span><span class="line">alpha2 <span class="o"><-</span> logit<span class="p">(</span>p_alpha2<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># interaction between covariate & treatment</span>
</span><span class="line">p_alpha3 <span class="o"><-</span> <span class="m">.2</span>
</span><span class="line">alpha3 <span class="o"><-</span> logit<span class="p">(</span>p_alpha3<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Site-specific covariates</span>
</span><span class="line">covariate <span class="o"><-</span> rnorm<span class="p">(</span>nsite<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">treat <span class="o"><-</span> rbinom<span class="p">(</span>nsite<span class="p">,</span> size<span class="o">=</span><span class="m">1</span><span class="p">,</span> prob<span class="o">=</span><span class="m">.5</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial occupancy states</span>
</span><span class="line">z0 <span class="o"><-</span> rbinom<span class="p">(</span>nsite<span class="p">,</span> <span class="m">1</span><span class="p">,</span> <span class="m">.5</span><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> 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> 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> 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>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">      <span class="c1"># logit psi</span>
</span><span class="line">      lpsi<span class="p">[</span>j<span class="p">,</span> t<span class="p">]</span> <span class="o"><-</span> beta <span class="o">+</span> rho <span class="o">*</span> z0<span class="p">[</span>j<span class="p">]</span> <span class="o">+</span>
</span><span class="line">          covariate<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> alpha1 <span class="o">+</span> treat<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> alpha2 <span class="o">+</span>
</span><span class="line">          alpha3 <span class="o">*</span> covariate<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> treat<span class="p">[</span>j<span class="p">]</span>
</span><span class="line">      psi<span class="p">[</span>j<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> t<span class="p">])</span>
</span><span class="line">      z<span class="p">[</span>j<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> t<span class="p">])</span> <span class="c1"># true occupancy state</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> t<span class="p">]</span> <span class="o"><-</span> beta <span class="o">+</span> rho <span class="o">*</span> z<span class="p">[</span>j<span class="p">,</span> t <span class="o">-</span> <span class="m">1</span><span class="p">]</span> <span class="o">+</span>
</span><span class="line">        covariate<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> alpha1 <span class="o">+</span> treat<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> alpha2 <span class="o">+</span>
</span><span class="line">        alpha3 <span class="o">*</span> covariate<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> treat<span class="p">[</span>j<span class="p">]</span>
</span><span class="line">      psi<span class="p">[</span>j<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> t<span class="p">])</span>
</span><span class="line">      z<span class="p">[</span>j<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> 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><span class="line"><span class="c1"># detection probability</span>
</span><span class="line">p <span class="o"><-</span> <span class="m">0.6</span>
</span><span class="line">
</span><span class="line"><span class="c1"># observations</span>
</span><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> 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>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> 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="o">*</span> z<span class="p">[</span>j<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><span class="line"><span class="c1"># Subset the data to exclude the final year, which we'll use for model validation</span>
</span><span class="line">xsub <span class="o"><-</span> x
</span><span class="line">xsub<span class="p">[,</span> nyear<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> <span class="kc">NA</span>
</span>

For illustration, I included a strong interaction between treatment
and the continuous site level covariate (could be elevation, area, etc). As such, a
measure of model fit such as AUC ought to identify a saturated model as
the best fitting. Handily, AUC is a derived parameter, and common
occupancy model parameters can be used to estimate a posterior.
To generate a posterior AUC, we need predicted occupancy probabilities
($\psi$) and realized occupancy states ($Z$) in the final year. Predicted occupancy
probabilities can be produced using data from previous years, and
realized occupancy states are assumed to be represented by the posterior
for $Z$ generated from a single-year model, fit to the data from the final
year of the study.

Fitting a saturated model

Begin by modeling occupancy probabilities as a function of both covariates
and their interaction, predicting $\psi$ in the final year:

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
<span class="line"><span class="c1"># model specification</span>
</span><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">    p_beta ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    beta <- log(p_beta / (1 - p_beta))</span>
</span><span class="line"><span class="s">    p_rho ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    rho <- log(p_rho / (1 - p_rho))</span>
</span><span class="line"><span class="s">    p_rho0 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    rho0 <- log(p_rho0 / (1 - p_rho0))</span>
</span><span class="line"><span class="s">    p ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa1 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa2 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa3~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    a1 <- log(pa1 / (1 - pa1))</span>
</span><span class="line"><span class="s">    a2 <- log(pa2 / (1 - pa2))</span>
</span><span class="line"><span class="s">    a3 <- log(pa3 / (1 - pa3))</span>
</span><span class="line"><span class="s">    </span>
</span><span class="line"><span class="s">    #### occupancy model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite) {</span>
</span><span class="line"><span class="s">      z0[j] ~ dbern(rho0)</span>
</span><span class="line"><span class="s">      logit(psi[j, 1]) <- beta + rho * z0[j] + </span>
</span><span class="line"><span class="s">        a1 * covariate[j] + a2 * treat[j] + </span>
</span><span class="line"><span class="s">        a3 * covariate[j] * treat[j]</span>
</span><span class="line"><span class="s">      z[j, 1] ~ dbern(psi[j, 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, t]) <- beta + rho * z[j, t-1] + </span>
</span><span class="line"><span class="s">          a1 * covariate[j] + a2 * treat[j] + </span>
</span><span class="line"><span class="s">          a3 * covariate[j] * treat[j]</span>
</span><span class="line"><span class="s">        z[j, t] ~ dbern(psi[j, 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">    #### observation model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite){</span>
</span><span class="line"><span class="s">      for (t in 1:nyear){</span>
</span><span class="line"><span class="s">        mu[j, t] <- z[j, t] * p</span>
</span><span class="line"><span class="s">        for (k in 1:nrep){</span>
</span><span class="line"><span class="s">          x[j, t, k] ~ dbern(mu[j, 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 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">"saturated.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># bundle data</span>
</span><span class="line">data <span class="o"><-</span> list<span class="p">(</span>x <span class="o">=</span> xsub<span class="p">,</span> nrep <span class="o">=</span> nrep<span class="p">,</span> nsite <span class="o">=</span> nsite<span class="p">,</span> nyear <span class="o">=</span> nyear<span class="p">,</span> covariate<span class="o">=</span>covariate<span class="p">,</span> treat<span class="o">=</span>treat<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial values</span>
</span><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> 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>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> t<span class="p">]</span> <span class="o"><-</span> max<span class="p">(</span>xsub<span class="p">[</span>j<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><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>z <span class="o">=</span> zinit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># parameters to monitor</span>
</span><span class="line">params <span class="o"><-</span> c<span class="p">(</span><span class="s">"p"</span><span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> <span class="s">"rho"</span><span class="p">,</span> <span class="s">"psi"</span><span class="p">,</span> <span class="s">"pa1"</span><span class="p">,</span> <span class="s">"pa2"</span><span class="p">,</span> <span class="s">"pa3"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># use JAGS</span>
</span><span class="line">require<span class="p">(</span>rjags<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># build model</span>
</span><span class="line">nchain <span class="o"><-</span> <span class="m">3</span>
</span><span class="line">ocmod <span class="o"><-</span> jags.model<span class="p">(</span>file <span class="o">=</span> <span class="s">"saturated.txt"</span><span class="p">,</span> inits <span class="o">=</span> inits<span class="p">,</span>
</span><span class="line">                    data <span class="o">=</span> data<span class="p">,</span> n.chains <span class="o">=</span> nchain<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># specify MCMC settings and start sampling</span>
</span><span class="line">nburn <span class="o"><-</span> <span class="m">2000</span>
</span><span class="line">niter <span class="o"><-</span> <span class="m">2000</span>
</span><span class="line">nthin <span class="o"><-</span> <span class="m">2</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> niter<span class="p">,</span>
</span><span class="line">                    thin<span class="o">=</span>nthin<span class="p">,</span> variable.names <span class="o">=</span> params<span class="p">)</span>
</span>

Now that we have our posteriors for $\psi$ at each site in the final
year, we can fit a single-year model to the final year’s data to
estimate $Z$.

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
<span class="line">xlast <span class="o"><-</span> x<span class="p">[,</span>nyear<span class="p">,]</span>
</span><span class="line">
</span><span class="line"><span class="c1"># model specification</span>
</span><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">    p_beta ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    beta <- log(p_beta / (1 - p_beta))</span>
</span><span class="line"><span class="s">    p ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa1 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa2 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa3~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    a1 <- log(pa1 / (1 - pa1))</span>
</span><span class="line"><span class="s">    a2 <- log(pa2 / (1 - pa2))</span>
</span><span class="line"><span class="s">    a3 <- log(pa3 / (1 - pa3))</span>
</span><span class="line"><span class="s">    </span>
</span><span class="line"><span class="s">    #### occupancy model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite) {</span>
</span><span class="line"><span class="s">      logit(psi[j]) <- beta + </span>
</span><span class="line"><span class="s">        a1 * covariate[j] + a2 * treat[j] + </span>
</span><span class="line"><span class="s">        a3 * covariate[j] * treat[j]</span>
</span><span class="line"><span class="s">      z[j] ~ dbern(psi[j])</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">      mu[j] <- z[j] * p</span>
</span><span class="line"><span class="s">      for (k in 1:nrep){</span>
</span><span class="line"><span class="s">        x[j, k] ~ dbern(mu[j])</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">"saturated2.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># bundle data</span>
</span><span class="line">data2 <span class="o"><-</span> list<span class="p">(</span>x <span class="o">=</span> xlast<span class="p">,</span> nrep <span class="o">=</span> nrep<span class="p">,</span>
</span><span class="line">              nsite <span class="o">=</span> nsite<span class="p">,</span> covariate<span class="o">=</span>covariate<span class="p">,</span>
</span><span class="line">              treat<span class="o">=</span>treat<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial values</span>
</span><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>
</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">    zinit<span class="p">[</span>j<span class="p">]</span> <span class="o"><-</span> max<span class="p">(</span>xlast<span class="p">[</span>j<span class="p">,</span> <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>z <span class="o">=</span> zinit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># parameters to monitor</span>
</span><span class="line">params2 <span class="o"><-</span> c<span class="p">(</span><span class="s">"p"</span><span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> <span class="s">"pa1"</span><span class="p">,</span> <span class="s">"pa2"</span><span class="p">,</span> <span class="s">"pa3"</span><span class="p">,</span> <span class="s">"z"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># build model</span>
</span><span class="line">ocmod2 <span class="o"><-</span> jags.model<span class="p">(</span>file <span class="o">=</span> <span class="s">"saturated2.txt"</span><span class="p">,</span> inits <span class="o">=</span> inits<span class="p">,</span>
</span><span class="line">                    data <span class="o">=</span> data2<span class="p">,</span> n.chains <span class="o">=</span> nchain<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># start sampling</span>
</span><span class="line">update<span class="p">(</span>ocmod2<span class="p">,</span> n.iter <span class="o">=</span> nburn<span class="p">)</span>
</span><span class="line">out2 <span class="o"><-</span> coda.samples<span class="p">(</span>ocmod2<span class="p">,</span> n.iter <span class="o">=</span> niter<span class="p">,</span>
</span><span class="line">                    thin<span class="o">=</span>nthin<span class="p">,</span> variable.names <span class="o">=</span> params2<span class="p">)</span>
</span>

To set up the data for AUC calculations, produce site by
iteration arrays for $\psi$ and $Z$:

1
2
3
4
5
6
7
8
9
10
11
12
<span class="line">require<span class="p">(</span>ggmcmc<span class="p">)</span>
</span><span class="line">ggd <span class="o"><-</span> ggs<span class="p">(</span>out<span class="p">)</span>
</span><span class="line">ggd2 <span class="o"><-</span> ggs<span class="p">(</span>out2<span class="p">)</span>
</span><span class="line">nstore <span class="o"><-</span> niter <span class="o">/</span> nthin
</span><span class="line">psi.pred <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nsite<span class="p">,</span> nstore<span class="o">*</span>nchain<span class="p">))</span>
</span><span class="line">Z.est <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nsite<span class="p">,</span> nstore<span class="o">*</span>nchain<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><span class="line">  indx <span class="o"><-</span> which<span class="p">(</span>ggd<span class="o">$</span>Parameter <span class="o">==</span> paste<span class="p">(</span><span class="s">"psi["</span><span class="p">,</span> j<span class="p">,</span> <span class="s">","</span><span class="p">,</span> nyear<span class="p">,</span> <span class="s">"]"</span><span class="p">,</span> sep<span class="o">=</span><span class="s">""</span><span class="p">))</span>
</span><span class="line">  psi.pred<span class="p">[</span>j<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> ggd<span class="o">$</span>value<span class="p">[</span>indx<span class="p">]</span>
</span><span class="line">  indx <span class="o"><-</span> which<span class="p">(</span>ggd2<span class="o">$</span>Parameter <span class="o">==</span> paste<span class="p">(</span><span class="s">"z["</span><span class="p">,</span> j<span class="p">,</span> <span class="s">"]"</span><span class="p">,</span> sep<span class="o">=</span><span class="s">""</span><span class="p">))</span>
</span><span class="line">  Z.est<span class="p">[</span>j<span class="p">,]</span> <span class="o"><-</span> ggd2<span class="o">$</span>value<span class="p">[</span>indx<span class="p">]</span>
</span><span class="line"><span class="p">}</span>
</span>

Now generate the posterior for AUC and store data on the true and false
positive rates to produce ROC curves.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
<span class="line">require<span class="p">(</span>ROCR<span class="p">)</span>
</span><span class="line">AUC1 <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> nstore <span class="o">*</span> nchain<span class="p">)</span>
</span><span class="line">fpr <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">,</span> nsite <span class="o">+</span> <span class="m">1</span><span class="p">))</span>
</span><span class="line">tpr <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">,</span> nsite <span class="o">+</span> <span class="m">1</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><span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">)){</span>
</span><span class="line">  psi.vals <span class="o"><-</span> psi.pred<span class="p">[,</span> i<span class="p">]</span>
</span><span class="line">  Z.vals <span class="o"><-</span> Z.est<span class="p">[,</span> i<span class="p">]</span>
</span><span class="line">  pred <span class="o"><-</span> prediction<span class="p">(</span>psi.vals<span class="p">,</span> factor<span class="p">(</span>Z.vals<span class="p">,</span> levels<span class="o">=</span>c<span class="p">(</span><span class="s">"0"</span><span class="p">,</span> <span class="s">"1"</span><span class="p">)))</span>
</span><span class="line">  perf <span class="o"><-</span> performance<span class="p">(</span>pred<span class="p">,</span> <span class="s">"auc"</span><span class="p">)</span>
</span><span class="line">  AUC1<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>y.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line">  perf <span class="o"><-</span> performance<span class="p">(</span>pred<span class="p">,</span> <span class="s">"tpr"</span><span class="p">,</span><span class="s">"fpr"</span><span class="p">)</span>
</span><span class="line">  fpr<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>x.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line">  tpr<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>y.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">require<span class="p">(</span>reshape2<span class="p">)</span>
</span><span class="line">fprm <span class="o"><-</span> melt<span class="p">(</span>fpr<span class="p">,</span> varnames<span class="o">=</span>c<span class="p">(</span><span class="s">"iter"</span><span class="p">,</span> <span class="s">"site"</span><span class="p">))</span>
</span><span class="line">tprm <span class="o"><-</span> melt<span class="p">(</span>tpr<span class="p">,</span> varnames <span class="o">=</span> c<span class="p">(</span><span class="s">"iter"</span><span class="p">,</span> <span class="s">"site"</span><span class="p">))</span>
</span><span class="line">ROC1 <span class="o"><-</span> data.frame<span class="p">(</span>fpr <span class="o">=</span> fprm<span class="o">$</span>value<span class="p">,</span>
</span><span class="line">                   tpr <span class="o">=</span> tprm<span class="o">$</span>value<span class="p">,</span>
</span><span class="line">                   iter <span class="o">=</span> rep<span class="p">(</span>fprm<span class="o">$</span>iter<span class="p">,</span> <span class="m">2</span><span class="p">))</span>
</span>

Fitting a simpler model

Having fitted a saturated model, we can now fit a simpler model that
includes only main effects:

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
<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">    p_beta ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    beta <- log(p_beta / (1 - p_beta))</span>
</span><span class="line"><span class="s">    p_rho ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    rho <- log(p_rho / (1 - p_rho))</span>
</span><span class="line"><span class="s">    p_rho0 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    rho0 <- log(p_rho0 / (1 - p_rho0))</span>
</span><span class="line"><span class="s">    p ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa1 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa2 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    a1 <- log(pa1 / (1 - pa1))</span>
</span><span class="line"><span class="s">    a2 <- log(pa2 / (1 - pa2))</span>
</span><span class="line"><span class="s">    </span>
</span><span class="line"><span class="s">    #### occupancy model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite) {</span>
</span><span class="line"><span class="s">      z0[j] ~ dbern(rho0)</span>
</span><span class="line"><span class="s">      logit(psi[j, 1]) <- beta + rho * z0[j] + </span>
</span><span class="line"><span class="s">        a1 * covariate[j] + a2 * treat[j] </span>
</span><span class="line"><span class="s">      z[j, 1] ~ dbern(psi[j, 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, t]) <- beta + rho * z[j, t-1] + </span>
</span><span class="line"><span class="s">          a1 * covariate[j] + a2 * treat[j] </span>
</span><span class="line"><span class="s">        z[j, t] ~ dbern(psi[j, 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">    #### observation model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite){</span>
</span><span class="line"><span class="s">      for (t in 1:nyear){</span>
</span><span class="line"><span class="s">        mu[j, t] <- z[j, t] * p</span>
</span><span class="line"><span class="s">        for (k in 1:nrep){</span>
</span><span class="line"><span class="s">          x[j, t, k] ~ dbern(mu[j, 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 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">"mainef.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">
</span><span class="line"><span class="c1"># bundle data</span>
</span><span class="line">data <span class="o"><-</span> list<span class="p">(</span>x <span class="o">=</span> xsub<span class="p">,</span> nrep <span class="o">=</span> nrep<span class="p">,</span> nsite <span class="o">=</span> nsite<span class="p">,</span> nyear <span class="o">=</span> nyear<span class="p">,</span> covariate<span class="o">=</span>covariate<span class="p">,</span> treat<span class="o">=</span>treat<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial values</span>
</span><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> 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>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> t<span class="p">]</span> <span class="o"><-</span> max<span class="p">(</span>xsub<span class="p">[</span>j<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><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>z <span class="o">=</span> zinit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># parameters to monitor</span>
</span><span class="line">params <span class="o"><-</span> c<span class="p">(</span><span class="s">"p"</span><span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> <span class="s">"rho"</span><span class="p">,</span> <span class="s">"psi"</span><span class="p">,</span> <span class="s">"pa1"</span><span class="p">,</span> <span class="s">"pa2"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># build model</span>
</span><span class="line">ocmod <span class="o"><-</span> jags.model<span class="p">(</span>file <span class="o">=</span> <span class="s">"mainef.txt"</span><span class="p">,</span> inits <span class="o">=</span> inits<span class="p">,</span>
</span><span class="line">                    data <span class="o">=</span> data<span class="p">,</span> n.chains <span class="o">=</span> nchain<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># specify MCMC settings and start sampling</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> niter<span class="p">,</span>
</span><span class="line">                    thin<span class="o">=</span>nthin<span class="p">,</span> variable.names <span class="o">=</span> params<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Fit a single-year model for the final year to estimate Z ##</span>
</span><span class="line"><span class="c1"># model specification</span>
</span><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">    p_beta ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    beta <- log(p_beta / (1 - p_beta))</span>
</span><span class="line"><span class="s">    p ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa1 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    pa2 ~ dbeta(1, 1)</span>
</span><span class="line"><span class="s">    a1 <- log(pa1 / (1 - pa1))</span>
</span><span class="line"><span class="s">    a2 <- log(pa2 / (1 - pa2))</span>
</span><span class="line"><span class="s">    </span>
</span><span class="line"><span class="s">    #### occupancy model</span>
</span><span class="line"><span class="s">    for (j in 1:nsite) {</span>
</span><span class="line"><span class="s">      logit(psi[j]) <- beta + </span>
</span><span class="line"><span class="s">        a1 * covariate[j] + a2 * treat[j]</span>
</span><span class="line"><span class="s">      z[j] ~ dbern(psi[j])</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">      mu[j] <- z[j] * p</span>
</span><span class="line"><span class="s">      for (k in 1:nrep){</span>
</span><span class="line"><span class="s">        x[j, k] ~ dbern(mu[j])</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">"mainef2.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># bundle data</span>
</span><span class="line">data2 <span class="o"><-</span> list<span class="p">(</span>x <span class="o">=</span> xlast<span class="p">,</span> nrep <span class="o">=</span> nrep<span class="p">,</span>
</span><span class="line">              nsite <span class="o">=</span> nsite<span class="p">,</span> covariate<span class="o">=</span>covariate<span class="p">,</span>
</span><span class="line">              treat<span class="o">=</span>treat<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># initial values</span>
</span><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>
</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">  zinit<span class="p">[</span>j<span class="p">]</span> <span class="o"><-</span> max<span class="p">(</span>xlast<span class="p">[</span>j<span class="p">,</span> <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>z <span class="o">=</span> zinit<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># parameters to monitor</span>
</span><span class="line">params2 <span class="o"><-</span> c<span class="p">(</span><span class="s">"p"</span><span class="p">,</span> <span class="s">"beta"</span><span class="p">,</span> <span class="s">"pa1"</span><span class="p">,</span> <span class="s">"pa2"</span><span class="p">,</span> <span class="s">"z"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># build model</span>
</span><span class="line">ocmod2 <span class="o"><-</span> jags.model<span class="p">(</span>file <span class="o">=</span> <span class="s">"mainef2.txt"</span><span class="p">,</span> inits <span class="o">=</span> inits<span class="p">,</span>
</span><span class="line">                     data <span class="o">=</span> data2<span class="p">,</span> n.chains <span class="o">=</span> nchain<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># start sampling</span>
</span><span class="line">update<span class="p">(</span>ocmod2<span class="p">,</span> n.iter <span class="o">=</span> nburn<span class="p">)</span>
</span><span class="line">out2 <span class="o"><-</span> coda.samples<span class="p">(</span>ocmod2<span class="p">,</span> n.iter <span class="o">=</span> niter<span class="p">,</span>
</span><span class="line">                     thin<span class="o">=</span>nthin<span class="p">,</span> variable.names <span class="o">=</span> params2<span class="p">)</span>
</span><span class="line">
</span><span class="line">
</span><span class="line"><span class="c1"># Generating AUC posterior...</span>
</span><span class="line">ggd <span class="o"><-</span> ggs<span class="p">(</span>out<span class="p">)</span>
</span><span class="line">ggd2 <span class="o"><-</span> ggs<span class="p">(</span>out2<span class="p">)</span>
</span><span class="line">nstore <span class="o"><-</span> niter <span class="o">/</span> nthin
</span><span class="line">psi.pred <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nsite<span class="p">,</span> nstore<span class="o">*</span>nchain<span class="p">))</span>
</span><span class="line">Z.est <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nsite<span class="p">,</span> nstore<span class="o">*</span>nchain<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><span class="line">  indx <span class="o"><-</span> which<span class="p">(</span>ggd<span class="o">$</span>Parameter <span class="o">==</span> paste<span class="p">(</span><span class="s">"psi["</span><span class="p">,</span> j<span class="p">,</span> <span class="s">","</span><span class="p">,</span> nyear<span class="p">,</span> <span class="s">"]"</span><span class="p">,</span> sep<span class="o">=</span><span class="s">""</span><span class="p">))</span>
</span><span class="line">  psi.pred<span class="p">[</span>j<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> ggd<span class="o">$</span>value<span class="p">[</span>indx<span class="p">]</span>
</span><span class="line">  indx <span class="o"><-</span> which<span class="p">(</span>ggd2<span class="o">$</span>Parameter <span class="o">==</span> paste<span class="p">(</span><span class="s">"z["</span><span class="p">,</span> j<span class="p">,</span> <span class="s">"]"</span><span class="p">,</span> sep<span class="o">=</span><span class="s">""</span><span class="p">))</span>
</span><span class="line">  Z.est<span class="p">[</span>j<span class="p">,]</span> <span class="o"><-</span> ggd2<span class="o">$</span>value<span class="p">[</span>indx<span class="p">]</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">AUC2 <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> nstore <span class="o">*</span> nchain<span class="p">)</span>
</span><span class="line">fpr <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">,</span> nsite <span class="o">+</span> <span class="m">1</span><span class="p">))</span>
</span><span class="line">tpr <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">,</span> nsite <span class="o">+</span> <span class="m">1</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><span class="p">(</span>nstore <span class="o">*</span> nchain<span class="p">)){</span>
</span><span class="line">  psi.vals <span class="o"><-</span> psi.pred<span class="p">[,</span> i<span class="p">]</span>
</span><span class="line">  Z.vals <span class="o"><-</span> Z.est<span class="p">[,</span> i<span class="p">]</span>
</span><span class="line">  pred <span class="o"><-</span> prediction<span class="p">(</span>psi.vals<span class="p">,</span> factor<span class="p">(</span>Z.vals<span class="p">,</span> levels<span class="o">=</span>c<span class="p">(</span><span class="s">"0"</span><span class="p">,</span> <span class="s">"1"</span><span class="p">)))</span>
</span><span class="line">  perf <span class="o"><-</span> performance<span class="p">(</span>pred<span class="p">,</span> <span class="s">"auc"</span><span class="p">)</span>
</span><span class="line">  AUC2<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>y.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line">  perf <span class="o"><-</span> performance<span class="p">(</span>pred<span class="p">,</span> <span class="s">"tpr"</span><span class="p">,</span><span class="s">"fpr"</span><span class="p">)</span>
</span><span class="line">  fpr<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>x.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line">  tpr<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> perf<span class="o">@</span>y.values<span class="p">[[</span><span class="m">1</span><span class="p">]]</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">fprm <span class="o"><-</span> melt<span class="p">(</span>fpr<span class="p">,</span> varnames<span class="o">=</span>c<span class="p">(</span><span class="s">"iter"</span><span class="p">,</span> <span class="s">"site"</span><span class="p">))</span>
</span><span class="line">tprm <span class="o"><-</span> melt<span class="p">(</span>tpr<span class="p">,</span> varnames <span class="o">=</span> c<span class="p">(</span><span class="s">"iter"</span><span class="p">,</span> <span class="s">"site"</span><span class="p">))</span>
</span><span class="line">ROC2 <span class="o"><-</span> data.frame<span class="p">(</span>fpr <span class="o">=</span> fprm<span class="o">$</span>value<span class="p">,</span>
</span><span class="line">                   tpr <span class="o">=</span> tprm<span class="o">$</span>value<span class="p">,</span>
</span><span class="line">                   iter <span class="o">=</span> rep<span class="p">(</span>fprm<span class="o">$</span>iter<span class="p">,</span> <span class="m">2</span><span class="p">))</span>
</span>

Comparing models

How well did our models predict occupancy in the final year of the study,
and was one better than the other?
We can answer this question by inspecting posteriors for AUC (larger
values are better), and the ROC curves.

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
<span class="line">p1 <span class="o"><-</span> ggplot<span class="p">(</span>ROC1<span class="p">,</span> aes<span class="p">(</span>x<span class="o">=</span>fpr<span class="p">,</span> y<span class="o">=</span>tpr<span class="p">,</span> group<span class="o">=</span>iter<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  geom_line<span class="p">(</span>alpha<span class="o">=</span><span class="m">0.05</span><span class="p">,</span> color<span class="o">=</span><span class="s">"blue"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_abline<span class="p">(</span>intercept <span class="o">=</span> <span class="m">0</span><span class="p">,</span> slope <span class="o">=</span> <span class="m">1</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">  ylab<span class="p">(</span><span class="s">"True positive rate"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"False positive rate"</span><span class="p">)</span> <span class="o">+</span> theme_bw<span class="p">()</span>
</span><span class="line">
</span><span class="line">p2 <span class="o"><-</span> ggplot<span class="p">(</span>ROC2<span class="p">,</span> aes<span class="p">(</span>x<span class="o">=</span>fpr<span class="p">,</span> y<span class="o">=</span>tpr<span class="p">,</span> group<span class="o">=</span>iter<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  geom_line<span class="p">(</span>alpha<span class="o">=</span><span class="m">0.05</span><span class="p">,</span> color<span class="o">=</span><span class="s">"red"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_abline<span class="p">(</span>intercept <span class="o">=</span> <span class="m">0</span><span class="p">,</span> slope <span class="o">=</span> <span class="m">1</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">  ylab<span class="p">(</span><span class="s">"True positive rate"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"False positive rate"</span><span class="p">)</span> <span class="o">+</span> theme_bw<span class="p">()</span>
</span><span class="line">
</span><span class="line">p3 <span class="o"><-</span> ggplot<span class="p">(</span>data.frame<span class="p">(</span>AUC <span class="o">=</span> c<span class="p">(</span>AUC1<span class="p">,</span> AUC2<span class="p">),</span>
</span><span class="line">              Model <span class="o">=</span> rep<span class="p">(</span>c<span class="p">(</span><span class="s">"Saturated model"</span><span class="p">,</span> <span class="s">"Main effects model "</span><span class="p">),</span>
</span><span class="line">                          each<span class="o">=</span>length<span class="p">(</span>AUC1<span class="p">))),</span>
</span><span class="line">       aes<span class="p">(</span>x<span class="o">=</span>AUC<span class="p">,</span> fill<span class="o">=</span>Model<span class="p">,</span> color<span class="o">=</span>Model<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  geom_density<span class="p">(</span>alpha<span class="o">=</span><span class="m">.8</span><span class="p">)</span> <span class="o">+</span> theme_bw<span class="p">()</span> <span class="o">+</span> ylab<span class="p">(</span><span class="s">"Posterior density"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"Area under receiver operating characteristic curve (AUC)"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  scale_colour_manual<span class="p">(</span>values<span class="o">=</span>c<span class="p">(</span><span class="s">"red"</span><span class="p">,</span> <span class="s">"blue"</span><span class="p">))</span> <span class="o">+</span>
</span><span class="line">  scale_fill_manual<span class="p">(</span>values<span class="o">=</span>c<span class="p">(</span><span class="s">"red"</span><span class="p">,</span> <span class="s">"blue"</span><span class="p">))</span> <span class="o">+</span>
</span><span class="line">  theme<span class="p">(</span>legend.position <span class="o">=</span> <span class="s">"top"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  theme<span class="p">(</span>legend.title<span class="o">=</span>element_blank<span class="p">())</span> <span class="o">+</span>
</span><span class="line">  ggtitle<span class="p">(</span><span class="s">"Posterior AUC and ROC curve comparison"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">require<span class="p">(</span>gridExtra<span class="p">)</span>
</span><span class="line">grid.arrange<span class="p">(</span>p3<span class="p">,</span> arrangeGrob<span class="p">(</span>p2<span class="p">,</span> p1<span class="p">,</span> ncol<span class="o">=</span><span class="m">2</span><span class="p">),</span> ncol<span class="o">=</span><span class="m">1</span><span class="p">)</span>
</span>

As expected, the model that generated the data fits better than the
model that excludes the strong interaction term. Note that AUC reflects
the accuracy of model predictions, and does not penalize model
complexity.

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)