Introduction

In this demo, we introduce Markov models and hidden Markov models following up on the survival example. Now we’ll be using longitudinal data, that is data collected over time on the same individuals. We will also briefly see how to simulate data. Simulating data often proves useful to better understand how your model works, check that you get the right answer by comparing the parameters you used to simulate the data and the estimates you get, and to communicate with others on the model hypotheses and limitations.

Simple survival example

First, we load nimble:

library(nimble)

We sart with a relatively simple model. We write down the code of our model in BUGS syntax:

naive.survival.model <- nimbleCode({
  # prior
  phi ~ dunif(0, 1)
  # likelihood
  y ~ dbinom(phi, n)
})

We read in the data and constants:

my.constants <- list(n = 57)
my.data <- list(y = 19)

We pecify initial values:

initial.values <- function() list(phi = runif(1,0,1))
initial.values()
$phi
[1] 0.808186

Which parameters to save?

parameters.to.save <- c("phi")

We specifiy MCMC details:

n.iter <- 5000
n.burnin <- 1000
n.chains <- 2
n.thin <- 1

And, we run our model, tadaa!

mcmc.output <- nimbleMCMC(code = naive.survival.model,     
                          data = my.data,  
                          constants = my.constants,
                          inits = initial.values,
                          monitors = parameters.to.save,
                          thin = n.thin,
                          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...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Let’s explore MCMC outputs:

str(mcmc.output)
List of 2
 $ chain1: num [1:4000, 1] 0.318 0.318 0.318 0.383 0.383 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr "phi"
 $ chain2: num [1:4000, 1] 0.414 0.262 0.268 0.268 0.268 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr "phi"
head(mcmc.output$chain1)
           phi
[1,] 0.3182470
[2,] 0.3182470
[3,] 0.3182470
[4,] 0.3826353
[5,] 0.3826353
[6,] 0.3826353

Let’s calculate some numerical summaries:

library(MCMCvis)
MCMCsummary(mcmc.output, round = 2)
    mean   sd 2.5%  50% 97.5% Rhat n.eff
phi 0.34 0.06 0.23 0.34  0.47    1  1954

Let’s get the trace and posterior density of survival:

MCMCtrace(mcmc.output,
          pdf = FALSE)

MCMCtrace(mcmc.output,
          pdf = FALSE,
          ind = TRUE,
          Rhat = TRUE,
          n.eff = TRUE)

Markov chain

Load some useful packages.

library(tidyverse)

We start with 57 individuals, and we monitor them over 5 occasions.

nind <- 57
nocc <- 5

For simplicity, we assume that we have a single cohort, that is all individuals enter the study at the same time, here on the first occasion.

first <- rep(1, nind) # single cohort

We set survival to 0.8, and define \(z\) as a matrix with dimensions the number of individuals (rows) and the number of occasions (columns).

phi <- 0.8
z <- matrix(NA, nrow = nind, ncol = nocc)

Now we simulate the fate of all individuals over time.

for (i in 1:nind){ # loop over individuals
  z[i,first[i]] <- 1 # all individuals are alive at first occasion
  for (t in (first[i]+1):nocc){ # loop over time
    z[i,t] <- rbinom(1, 1, phi * z[i,t-1]) # if z[i,t-1] = 1, then z[i,t] ~ dbern(phi)
                                           # if z[i,t-1] = 0, then z[i,t] ~ dbern(0) = 0 (once you're dead, you remain dead)
  }
}

The zeros are replaced by twos to match the coding proposed in the lecture. One is for alive, two is for dead.

z[z==0] <- 2 # 2 = dead, 1 = alive

Name the columns.

colnames(z) <- paste0("winter ", 1:nocc)

Display the matrix \(z\).

z %>% 
  as_tibble() %>% 
  add_column(id = 1:nind, .before = "winter 1") %>%
  kableExtra::kable() %>%
  kableExtra::scroll_box(width = "100%", height = "400px")
id winter 1 winter 2 winter 3 winter 4 winter 5
1 1 2 2 2 2
2 1 1 1 1 1
3 1 1 1 1 2
4 1 1 1 1 1
5 1 1 1 1 1
6 1 2 2 2 2
7 1 2 2 2 2
8 1 1 2 2 2
9 1 2 2 2 2
10 1 1 1 1 1
11 1 1 1 2 2
12 1 1 1 1 2
13 1 1 1 1 1
14 1 2 2 2 2
15 1 2 2 2 2
16 1 1 1 1 1
17 1 2 2 2 2
18 1 1 1 1 2
19 1 1 2 2 2
20 1 2 2 2 2
21 1 1 2 2 2
22 1 2 2 2 2
23 1 1 2 2 2
24 1 1 2 2 2
25 1 1 1 1 1
26 1 1 1 1 2
27 1 1 1 1 1
28 1 1 1 1 2
29 1 1 2 2 2
30 1 1 1 1 1
31 1 1 1 1 1
32 1 1 1 2 2
33 1 1 1 1 1
34 1 2 2 2 2
35 1 1 1 1 1
36 1 1 1 1 1
37 1 1 1 1 1
38 1 1 1 1 2
39 1 1 1 1 1
40 1 1 1 2 2
41 1 1 1 1 2
42 1 1 1 1 2
43 1 1 1 1 1
44 1 1 1 1 1
45 1 1 1 1 1
46 1 1 1 1 1
47 1 2 2 2 2
48 1 1 1 1 2
49 1 1 1 1 1
50 1 1 2 2 2
51 1 1 1 1 1
52 1 2 2 2 2
53 1 1 2 2 2
54 1 1 1 1 1
55 1 1 1 1 1
56 1 1 1 1 2
57 1 2 2 2 2

Now we’re going to fit a Markov model to the data we’ve just simulated. We load nimble.

library(nimble)

Then we build the model.

markov.survival <- nimbleCode({
  phi ~ dunif(0, 1) # prior
  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
  # likelihood
  for (i in 1:N){
    z[i,1] ~ dcat(delta[1:2])
    for (j in 2:T){
      z[i,j] ~ dcat(gamma[z[i,j-1], 1:2])
    }
  }})

]

We put the constants in a list.

my.constants <- list(N = 57, T = 5)
my.constants
$N
[1] 57

$T
[1] 5

We put the data in a list.

my.data <- list(z = z)
my.data
$z
      winter 1 winter 2 winter 3 winter 4 winter 5
 [1,]        1        2        2        2        2
 [2,]        1        1        1        1        1
 [3,]        1        1        1        1        2
 [4,]        1        1        1        1        1
 [5,]        1        1        1        1        1
 [6,]        1        2        2        2        2
 [7,]        1        2        2        2        2
 [8,]        1        1        2        2        2
 [9,]        1        2        2        2        2
[10,]        1        1        1        1        1
[11,]        1        1        1        2        2
[12,]        1        1        1        1        2
[13,]        1        1        1        1        1
[14,]        1        2        2        2        2
[15,]        1        2        2        2        2
[16,]        1        1        1        1        1
[17,]        1        2        2        2        2
[18,]        1        1        1        1        2
[19,]        1        1        2        2        2
[20,]        1        2        2        2        2
[21,]        1        1        2        2        2
[22,]        1        2        2        2        2
[23,]        1        1        2        2        2
[24,]        1        1        2        2        2
[25,]        1        1        1        1        1
[26,]        1        1        1        1        2
[27,]        1        1        1        1        1
[28,]        1        1        1        1        2
[29,]        1        1        2        2        2
[30,]        1        1        1        1        1
[31,]        1        1        1        1        1
[32,]        1        1        1        2        2
[33,]        1        1        1        1        1
[34,]        1        2        2        2        2
[35,]        1        1        1        1        1
[36,]        1        1        1        1        1
[37,]        1        1        1        1        1
[38,]        1        1        1        1        2
[39,]        1        1        1        1        1
[40,]        1        1        1        2        2
[41,]        1        1        1        1        2
[42,]        1        1        1        1        2
[43,]        1        1        1        1        1
[44,]        1        1        1        1        1
[45,]        1        1        1        1        1
[46,]        1        1        1        1        1
[47,]        1        2        2        2        2
[48,]        1        1        1        1        2
[49,]        1        1        1        1        1
[50,]        1        1        2        2        2
[51,]        1        1        1        1        1
[52,]        1        2        2        2        2
[53,]        1        1        2        2        2
[54,]        1        1        1        1        1
[55,]        1        1        1        1        1
[56,]        1        1        1        1        2
[57,]        1        2        2        2        2

We write a function that generates initial values for survival.

initial.values <- function() list(phi = runif(1,0,1))
initial.values()
$phi
[1] 0.4641216

We specify that we’d like to monitor survival.

parameters.to.save <- c("phi")
parameters.to.save
[1] "phi"

Let’s specify a few details about the chains.

n.iter <- 5000
n.burnin <- 1000
n.chains <- 2

And now, run nimble.

mcmc.output <- nimbleMCMC(code = markov.survival, 
                          constants = my.constants,
                          data = my.data,              
                          inits = initial.values,
                          monitors = parameters.to.save,
                          niter = n.iter, 
                          nburnin = n.burnin, 
                          nchains = n.chains)
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
library(MCMCvis)
MCMCsummary(mcmc.output, round = 2)
    mean   sd 2.5% 50% 97.5% Rhat n.eff
phi  0.8 0.03 0.73 0.8  0.85    1  1700

The posterior mean is close to the value of survival we used to simulate the data \(\phi = 0.8\). Great!

Note that you should be able to write the model in a more efficient way with matrices and vectors.

markov.survival <- nimbleCode({
  phi ~ dunif(0, 1) # prior
  gamma[1:2,1:2] <- matrix(c(phi, 0, 1 - phi, 1), nrow = 2)
  delta[1:2] <- c(1, 0)
  # likelihood
  for (i in 1:N){
    z[i,1] ~ dcat(delta[1:2])
    for (j in 2:T){
      z[i,j] ~ dcat(gamma[z[i,j-1], 1:2])
    }
  }})

]

Hiden Markov chain

On top of the matrix of alive/dead states, we add the detection process. First we set the detection probability to \(p = 0.6\).

p <- 0.6

Then we say for now that the matrix of detections and non-detections is just \(z\) in which dead individuals are not detected.

y <- z
y[z==2] <- 0

Now for alive individuals, those for which \(y = z = 1\), we say they’re detected with probability 1.

y[y==1] <- rbinom(n = sum(y==1), 1, p)

Let’s get the number of individuals that have at least a detection.

nobs <- sum(apply(y,1,sum) != 0)
nobs
[1] 51

Before going any further, we need to get rid of the 6 individuals that have never been detected.

y <- y[apply(y,1,sum)!=0, ]

Let’s get the occasion of first capture for each individual.

first <- apply(y, 1, function(x) min(which(x != 0)))

For convenience, we will replace the 0s before first detection by NAs.

for (i in 1:nobs){
  if(first[i] > 1) y[i, 1:(first[i]-1)] <- NA
}
y %>%
  as_tibble() %>%
  add_column(id = 1:nobs, .before = "winter 1") %>%
  kableExtra::kable() %>%
  kableExtra::scroll_box(width = "100%", height = "300px")
id winter 1 winter 2 winter 3 winter 4 winter 5
1 1 1 0 0 1
2 1 0 1 1 0
3 NA NA 1 0 1
4 1 0 0 0 1
5 1 0 0 0 0
6 1 0 0 0 0
7 NA 1 0 0 0
8 NA 1 1 1 0
9 NA 1 1 0 0
10 1 0 1 1 0
11 1 1 1 0 0
12 1 0 0 0 0
13 NA 1 1 1 1
14 1 0 0 0 0
15 NA 1 1 1 0
16 1 1 0 0 0
17 1 0 0 0 0
18 NA 1 0 0 0
19 1 0 0 0 0
20 1 1 1 0 0
21 1 1 1 1 0
22 1 0 0 1 0
23 1 0 1 1 0
24 1 0 0 0 0
25 1 1 1 1 1
26 NA 1 1 0 0
27 NA 1 1 0 0
28 1 1 0 1 1
29 1 0 0 0 0
30 1 1 1 1 0
31 NA 1 1 1 1
32 NA NA 1 1 0
33 1 1 0 1 0
34 1 1 0 1 1
35 1 0 0 0 0
36 1 0 1 1 0
37 1 0 1 1 0
38 NA NA NA NA 1
39 NA 1 1 0 1
40 NA NA NA NA 1
41 1 0 1 0 0
42 1 0 0 0 0
43 1 1 1 1 0
44 1 1 1 0 1
45 1 1 0 0 0
46 1 1 1 1 1
47 1 0 0 0 0
48 NA 1 1 1 0
49 1 1 1 0 0
50 1 1 1 0 0
51 1 0 0 0 0

Now we’re ready to fit our HMM to the data we simulated. As usual, we first define the model.

hmm.survival <- 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]) # initial state prob
    for (j in (first[i]+1):T){
      z[i,j] ~ dcat(gamma[z[i,j-1], 1:2]) # z_t given z_(t-1)
      y[i,j] ~ dcat(omega[z[i,j], 1:2])   # y_t given z_t
    }
  }
})

Then put the constants in a list.

my.constants <- list(N = nrow(y),      # nb of individuals
                     T = 5,            # nb of occasions
                     first = first)    # occasions of first capture
my.constants
$N
[1] 51

$T
[1] 5

$first
 [1] 1 1 3 1 1 1 2 2 2 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 2 2 1 1 1 2 3 1 1 1 1 1 5
[39] 2 5 1 1 1 1 1 1 1 2 1 1 1

Now 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)

Specify the parameters we’d like to monitor.

parameters.to.save <- c("phi", "p")
parameters.to.save
[1] "phi" "p"  

MCMC details.

n.iter <- 5000
n.burnin <- 1000
n.chains <- 2

At last, let’s run nimble.

mcmc.output <- nimbleMCMC(code = hmm.survival,               # model code 
                          constants = my.constants,          # constants
                          data = my.data,                    # data
                          inits = initial.values,            # initial values
                          monitors = parameters.to.save,     # parameters to monitor
                          niter = n.iter,                    # nb of iterations
                          nburnin = n.burnin,                # length of the burn-in period
                          nchains = n.chains)                # nb of chains
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

We display the results. The posterior means of detection and survival are close to the parameters we used to simulate the data.

library(MCMCvis)
MCMCsummary(mcmc.output, round = 2)
    mean   sd 2.5%  50% 97.5% Rhat n.eff
p   0.63 0.06 0.51 0.63  0.75 1.01   613
phi 0.81 0.04 0.72 0.81  0.89 1.01   674