Class 5 live demo: Survival estimation
The team
last updated: 2021-05-15
Introduction
We are going to analyze real capture-recapture data. The data concern the European Dipper (Cinclus cinclus). Captures were carried out for 7 years (1981-1987) in eastern France by Gilbert Marzolin who kindly provided us with the data. The data consist of initial markings and recaptures of almost 300 breeding adults each year during the March-June period. Birds were at least 1 year old when initially banded.
Load nimble first.
library(nimble)Data
Read in the data.
dipper <- read_csv("dipper.csv")
dipper %>%
kableExtra::kable() %>%
kableExtra::scroll_box(width = "100%", height = "400px")| year_1981 | year_1982 | year_1983 | year_1984 | year_1985 | year_1986 | year_1987 | sex | wing_length |
|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 1 | 0 | M | 95 |
| 1 | 1 | 1 | 1 | 1 | 0 | 0 | F | 88 |
| 1 | 1 | 1 | 1 | 0 | 0 | 0 | M | 94 |
| 1 | 1 | 1 | 1 | 0 | 0 | 0 | F | 85 |
| 1 | 1 | 0 | 1 | 1 | 1 | 0 | F | 86 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | M | 97 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | M | 98 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | F | 89 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | F | 86 |
| 1 | 0 | 1 | 0 | 0 | 0 | 0 | M | 98 |
| 1 | 0 | 1 | 0 | 0 | 0 | 0 | F | 92 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | M | 97 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | M | 95 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | M | 98 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | F | 91 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | F | 89 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | F | 87 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | F | 90 |
| 0 | 1 | 1 | 1 | 1 | 1 | 1 | F | 87 |
| 0 | 1 | 1 | 1 | 1 | 1 | 1 | F | 86 |
| 0 | 1 | 1 | 1 | 1 | 1 | 0 | F | 88 |
| 0 | 1 | 1 | 1 | 1 | 0 | 0 | M | 99 |
| 0 | 1 | 1 | 1 | 1 | 0 | 0 | F | 84 |
| 0 | 1 | 1 | 1 | 1 | 0 | 0 | F | 87 |
| 0 | 1 | 1 | 1 | 0 | 0 | 0 | M | 96 |
| 0 | 1 | 1 | 1 | 0 | 0 | 0 | F | 89 |
| 0 | 1 | 1 | 0 | 1 | 1 | 0 | F | 89 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 91 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 98 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 94 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 90 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | M | 93 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | F | 91 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | F | 84 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 1 | 1 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 94 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 93 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 92 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 94 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 98 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 96 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | M | 98 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 93 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 85 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 84 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 87 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 1 | 1 | 1 | 1 | 1 | F | 89 |
| 0 | 0 | 1 | 1 | 1 | 1 | 1 | F | 87 |
| 0 | 0 | 1 | 1 | 1 | 1 | 0 | M | 96 |
| 0 | 0 | 1 | 1 | 1 | 1 | 0 | F | 85 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | M | 98 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | M | 99 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | F | 90 |
| 0 | 0 | 1 | 1 | 1 | 0 | 0 | F | 86 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 96 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 98 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 93 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 94 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 98 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | M | 96 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | F | 90 |
| 0 | 0 | 1 | 1 | 0 | 0 | 0 | F | 87 |
| 0 | 0 | 1 | 0 | 1 | 1 | 0 | M | 95 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 93 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 99 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 96 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 94 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 95 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 100 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 92 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 96 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 91 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 90 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 91 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 87 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 90 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 84 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 87 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 85 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 90 |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | F | 86 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 98 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 100 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 94 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 97 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | M | 95 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | F | 85 |
| 0 | 0 | 0 | 1 | 1 | 1 | 1 | F | 85 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | M | 93 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | M | 99 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | M | 95 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | F | 86 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | F | 85 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | F | 92 |
| 0 | 0 | 0 | 1 | 1 | 1 | 0 | F | 89 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 97 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 92 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 96 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 95 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | M | 95 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | F | 87 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | F | 87 |
| 0 | 0 | 0 | 1 | 1 | 0 | 0 | F | 87 |
| 0 | 0 | 0 | 1 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 1 | 0 | 0 | 1 | M | 96 |
| 0 | 0 | 0 | 1 | 0 | 0 | 1 | F | 87 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 94 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 92 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 96 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 95 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 97 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | M | 95 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 92 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 89 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 87 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 100 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | M | 98 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 85 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 89 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 1 | 1 | 1 | F | 91 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | M | 96 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | M | 97 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | M | 95 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 86 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 87 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 90 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 90 |
| 0 | 0 | 0 | 0 | 1 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 96 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 98 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 98 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 97 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 97 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 98 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | M | 94 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 89 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 89 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 92 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 91 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 85 |
| 0 | 0 | 0 | 0 | 1 | 0 | 0 | F | 85 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 99 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 92 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 90 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 93 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 86 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 1 | 1 | F | 91 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 98 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 99 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 98 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 91 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 99 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | M | 91 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 91 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 84 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 86 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 89 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 87 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 92 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 98 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 94 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 96 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 97 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | M | 95 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 86 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 87 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 91 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 86 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 92 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 86 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 89 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 88 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 90 |
| 0 | 0 | 0 | 0 | 0 | 0 | 1 | F | 93 |
Format the data.
y <- dipper %>%
select(year_1981:year_1987) %>%
as.matrix()
head(y) year_1981 year_1982 year_1983 year_1984 year_1985 year_1986 year_1987
[1,] 1 1 1 1 1 1 0
[2,] 1 1 1 1 1 0 0
[3,] 1 1 1 1 0 0 0
[4,] 1 1 1 1 0 0 0
[5,] 1 1 0 1 1 1 0
[6,] 1 1 0 0 0 0 0CJS models
We start the session by fitting models with or without a time effect on survival and recapture probabilities.
A model with constant parameters.
hmm.phip <- nimbleCode({
phi ~ dunif(0, 1) # prior survival
p ~ dunif(0, 1) # prior detection
# likelihood
gamma[1,1] <- phi # Pr(alive t -> alive t+1)
gamma[1,2] <- 1 - phi # Pr(alive t -> dead t+1)
gamma[2,1] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2] <- 1 # Pr(dead t -> dead t+1)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})Get the occasion of first capture for all individuals.
first <- apply(y, 1, function(x) min(which(x !=0)))A list with constants.
my.constants <- list(N = nrow(y),
T = ncol(y),
first = first)
my.constants$N
[1] 294
$T
[1] 7
$first
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3
[75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[112] 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[149] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[186] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6
[223] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7
[260] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7Now the data in a list. Note that we add 1 to the data to have 1 for non-detections and 2 for detections. You may use the coding you prefer of course, you will just need to adjust the \(\Omega\) and \(\Gamma\) matrices in the model above.
my.data <- list(y = y + 1)Specify initial values. For the latent states, we go for the easy way, and say that all individuals are alive through the study period.
zinits <- y + 1 # non-detection -> alive
zinits[zinits == 2] <- 1 # dead -> alive
initial.values <- function() list(phi = runif(1,0,1),
p = runif(1,0,1),
z = zinits)
initial.values()$phi
[1] 0.03299212
$p
[1] 0.1289726
$z
year_1981 year_1982 year_1983 year_1984 year_1985 year_1986 year_1987
[1,] 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1
[6,] 1 1 1 1 1 1 1
[7,] 1 1 1 1 1 1 1
[8,] 1 1 1 1 1 1 1
[9,] 1 1 1 1 1 1 1
[10,] 1 1 1 1 1 1 1
[11,] 1 1 1 1 1 1 1
[12,] 1 1 1 1 1 1 1
[13,] 1 1 1 1 1 1 1
[14,] 1 1 1 1 1 1 1
[15,] 1 1 1 1 1 1 1
[16,] 1 1 1 1 1 1 1
[17,] 1 1 1 1 1 1 1
[18,] 1 1 1 1 1 1 1
[19,] 1 1 1 1 1 1 1
[20,] 1 1 1 1 1 1 1
[21,] 1 1 1 1 1 1 1
[22,] 1 1 1 1 1 1 1
[23,] 1 1 1 1 1 1 1
[24,] 1 1 1 1 1 1 1
[25,] 1 1 1 1 1 1 1
[26,] 1 1 1 1 1 1 1
[27,] 1 1 1 1 1 1 1
[28,] 1 1 1 1 1 1 1
[29,] 1 1 1 1 1 1 1
[30,] 1 1 1 1 1 1 1
[31,] 1 1 1 1 1 1 1
[32,] 1 1 1 1 1 1 1
[33,] 1 1 1 1 1 1 1
[34,] 1 1 1 1 1 1 1
[35,] 1 1 1 1 1 1 1
[36,] 1 1 1 1 1 1 1
[37,] 1 1 1 1 1 1 1
[38,] 1 1 1 1 1 1 1
[39,] 1 1 1 1 1 1 1
[40,] 1 1 1 1 1 1 1
[41,] 1 1 1 1 1 1 1
[42,] 1 1 1 1 1 1 1
[43,] 1 1 1 1 1 1 1
[44,] 1 1 1 1 1 1 1
[45,] 1 1 1 1 1 1 1
[46,] 1 1 1 1 1 1 1
[47,] 1 1 1 1 1 1 1
[48,] 1 1 1 1 1 1 1
[49,] 1 1 1 1 1 1 1
[50,] 1 1 1 1 1 1 1
[51,] 1 1 1 1 1 1 1
[52,] 1 1 1 1 1 1 1
[53,] 1 1 1 1 1 1 1
[54,] 1 1 1 1 1 1 1
[55,] 1 1 1 1 1 1 1
[56,] 1 1 1 1 1 1 1
[57,] 1 1 1 1 1 1 1
[58,] 1 1 1 1 1 1 1
[59,] 1 1 1 1 1 1 1
[60,] 1 1 1 1 1 1 1
[61,] 1 1 1 1 1 1 1
[62,] 1 1 1 1 1 1 1
[63,] 1 1 1 1 1 1 1
[64,] 1 1 1 1 1 1 1
[65,] 1 1 1 1 1 1 1
[66,] 1 1 1 1 1 1 1
[67,] 1 1 1 1 1 1 1
[68,] 1 1 1 1 1 1 1
[69,] 1 1 1 1 1 1 1
[70,] 1 1 1 1 1 1 1
[71,] 1 1 1 1 1 1 1
[72,] 1 1 1 1 1 1 1
[73,] 1 1 1 1 1 1 1
[74,] 1 1 1 1 1 1 1
[75,] 1 1 1 1 1 1 1
[76,] 1 1 1 1 1 1 1
[77,] 1 1 1 1 1 1 1
[78,] 1 1 1 1 1 1 1
[79,] 1 1 1 1 1 1 1
[80,] 1 1 1 1 1 1 1
[81,] 1 1 1 1 1 1 1
[82,] 1 1 1 1 1 1 1
[83,] 1 1 1 1 1 1 1
[84,] 1 1 1 1 1 1 1
[85,] 1 1 1 1 1 1 1
[86,] 1 1 1 1 1 1 1
[87,] 1 1 1 1 1 1 1
[88,] 1 1 1 1 1 1 1
[89,] 1 1 1 1 1 1 1
[90,] 1 1 1 1 1 1 1
[91,] 1 1 1 1 1 1 1
[92,] 1 1 1 1 1 1 1
[93,] 1 1 1 1 1 1 1
[94,] 1 1 1 1 1 1 1
[95,] 1 1 1 1 1 1 1
[96,] 1 1 1 1 1 1 1
[97,] 1 1 1 1 1 1 1
[98,] 1 1 1 1 1 1 1
[99,] 1 1 1 1 1 1 1
[100,] 1 1 1 1 1 1 1
[101,] 1 1 1 1 1 1 1
[102,] 1 1 1 1 1 1 1
[103,] 1 1 1 1 1 1 1
[104,] 1 1 1 1 1 1 1
[105,] 1 1 1 1 1 1 1
[106,] 1 1 1 1 1 1 1
[107,] 1 1 1 1 1 1 1
[108,] 1 1 1 1 1 1 1
[109,] 1 1 1 1 1 1 1
[110,] 1 1 1 1 1 1 1
[111,] 1 1 1 1 1 1 1
[112,] 1 1 1 1 1 1 1
[113,] 1 1 1 1 1 1 1
[114,] 1 1 1 1 1 1 1
[115,] 1 1 1 1 1 1 1
[116,] 1 1 1 1 1 1 1
[117,] 1 1 1 1 1 1 1
[118,] 1 1 1 1 1 1 1
[119,] 1 1 1 1 1 1 1
[120,] 1 1 1 1 1 1 1
[121,] 1 1 1 1 1 1 1
[122,] 1 1 1 1 1 1 1
[123,] 1 1 1 1 1 1 1
[124,] 1 1 1 1 1 1 1
[125,] 1 1 1 1 1 1 1
[126,] 1 1 1 1 1 1 1
[127,] 1 1 1 1 1 1 1
[128,] 1 1 1 1 1 1 1
[129,] 1 1 1 1 1 1 1
[130,] 1 1 1 1 1 1 1
[131,] 1 1 1 1 1 1 1
[132,] 1 1 1 1 1 1 1
[133,] 1 1 1 1 1 1 1
[134,] 1 1 1 1 1 1 1
[135,] 1 1 1 1 1 1 1
[136,] 1 1 1 1 1 1 1
[137,] 1 1 1 1 1 1 1
[138,] 1 1 1 1 1 1 1
[139,] 1 1 1 1 1 1 1
[140,] 1 1 1 1 1 1 1
[141,] 1 1 1 1 1 1 1
[142,] 1 1 1 1 1 1 1
[143,] 1 1 1 1 1 1 1
[144,] 1 1 1 1 1 1 1
[145,] 1 1 1 1 1 1 1
[146,] 1 1 1 1 1 1 1
[147,] 1 1 1 1 1 1 1
[148,] 1 1 1 1 1 1 1
[149,] 1 1 1 1 1 1 1
[150,] 1 1 1 1 1 1 1
[151,] 1 1 1 1 1 1 1
[152,] 1 1 1 1 1 1 1
[153,] 1 1 1 1 1 1 1
[154,] 1 1 1 1 1 1 1
[155,] 1 1 1 1 1 1 1
[156,] 1 1 1 1 1 1 1
[157,] 1 1 1 1 1 1 1
[158,] 1 1 1 1 1 1 1
[159,] 1 1 1 1 1 1 1
[160,] 1 1 1 1 1 1 1
[161,] 1 1 1 1 1 1 1
[162,] 1 1 1 1 1 1 1
[163,] 1 1 1 1 1 1 1
[164,] 1 1 1 1 1 1 1
[165,] 1 1 1 1 1 1 1
[166,] 1 1 1 1 1 1 1
[167,] 1 1 1 1 1 1 1
[168,] 1 1 1 1 1 1 1
[169,] 1 1 1 1 1 1 1
[170,] 1 1 1 1 1 1 1
[171,] 1 1 1 1 1 1 1
[172,] 1 1 1 1 1 1 1
[173,] 1 1 1 1 1 1 1
[174,] 1 1 1 1 1 1 1
[175,] 1 1 1 1 1 1 1
[176,] 1 1 1 1 1 1 1
[177,] 1 1 1 1 1 1 1
[178,] 1 1 1 1 1 1 1
[179,] 1 1 1 1 1 1 1
[180,] 1 1 1 1 1 1 1
[181,] 1 1 1 1 1 1 1
[182,] 1 1 1 1 1 1 1
[183,] 1 1 1 1 1 1 1
[184,] 1 1 1 1 1 1 1
[185,] 1 1 1 1 1 1 1
[186,] 1 1 1 1 1 1 1
[187,] 1 1 1 1 1 1 1
[188,] 1 1 1 1 1 1 1
[189,] 1 1 1 1 1 1 1
[190,] 1 1 1 1 1 1 1
[191,] 1 1 1 1 1 1 1
[192,] 1 1 1 1 1 1 1
[193,] 1 1 1 1 1 1 1
[194,] 1 1 1 1 1 1 1
[195,] 1 1 1 1 1 1 1
[196,] 1 1 1 1 1 1 1
[197,] 1 1 1 1 1 1 1
[198,] 1 1 1 1 1 1 1
[199,] 1 1 1 1 1 1 1
[200,] 1 1 1 1 1 1 1
[201,] 1 1 1 1 1 1 1
[202,] 1 1 1 1 1 1 1
[203,] 1 1 1 1 1 1 1
[204,] 1 1 1 1 1 1 1
[205,] 1 1 1 1 1 1 1
[206,] 1 1 1 1 1 1 1
[207,] 1 1 1 1 1 1 1
[208,] 1 1 1 1 1 1 1
[209,] 1 1 1 1 1 1 1
[210,] 1 1 1 1 1 1 1
[211,] 1 1 1 1 1 1 1
[212,] 1 1 1 1 1 1 1
[213,] 1 1 1 1 1 1 1
[214,] 1 1 1 1 1 1 1
[215,] 1 1 1 1 1 1 1
[216,] 1 1 1 1 1 1 1
[217,] 1 1 1 1 1 1 1
[218,] 1 1 1 1 1 1 1
[219,] 1 1 1 1 1 1 1
[220,] 1 1 1 1 1 1 1
[221,] 1 1 1 1 1 1 1
[222,] 1 1 1 1 1 1 1
[223,] 1 1 1 1 1 1 1
[224,] 1 1 1 1 1 1 1
[225,] 1 1 1 1 1 1 1
[226,] 1 1 1 1 1 1 1
[227,] 1 1 1 1 1 1 1
[228,] 1 1 1 1 1 1 1
[229,] 1 1 1 1 1 1 1
[230,] 1 1 1 1 1 1 1
[231,] 1 1 1 1 1 1 1
[232,] 1 1 1 1 1 1 1
[233,] 1 1 1 1 1 1 1
[234,] 1 1 1 1 1 1 1
[235,] 1 1 1 1 1 1 1
[236,] 1 1 1 1 1 1 1
[237,] 1 1 1 1 1 1 1
[238,] 1 1 1 1 1 1 1
[239,] 1 1 1 1 1 1 1
[240,] 1 1 1 1 1 1 1
[241,] 1 1 1 1 1 1 1
[242,] 1 1 1 1 1 1 1
[243,] 1 1 1 1 1 1 1
[244,] 1 1 1 1 1 1 1
[245,] 1 1 1 1 1 1 1
[246,] 1 1 1 1 1 1 1
[247,] 1 1 1 1 1 1 1
[248,] 1 1 1 1 1 1 1
[249,] 1 1 1 1 1 1 1
[250,] 1 1 1 1 1 1 1
[251,] 1 1 1 1 1 1 1
[252,] 1 1 1 1 1 1 1
[253,] 1 1 1 1 1 1 1
[254,] 1 1 1 1 1 1 1
[255,] 1 1 1 1 1 1 1
[256,] 1 1 1 1 1 1 1
[257,] 1 1 1 1 1 1 1
[258,] 1 1 1 1 1 1 1
[259,] 1 1 1 1 1 1 1
[260,] 1 1 1 1 1 1 1
[261,] 1 1 1 1 1 1 1
[262,] 1 1 1 1 1 1 1
[263,] 1 1 1 1 1 1 1
[264,] 1 1 1 1 1 1 1
[265,] 1 1 1 1 1 1 1
[266,] 1 1 1 1 1 1 1
[267,] 1 1 1 1 1 1 1
[268,] 1 1 1 1 1 1 1
[269,] 1 1 1 1 1 1 1
[270,] 1 1 1 1 1 1 1
[271,] 1 1 1 1 1 1 1
[272,] 1 1 1 1 1 1 1
[273,] 1 1 1 1 1 1 1
[274,] 1 1 1 1 1 1 1
[275,] 1 1 1 1 1 1 1
[276,] 1 1 1 1 1 1 1
[277,] 1 1 1 1 1 1 1
[278,] 1 1 1 1 1 1 1
[279,] 1 1 1 1 1 1 1
[280,] 1 1 1 1 1 1 1
[281,] 1 1 1 1 1 1 1
[282,] 1 1 1 1 1 1 1
[283,] 1 1 1 1 1 1 1
[284,] 1 1 1 1 1 1 1
[285,] 1 1 1 1 1 1 1
[286,] 1 1 1 1 1 1 1
[287,] 1 1 1 1 1 1 1
[288,] 1 1 1 1 1 1 1
[289,] 1 1 1 1 1 1 1
[290,] 1 1 1 1 1 1 1
[291,] 1 1 1 1 1 1 1
[292,] 1 1 1 1 1 1 1
[293,] 1 1 1 1 1 1 1
[294,] 1 1 1 1 1 1 1Some information that we now pass as initial value info (observations of alive) are actually known states, and could also be passed as data – in which case the initial values have to be 0.
Specify the parameters we wish to monitor.
parameters.to.save <- c("phi", "p")
parameters.to.save[1] "phi" "p" MCMC details.
n.iter <- 2500
n.burnin <- 1000
n.chains <- 2At last, let’s run nimble.
mcmc.phip <- nimbleMCMC(code = hmm.phip,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Examine the results.
library(MCMCvis)
MCMCsummary(mcmc.phip, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
p 0.89 0.03 0.83 0.90 0.94 1.01 262
phi 0.56 0.03 0.52 0.56 0.61 1.01 489Now a model with time-varying survival probabilities, and constant detection.
hmm.phitp <- nimbleCode({
for (t in 1:(T-1)){
phi[t] ~ dunif(0, 1) # prior survival
gamma[1,1,t] <- phi[t] # Pr(alive t -> alive t+1)
gamma[1,2,t] <- 1 - phi[t] # Pr(alive t -> dead t+1)
gamma[2,1,t] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,t] <- 1 # Pr(dead t -> dead t+1)
}
p ~ dunif(0, 1) # prior detection
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, j-1])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})The initial values.
initial.values <- function() list(phi = runif(my.constants$T-1,0,1),
p = runif(1,0,1),
z = zinits)Run nimble.
mcmc.phitp <- nimbleMCMC(code = hmm.phitp,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Display the results.
MCMCsummary(object = mcmc.phitp, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
p 0.89 0.03 0.82 0.89 0.94 1.02 255
phi[1] 0.62 0.11 0.41 0.62 0.81 1.01 567
phi[2] 0.46 0.06 0.34 0.46 0.59 1.00 652
phi[3] 0.49 0.06 0.37 0.49 0.60 1.00 561
phi[4] 0.62 0.06 0.51 0.62 0.73 1.01 572
phi[5] 0.61 0.05 0.50 0.61 0.71 1.01 546
phi[6] 0.59 0.06 0.48 0.58 0.70 1.02 428Now a model with time-varying detection and constant survival.
hmm.phipt <- nimbleCode({
phi ~ dunif(0, 1) # prior survival
gamma[1,1] <- phi # Pr(alive t -> alive t+1)
gamma[1,2] <- 1 - phi # Pr(alive t -> dead t+1)
gamma[2,1] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2] <- 1 # Pr(dead t -> dead t+1)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
for (t in 1:(T-1)){
p[t] ~ dunif(0, 1) # prior detection
omega[1,1,t] <- 1 - p[t] # Pr(alive t -> non-detected t)
omega[1,2,t] <- p[t] # Pr(alive t -> detected t)
omega[2,1,t] <- 1 # Pr(dead t -> non-detected t)
omega[2,2,t] <- 0 # Pr(dead t -> detected t)
}
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2])
y[i,j] ~ dcat(omega[z[i,j], 1:2, j-1])
}
}
})Initial values.
initial.values <- function() list(phi = runif(1,0,1),
p = runif(my.constants$T-1,0,1),
z = zinits)Run nimble.
mcmc.phipt <- nimbleMCMC(code = hmm.phipt,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Display the results.
MCMCsummary(object = mcmc.phipt, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
p[1] 0.75 0.11 0.50 0.76 0.93 1.03 540
p[2] 0.86 0.08 0.67 0.87 0.98 1.02 393
p[3] 0.84 0.07 0.70 0.85 0.97 1.01 400
p[4] 0.89 0.05 0.78 0.89 0.96 1.01 432
p[5] 0.91 0.04 0.82 0.92 0.98 1.04 381
p[6] 0.90 0.07 0.75 0.92 0.99 1.01 161
phi 0.56 0.03 0.51 0.56 0.61 1.00 522Eventually, the CJS model with both time-varying survival and recapture probabilities.
hmm.phitpt <- nimbleCode({
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
for (t in 1:(T-1)){ # loop over time
phi[t] ~ dunif(0, 1) # prior survival
gamma[1,1,t] <- phi[t] # Pr(alive t -> alive t+1)
gamma[1,2,t] <- 1 - phi[t] # Pr(alive t -> dead t+1)
gamma[2,1,t] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,t] <- 1 # Pr(dead t -> dead t+1)
p[t] ~ dunif(0, 1) # prior detection
omega[1,1,t] <- 1 - p[t] # Pr(alive t -> non-detected t)
omega[1,2,t] <- p[t] # Pr(alive t -> detected t)
omega[2,1,t] <- 1 # Pr(dead t -> non-detected t)
omega[2,2,t] <- 0 # Pr(dead t -> detected t)
}
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, j-1])
y[i,j] ~ dcat(omega[z[i,j], 1:2, j-1])
}
}
})Initial values.
initial.values <- function() list(phi = runif(my.constants$T-1,0,1),
p = runif(my.constants$T-1,0,1),
z = zinits)Run nimble.
mcmc.phitpt <- nimbleMCMC(code = hmm.phitpt,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Display the numerical summaries. Note the small effective sample size for the last survival and recapture probabilities.
MCMCsummary(object = mcmc.phitpt, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
p[1] 0.67 0.14 0.39 0.68 0.91 1.01 280
p[2] 0.87 0.08 0.67 0.88 0.98 1.01 269
p[3] 0.88 0.06 0.73 0.90 0.98 1.02 294
p[4] 0.88 0.06 0.75 0.88 0.96 1.01 305
p[5] 0.91 0.05 0.80 0.92 0.98 1.02 234
p[6] 0.77 0.14 0.51 0.77 0.98 2.13 40
phi[1] 0.72 0.13 0.46 0.72 0.99 1.01 246
phi[2] 0.45 0.07 0.31 0.44 0.60 1.01 457
phi[3] 0.48 0.06 0.37 0.48 0.60 1.00 470
phi[4] 0.63 0.06 0.51 0.63 0.75 1.00 492
phi[5] 0.60 0.06 0.48 0.59 0.71 1.03 457
phi[6] 0.69 0.14 0.47 0.67 0.98 2.10 36Caterpillar plot of the estimates.
MCMCplot(object = mcmc.phitpt)
Let’s focus for a minute on the last survival probability. See how mixing is bad and the overlap with the prior is big. This parameter is redundant, and it can be shown that only the product of \(\phi_6\) and \(p_7\) can be estimated.
priors <- runif(3000, 0, 1)
MCMCtrace(object = mcmc.phitpt,
ISB = FALSE,
exact = TRUE,
params = c("phi[6]"),
pdf = FALSE,
priors = priors)
Model selection with WAIC
We re-run the four models above, but now we make sure we monitor the \(z\) and set the WAIC argument to TRUE in nimbleMCMC().
my.constants <- list(N = nrow(y), T = ncol(y), first = first)
my.data <- list(y = y + 1)
initial.values <- function() list(phi = runif(1,0,1),
p = runif(1,0,1),
z = zinits)
parameters.to.save <- c("phi", "p", "z")
mcmc.phip <- nimbleMCMC(code = hmm.phip,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
WAIC = TRUE)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
Monitored nodes are valid for WAIC.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|initial.values <- function() list(phi = runif(my.constants$T-1,0,1),
p = runif(1,0,1),
z = zinits)
mcmc.phitp <- nimbleMCMC(code = hmm.phitp,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
WAIC = TRUE)defining model...
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
Monitored nodes are valid for WAIC.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|initial.values <- function() list(phi = runif(1,0,1),
p = runif(my.constants$T-1,0,1),
z = zinits)
mcmc.phipt <- nimbleMCMC(code = hmm.phipt,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
WAIC = TRUE)defining model...
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
Monitored nodes are valid for WAIC.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|initial.values <- function() list(phi = runif(my.constants$T-1,0,1),
p = runif(my.constants$T-1,0,1),
z = zinits)
mcmc.phitpt <- nimbleMCMC(code = hmm.phitpt,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
WAIC = TRUE)defining model...
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
Monitored nodes are valid for WAIC.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Now we report model ranking.
data.frame(model = c("(phi,p)",
"(phit,p)",
"(phi,pt)",
"(phit,pt)"),
WAIC = c(mcmc.phip$WAIC,
mcmc.phitp$WAIC,
mcmc.phipt$WAIC,
mcmc.phitpt$WAIC)) model WAIC
1 (phi,p) 267.4950
2 (phit,p) 274.3351
3 (phi,pt) 269.4981
4 (phit,pt) 311.2689Add a temporal covariate
Now we’d like to add a temporal covariate to try and explain annual variation in survival. We pick water flow in river. We specify the relationship on the logit scale.
hmm.phiflowp <- nimbleCode({
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
for (t in 1:(T-1)){
logit(phi[t]) <- beta[1] + beta[2] * flow[t]
gamma[1,1,t] <- phi[t] # Pr(alive t -> alive t+1)
gamma[1,2,t] <- 1 - phi[t] # Pr(alive t -> dead t+1)
gamma[2,1,t] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,t] <- 1 # Pr(dead t -> dead t+1)
}
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
beta[1] ~ dnorm(0, 1.5) # prior intercept
beta[2] ~ dnorm(0, 1.5) # prior slope
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, j-1])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})We only take the values we need, and standardize the covariate.
# water flow in L/s
water_flow <- c(443, 1114, 529, 434, 627, 466, 730) # 1981, 1982, ..., 1987
water_flow <- water_flow[-7]
water_flow_st <- (water_flow - mean(water_flow))/sd(water_flow)Constants in a list.
my.constants <- list(N = nrow(y),
T = ncol(y),
first = first,
flow = water_flow_st)
my.constants$N
[1] 294
$T
[1] 7
$first
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3
[75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[112] 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[149] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[186] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6
[223] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7
[260] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
$flow
[1] -0.61028761 1.96250602 -0.28054059 -0.64479602 0.09521765 -0.52209945Initial values.
initial.values <- function() list(beta = rnorm(2,0,1),
p = runif(1,0,1),
z = zinits)Parameters to be monitored.
parameters.to.save <- c("beta", "p", "phi")
parameters.to.save[1] "beta" "p" "phi" Run nimble.
mcmc.phiflowp <- nimbleMCMC(code = hmm.phiflowp,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Caterpillar plot of the regression parameters. The posterior distribution of the slope is centered on negative values, suggesting the as water flow increases, survival decreases.
MCMCplot(object = mcmc.phiflowp, params = "beta", ISB = TRUE)
Caterpillar plot of the survival estimates. Survival between 1982 and 1983 seems to have been affected highly by a huge water flow compared to the other years.
MCMCplot(object = mcmc.phiflowp, params = "phi", ISB = TRUE)
Yearly random effects
We may wish to allow for extra variation in the survival vs. water flow relationship. To do so, we consider a yearly random effect. The prior on the standard deviation of the random effect is uniform between 0 and 10.
hmm.phiflowREpt <- nimbleCode({
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
for (t in 1:(T-1)){
logit(phi[t]) <- beta[1] + beta[2] * flow[t] + eps[t] # eps is random effect
eps[t] ~ dnorm(0, sd = sdeps)
gamma[1,1,t] <- phi[t] # Pr(alive t -> alive t+1)
gamma[1,2,t] <- 1 - phi[t] # Pr(alive t -> dead t+1)
gamma[2,1,t] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,t] <- 1 # Pr(dead t -> dead t+1)
p[t] ~ dunif(0, 1) # prior detection
omega[1,1,t] <- 1 - p[t] # Pr(alive t -> non-detected t)
omega[1,2,t] <- p[t] # Pr(alive t -> detected t)
omega[2,1,t] <- 1 # Pr(dead t -> non-detected t)
omega[2,2,t] <- 0 # Pr(dead t -> detected t)
}
beta[1] ~ dnorm(0, 1.5) # prior intercept
beta[2] ~ dnorm(0, 1.5) # prior slope
sdeps ~ dunif(0,10)
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, j-1])
y[i,j] ~ dcat(omega[z[i,j], 1:2, j-1])
}
}
})Initial values.
initial.values <- function() list(beta = rnorm(2,0,1),
p = runif(my.constants$T-1,0,1),
sdeps = runif(1,0,3),
z = zinits)Parameters to be monitored.
parameters.to.save <- c("beta", "p", "phi", "sdeps")MCMC details. Note that we’ve increased the number of iterations and the length of the burn-in period.
n.iter <- 10000
n.burnin <- 5000
n.chains <- 2Run nimble.
mcmc.phiflowREpt <- nimbleMCMC(code = hmm.phiflowREpt,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
checking model calculations...NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Display outputs. Seems that the water flow effect is not so important anymore.
MCMCsummary(object = mcmc.phiflowREpt, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.31 0.22 -0.12 0.30 0.79 1.16 210
beta[2] -0.25 0.23 -0.77 -0.24 0.19 1.04 314
p[1] 0.70 0.13 0.44 0.71 0.92 1.00 1438
p[2] 0.88 0.08 0.68 0.89 0.98 1.05 982
p[3] 0.87 0.07 0.71 0.88 0.97 1.01 1078
p[4] 0.88 0.05 0.76 0.88 0.96 1.00 1224
p[5] 0.91 0.05 0.79 0.91 0.98 1.00 1012
p[6] 0.84 0.10 0.62 0.85 0.99 1.04 223
phi[1] 0.64 0.09 0.50 0.63 0.84 1.08 353
phi[2] 0.45 0.07 0.32 0.45 0.59 1.02 1546
phi[3] 0.53 0.06 0.41 0.54 0.63 1.16 465
phi[4] 0.62 0.05 0.52 0.62 0.72 1.02 1315
phi[5] 0.59 0.05 0.50 0.58 0.69 1.06 1122
phi[6] 0.62 0.08 0.50 0.61 0.81 1.05 237
sdeps 0.37 0.37 0.02 0.28 1.21 1.33 89Trace plots for the standard deviation of the random effect.
MCMCtrace(object = mcmc.phiflowREpt, params = "sdeps", pdf = FALSE)
Add a discrete individual covariate (aka group)
OK now we’re gonna illustrate how to have a discrete individual covariate, which we often call group. Here we will consider the sex of individuals. There are two methods to include sex as a covariate.
First, let us define the covariate sex that takes value 0 if the individual is a male, and 1 if it is a female. We put these values in a vector sex.
sex <- if_else(dipper$sex == "M", 0, 1)
sex [1] 0 1 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 0 1 1 1 1 1 1 1 0 1 1 0 1 1 0 0 0 0 0 0
[38] 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
[75] 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1
[149] 1 1 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0
[186] 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1
[223] 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0
[260] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Now we write the model, and write the survival as a function of the covariate sex on the logit scale \(\text{logit}(\phi_i) = \beta_1 + \beta_2 * \text{sex}_i\). We need to make the matrix \(\Gamma\) individual specific. For convenience we also define \(\phi_{\text{male}} = \beta_1\) with \(\text{sex}_i = 0\) and \(\phi_{\text{female}} = \beta_1 + \beta_2\) with \(\text{sex}_i = 1\).
hmm.phisexp <- nimbleCode({
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){ # loop over individuals
logit(phi[i]) <- beta[1] + beta[2] * sex[i]
gamma[1,1,i] <- phi[i] # Pr(alive t -> alive t+1)
gamma[1,2,i] <- 1 - phi[i] # Pr(alive t -> dead t+1)
gamma[2,1,i] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,i] <- 1 # Pr(dead t -> dead t+1)
}
beta[1] ~ dnorm(mean = 0, sd = 1.5)
beta[2] ~ dnorm(mean = 0, sd = 1.5)
phi_male <- 1/(1+exp(-beta[1]))
phi_female <- 1/(1+exp(-(beta[1]+beta[2])))
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, i])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})Constants in a list.
my.constants <- list(N = nrow(y),
T = ncol(y),
first = first,
sex = sex)Data in a list.
my.data <- list(y = y + 1)Initial values.
initial.values <- function() list(beta = rnorm(2,0,1),
p = runif(1,0,1),
z = zinits)Parameters to monitor.
parameters.to.save <- c("beta", "p", "phi_male", "phi_female")MCMC details.
n.iter <- 2500
n.burnin <- 1000
n.chains <- 2Run nimble.
mcmc.phisexp <- nimbleMCMC(code = hmm.phisexp,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Display results.
MCMCsummary(object = mcmc.phisexp, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.29 0.15 0.00 0.28 0.58 1.07 190
beta[2] -0.07 0.20 -0.47 -0.07 0.32 1.05 227
p 0.89 0.03 0.83 0.90 0.94 1.00 230
phi_female 0.55 0.03 0.48 0.55 0.62 1.00 678
phi_male 0.57 0.04 0.50 0.57 0.64 1.07 191Another method to include a group effect is to use nested indexing. Let’s use a covariate \(\text{sex}\) that contains 1s and 2s, indicating the sex of each individual: 1 if male, and 2 if female. E.g. for individual \(i = 2\), beta[sex[i]] gives beta[sex[2]] which will be beta[1] or beta[2] depending on whether sex[2] is 1 or 2.
sex <- if_else(dipper$sex == "M", 1, 2)
sex [1] 1 2 1 2 2 1 1 1 1 2 2 1 2 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 1 2 2 1 1 1 1 1 1
[38] 1 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1
[75] 2 1 1 1 1 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2
[112] 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 2 2 1 1 1 2 2 2 2 1 1 1 1 1 1 2 2 2 2
[149] 2 2 1 2 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1
[186] 1 1 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2
[223] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1
[260] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2Write the model.
hmm.phisexpt.ni <- nimbleCode({
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){
phi[i] <- beta[sex[i]]
gamma[1,1,i] <- phi[i] # Pr(alive t -> alive t+1)
gamma[1,2,i] <- 1 - phi[i] # Pr(alive t -> dead t+1)
gamma[2,1,i] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,i] <- 1 # Pr(dead t -> dead t+1)
}
beta[1] ~ dunif(0,1)
beta[2] ~ dunif(0,1)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, i])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})My constants.
my.constants <- list(N = nrow(y),
T = ncol(y),
first = first,
sex = sex) # beta[1] male survival
# beta[2] female survivalParameters to monitor.
parameters.to.save <- c("beta", "p")Initial values.
initial.values <- function() list(beta = runif(2,0,1),
p = runif(1,0,1),
z = zinits)Run nimble.
mcmc.phisexp.ni <- nimbleMCMC(code = hmm.phisexpt.ni,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Dislpay results. Compare with the other method above, the estimates are very similar.
MCMCsummary(object = mcmc.phisexp.ni, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.57 0.04 0.50 0.57 0.64 1.01 615
beta[2] 0.55 0.03 0.48 0.55 0.62 1.01 575
p 0.90 0.03 0.83 0.90 0.95 1.01 258Add a continuous individual covariate
Besides discrete individual covariates, you might want to have continuous individual covariates, e.g. wing length in the dipper case study. Note that we’re considering an individual trait that takes the same value whatever the occasion. If we were to have time-varying individual covariate in the model, we would have to do something about missing values of the covariate when an individual is not recaptured. The easiest way to cope with time-varying individual covariate is to discretize and treat levels of the covariates as states. More in the next live demo. Back to wing length. We first standardize the covariate.
wing.length.st <- as.vector(scale(dipper$wing_length))
head(wing.length.st)[1] 0.7581335 -0.8671152 0.5259551 -1.5636504 -1.3314720 1.2224903Now we write the model. Basically we replace sex by wing length in the first method we used in the previous section. Easy.
hmm.phiwlp <- nimbleCode({
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){
logit(phi[i]) <- beta[1] + beta[2] * winglength[i]
gamma[1,1,i] <- phi[i] # Pr(alive t -> alive t+1)
gamma[1,2,i] <- 1 - phi[i] # Pr(alive t -> dead t+1)
gamma[2,1,i] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,i] <- 1 # Pr(dead t -> dead t+1)
}
beta[1] ~ dnorm(mean = 0, sd = 1.5)
beta[2] ~ dnorm(mean = 0, sd = 1.5)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, i])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})Constants in a list.
my.constants <- list(N = nrow(y),
T = ncol(y),
first = first,
winglength = wing.length.st)Initial values.
initial.values <- function() list(beta = rnorm(2,0,1),
p = runif(1,0,1),
z = zinits)Run nimble.
mcmc.phiwlp <- nimbleMCMC(code = hmm.phiwlp,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions...
checking model calculations...
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Numerical summaries. Wing length does not seem to explain much individual-to-individual variation in survival.
MCMCsummary(mcmc.phiwlp, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.24 0.10 0.02 0.25 0.45 1.01 528
beta[2] -0.02 0.09 -0.21 -0.02 0.15 1.00 674
p 0.90 0.03 0.84 0.90 0.95 1.00 259Let’s plot the relationship. First, we gather the values generated from the posterior distribution of the regression parameters in the two chains.
beta1 <- c(mcmc.phiwlp$chain1[,'beta[1]'], mcmc.phiwlp$chain2[,'beta[1]'])
beta2 <- c(mcmc.phiwlp$chain1[,'beta[2]'], mcmc.phiwlp$chain2[,'beta[2]'])Then we define a grid of values for wing length, and predict survival for each MCMC iteration.
predicted_survival <- matrix(NA, nrow = length(beta1), ncol = length(my.constants$winglength))
for (i in 1:length(beta1)){
for (j in 1:length(my.constants$winglength)){
predicted_survival[i,j] <- plogis(beta1[i] + beta2[i] * my.constants$winglength[j])
}
}Now we calculate posterior mean and the credible interval. Note the ordering.
mean_survival <- apply(predicted_survival, 2, mean)
lci <- apply(predicted_survival, 2, quantile, prob = 2.5/100)
uci <- apply(predicted_survival, 2, quantile, prob = 97.5/100)
ord <- order(my.constants$winglength)
df <- data.frame(wing_length = my.constants$winglength[ord],
survival = mean_survival[ord],
lci = lci[ord],
uci = uci[ord])Now time to visualize.
df %>%
ggplot() +
aes(x = wing_length, y = survival) +
geom_line() +
geom_ribbon(aes(ymin = lci, ymax = uci), fill = "grey70", alpha = 0.5) +
ylim(0,1) +
labs(x = "wing length", y = "estimated survival")
Add an individual random effect on top of the individual covariate
We add an individual random effect.
hmm.phiwlrep <- nimbleCode({
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){
logit(phi[i]) <- beta[1] + beta[2] * winglength[i] + eps[i]
eps[i] ~ dnorm(mean = 0, sd = sdeps)
gamma[1,1,i] <- phi[i] # Pr(alive t -> alive t+1)
gamma[1,2,i] <- 1 - phi[i] # Pr(alive t -> dead t+1)
gamma[2,1,i] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,i] <- 1 # Pr(dead t -> dead t+1)
}
beta[1] ~ dnorm(mean = 0, sd = 1.5)
beta[2] ~ dnorm(mean = 0, sd = 1.5)
sdeps ~ dunif(0, 10)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, i])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})Initial values.
initial.values <- function() list(beta = rnorm(2,0,1.5),
sdeps = runif(1,0,3),
p = runif(1,0,1),
z = zinits)Parameters to be monitored.
parameters.to.save <- c("beta", "sdeps", "p")MCMC details. Note that we increase the number of iterations and the length of the burn-in period.
n.iter <- 10000
n.burnin <- 5000
n.chains <- 2Run nimble.
mcmc.phiwlrep <- nimbleMCMC(code = hmm.phiwlrep,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
checking model calculations...NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Numerical summaries.
MCMCsummary(mcmc.phiwlrep, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.18 0.12 -0.07 0.19 0.41 1.07 700
beta[2] 0.00 0.11 -0.21 0.00 0.21 1.02 1626
p 0.90 0.03 0.83 0.90 0.95 1.00 863
sdeps 0.55 0.22 0.12 0.53 1.03 2.07 16Let’s try something else. We reparameterize by non-centering.
hmm.phiwlrep <- nimbleCode({
p ~ dunif(0, 1) # prior detection
omega[1,1] <- 1 - p # Pr(alive t -> non-detected t)
omega[1,2] <- p # Pr(alive t -> detected t)
omega[2,1] <- 1 # Pr(dead t -> non-detected t)
omega[2,2] <- 0 # Pr(dead t -> detected t)
for (i in 1:N){
logit(phi[i]) <- beta[1] + beta[2] * winglength[i] + sdeps * eps[i]
eps[i] ~ dnorm(mean = 0, sd = 1)
gamma[1,1,i] <- phi[i] # Pr(alive t -> alive t+1)
gamma[1,2,i] <- 1 - phi[i] # Pr(alive t -> dead t+1)
gamma[2,1,i] <- 0 # Pr(dead t -> alive t+1)
gamma[2,2,i] <- 1 # Pr(dead t -> dead t+1)
}
beta[1] ~ dnorm(mean = 0, sd = 1.5)
beta[2] ~ dnorm(mean = 0, sd = 1.5)
sdeps ~ dunif(0, 10)
delta[1] <- 1 # Pr(alive t = 1) = 1
delta[2] <- 0 # Pr(dead t = 1) = 0
# likelihood
for (i in 1:N){
z[i,first[i]] ~ dcat(delta[1:2])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, i])
y[i,j] ~ dcat(omega[z[i,j], 1:2])
}
}
})Initial values.
initial.values <- function() list(beta = rnorm(2,0,1.5),
sdeps = runif(1,0,3),
p = runif(1,0,1),
z = zinits)Parameters to be monitored.
parameters.to.save <- c("beta", "sdeps", "p")MCMC details. Note that we increase the number of iterations and the length of the burn-in period.
n.iter <- 10000
n.burnin <- 5000
n.chains <- 2Run nimble.
mcmc.phiwlrep <- nimbleMCMC(code = hmm.phiwlrep,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)defining model...building model...setting data and initial values...running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
checking model calculations...NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|Numerical summaries. Much better.
MCMCsummary(mcmc.phiwlrep, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
beta[1] 0.21 0.12 -0.04 0.21 0.44 1.00 575
beta[2] -0.01 0.10 -0.21 -0.02 0.19 1.01 1729
p 0.90 0.03 0.83 0.90 0.95 1.00 811
sdeps 0.40 0.27 0.01 0.36 1.01 1.01 165Let’s plot the posterior distribution of the standard deviation of the individual random effect.
sdeps <- c(mcmc.phiwlrep$chain1[,"sdeps"], mcmc.phiwlrep$chain2[,"sdeps"])
sdeps %>%
as_tibble() %>%
ggplot() +
aes(x = value) +
geom_histogram(color = "white", binwidth = .03, fill = "gray70") +
geom_density(aes(y = .03 * ..count..))