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
|
# ubms: Unmarked Bayesian Models with Stan
<!-- badges: start -->
[](https://github.com/kenkellner/ubms/actions)
[](https://codecov.io/gh/kenkellner/ubms?branch=master)
[](https://cran.r-project.org/web/packages/ubms/index.html)
<!-- badges: end -->
`ubms` is an `R` package for fitting Bayesian hierarchical models of
animal occurrence and abundance. The package has a formula-based
interface compatible with
[unmarked](https://cran.r-project.org/web/packages/unmarked/index.html),
but the model is fit using MCMC with [Stan](https://mc-stan.org/)
instead of using maximum likelihood. Currently there are Stan versions
of unmarked functions `occu`, `occuRN`, `colext`, `occuTTD`, `pcount`,
`distsamp`, and `multinomPois`. These functions follow the `stan_`
prefix naming format established by
[rstanarm](https://cran.r-project.org/web/packages/rstanarm/index.html).
For example, the Stan version of the `unmarked` function `occu` is
`stan_occu`.
Advantages compared to `unmarked`:
1. Obtain posterior distributions of parameters and derived parameters
2. Include random effects in parameter formulas (same syntax as `lme4`)
3. Assess model fit using WAIC and LOO via the
[loo](https://cran.r-project.org/web/packages/loo/index.html)
package
Disadvantages compared to `unmarked`:
1. MCMC is slower than maximum likelihood
2. Not all model types are supported
3. Potential for convergence issues
## Installation
`ubms` is on
[CRAN](https://cran.r-project.org/web/packages/ubms/index.html):
``` r
install.packages("ubms")
```
Alternatively, the latest development version can be installed from
Github:
``` r
# install.packages("devtools")
devtools::install_github("kenkellner/ubms")
```
## Example
Simulate occupancy data including a random effect on occupancy:
``` r
library(ubms)
set.seed(123)
dat_occ <- data.frame(x1=rnorm(500))
dat_p <- data.frame(x2=rnorm(500*5))
y <- matrix(NA, 500, 5)
z <- rep(NA, 500)
b <- c(0.4, -0.5, 0.3, 0.5)
re_fac <- factor(sample(letters[1:26], 500, replace=T))
dat_occ$group <- re_fac
re <- rnorm(26, 0, 1.2)
re_idx <- as.numeric(re_fac)
idx <- 1
for (i in 1:500){
z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
for (j in 1:5){
y[i,j] <- z[i]*rbinom(1,1,
plogis(b[3] + b[4]*dat_p$x2[idx]))
idx <- idx + 1
}
}
```
Create `unmarked` frame:
``` r
umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
```
Fit a model with a random intercept, using 3 parallel cores:
``` r
(fm <- stan_occu(~x2 ~x1 + (1|group), umf, chains=3, cores=3, seed=123))
```
##
## Call:
## stan_occu(formula = ~x2 ~ x1 + (1 | group), data = umf, chains = 3,
## cores = 3, refresh = 0, seed = 123)
##
## Occupancy (logit-scale):
## Estimate SD 2.5% 97.5% n_eff Rhat
## (Intercept) 0.318 0.290 -0.254 0.887 926 1.001
## x1 -0.461 0.114 -0.683 -0.237 4428 0.999
## sigma [1|group] 1.337 0.261 0.911 1.924 2112 1.002
##
## Detection (logit-scale):
## Estimate SD 2.5% 97.5% n_eff Rhat
## (Intercept) 0.383 0.0607 0.265 0.503 4539 1
## x2 0.586 0.0591 0.470 0.704 4963 1
##
## LOOIC: 2267.291
## Runtime: 27.948 sec
Examine residuals for occupancy and detection submodels (following
[Wright et al. 2019](https://doi.org/10.1002/ecy.2703)). Each panel
represents a draw from the posterior distribution.
``` r
plot(fm)
```
<!-- -->
Assess model goodness-of-fit with a posterior predictive check, using
the MacKenzie-Bailey chi-square test:
``` r
fm_fit <- gof(fm, draws=500)
plot(fm_fit)
```
<!-- -->
Look at the marginal effect of `x2` on detection:
``` r
plot_effects(fm, "det")
```
<!-- -->
A more detailed vignette for the package is available
[here](https://kenkellner.com/blog/ubms-vignette.html).
|
