Class 5 live demo: Survival estimation
last updated: 2023-06-25
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.5275765
$p
[1] 0.554956
$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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.02 247
phi 0.56 0.03 0.51 0.56 0.61 1.00 524Now 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.83 0.90 0.95 1.04 240
phi[1] 0.62 0.10 0.43 0.63 0.81 1.01 612
phi[2] 0.45 0.06 0.33 0.45 0.58 1.00 588
phi[3] 0.48 0.05 0.38 0.48 0.59 1.00 626
phi[4] 0.63 0.05 0.52 0.63 0.74 1.00 646
phi[5] 0.61 0.05 0.51 0.61 0.72 1.00 489
phi[6] 0.58 0.06 0.48 0.58 0.70 1.01 453Now 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.49 0.76 0.92 1.00 508
p[2] 0.85 0.08 0.66 0.86 0.98 1.00 292
p[3] 0.84 0.07 0.69 0.84 0.96 1.01 346
p[4] 0.88 0.05 0.77 0.89 0.96 1.00 386
p[5] 0.91 0.05 0.80 0.92 0.98 1.02 288
p[6] 0.90 0.06 0.75 0.91 0.99 1.06 184
phi 0.56 0.03 0.51 0.56 0.61 1.02 449Eventually, 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.66 0.14 0.39 0.67 0.90 1.04 299
p[2] 0.87 0.08 0.68 0.88 0.98 1.04 296
p[3] 0.87 0.07 0.70 0.88 0.97 1.01 239
p[4] 0.87 0.06 0.75 0.88 0.97 1.00 343
p[5] 0.89 0.06 0.77 0.90 0.98 1.01 280
p[6] 0.80 0.13 0.55 0.81 1.00 1.03 32
phi[1] 0.73 0.13 0.47 0.73 0.97 1.05 170
phi[2] 0.45 0.07 0.32 0.44 0.60 1.01 390
phi[3] 0.49 0.06 0.38 0.48 0.62 1.01 429
phi[4] 0.62 0.06 0.50 0.62 0.74 1.03 382
phi[5] 0.61 0.06 0.50 0.60 0.72 1.03 451
phi[6] 0.66 0.12 0.46 0.64 0.92 1.02 43Caterpillar 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
[Warning] There are individual pWAIC values that are greater than 0.4. This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017), at least in cases without grouping of data nodes or multivariate data nodes.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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
[Warning] There are individual pWAIC values that are greater than 0.4. This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017), at least in cases without grouping of data nodes or multivariate data nodes.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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
[Warning] There are individual pWAIC values that are greater than 0.4. This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017), at least in cases without grouping of data nodes or multivariate data nodes.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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.running chain 1...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|running chain 2...|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
[Warning] There are individual pWAIC values that are greater than 0.4. This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017), at least in cases without grouping of data nodes or multivariate data nodes.Now we report model ranking.
data.frame(model = c("(phi,p)",
"(phit,p)",
"(phi,pt)",
"(phit,pt)"),
WAIC = c(mcmc.phip$WAIC$WAIC,
mcmc.phitp$WAIC$WAIC,
mcmc.phipt$WAIC$WAIC,
mcmc.phitpt$WAIC$WAIC)) model WAIC
1 (phi,p) 269.7288
2 (phit,p) 275.2350
3 (phi,pt) 271.6999
4 (phit,pt) 316.4553Add 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensions [Note] 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[Note] NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.Compiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.29 0.23 -0.12 0.27 0.86 1.00 170
beta[2] -0.24 0.24 -0.70 -0.24 0.26 1.02 255
p[1] 0.71 0.13 0.45 0.72 0.92 1.00 1401
p[2] 0.88 0.08 0.69 0.89 0.98 1.01 968
p[3] 0.87 0.07 0.71 0.88 0.97 1.01 889
p[4] 0.88 0.05 0.75 0.89 0.96 1.01 1248
p[5] 0.91 0.05 0.80 0.92 0.98 1.00 975
p[6] 0.84 0.12 0.56 0.85 1.00 1.00 172
phi[1] 0.64 0.08 0.49 0.63 0.84 1.00 528
phi[2] 0.45 0.07 0.33 0.45 0.59 1.00 1986
phi[3] 0.53 0.06 0.41 0.54 0.63 1.01 394
phi[4] 0.62 0.05 0.52 0.62 0.72 1.00 1342
phi[5] 0.59 0.05 0.50 0.58 0.69 1.01 1008
phi[6] 0.63 0.09 0.50 0.61 0.90 1.00 140
sdeps 0.37 0.30 0.03 0.28 1.20 1.01 105Trace 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.01 0.29 0.60 1.01 224
beta[2] -0.08 0.21 -0.54 -0.08 0.31 1.02 189
p 0.90 0.03 0.84 0.90 0.95 1.01 268
phi_female 0.55 0.04 0.48 0.55 0.62 1.00 444
phi_male 0.57 0.04 0.50 0.57 0.65 1.01 226Another 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.03 0.50 0.57 0.63 1 770
beta[2] 0.55 0.03 0.49 0.55 0.61 1 599
p 0.90 0.03 0.84 0.90 0.95 1 261Add 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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensionsChecking model calculationsCompiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.25 0.10 0.04 0.25 0.45 1.01 533
beta[2] -0.02 0.10 -0.21 -0.02 0.18 1.02 683
p 0.89 0.03 0.83 0.90 0.94 1.01 267Let’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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensions [Note] 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[Note] NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.Compiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.06 0.17 0.42 1.00 1140
beta[2] 0.00 0.11 -0.20 0.00 0.21 1.00 1584
p 0.90 0.03 0.84 0.90 0.95 1.00 916
sdeps 0.60 0.16 0.35 0.58 0.93 1.03 41Let’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 modelBuilding modelSetting data and initial valuesRunning calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.Checking model sizes and dimensions [Note] 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[Note] NAs were detected in model variables: eps, logProb_eps, phi, gamma, logProb_z.Compiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.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.20 0.12 -0.05 0.20 0.43 1.00 867
beta[2] 0.00 0.10 -0.20 0.00 0.19 1.00 1694
p 0.90 0.03 0.83 0.90 0.95 1.02 770
sdeps 0.43 0.27 0.02 0.41 0.97 1.00 180Let’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..))Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.