4 Longitudinal data

Before diving into this vignette, we recommend reading the vignettes Introduction to LaMa and Finer control.

library(LaMa)
#> Loading required package: RTMB

In real-data applications, one will often be faced with a data set consisting of several measurement tracks that can reasonably be assumed to be mutually independent. Examples of such a longitudinal structure include GPS tracks of several individuals (or several tracks, e.g. days, of one individual), or when analysing sports data, time series for separate games. In such settings, the researcher must decide whether to pool parameters across tracks or not. Here, we provide brief examples for complete and partial pooling.

When the data consists of K independent tracks, the joint likelihood is L(\boldsymbol{\theta}) = \prod_{k=1}^K L_k(\boldsymbol{\theta}), where L_k(\boldsymbol{\theta}) is the usual HMM likelihood for the k-th track, so the log-likelihood is a sum over K tracks. In LaMa, this is handled efficiently by passing a trackID vector to forward(), which specifies the track index of each observation. The forward algorithm then resets to the initial distribution at the start of each new track and sums up the individual log-likelihood contributions automatically. Importantly, always using the trackID argument rather than calling forward() separately for each track ensures that the automatic reporting mechanism inside forward() — on which viterbi(), stateprobs(), and other downstream functions rely — works correctly.

Complete pooling

With complete pooling, all parameters are shared across individuals.

Writing the negative log-likelihood function

The likelihood function is almost identical to the single-track case — the only difference is passing trackID as an additional argument to forward(). With delta = "stationary", each track starts in the stationary distribution.

nll_pool = function(par) {
  getAll(par, dat)
  Gamma = tpm(eta)
  mu = exp(log_mu); REPORT(mu)
  sigma = exp(log_sigma); REPORT(sigma)
  allprobs = statedens(x, dnorm, mu, sigma)
  -forward("stationary", Gamma, allprobs, trackID)
}

Simulating data

We can use the same likelihood function to simulate data. lama() does not need a fitted model: called on an object from MakeADFun() before optimisation, it reports at the parameter values the object was built with – here the true parameters. simulate() then draws states and observations from this model on the same design, i.e. with the tracks given by trackID. As there are no data yet, we pass an all-NA placeholder for x, which tells simulate() to fill in the whole stream. We simulate K = 200 tracks of length n = 50:

K = 200 # number of individuals
n = 50  # observations per individual

# true parameters, on the same scale as in the likelihood
par_true = list(
  eta = c(-3, -2),           # t.p.m. parameters (logit scale)
  log_mu = log(c(15, 60)),   # state-dependent means (log scale)
  log_sigma = log(c(10, 40)) # state-dependent sds (log scale)
)
dat = list(
  x = rep(NA, n * K),        # placeholder for the observations
  trackID = rep(1:K, each = n)
)

sim = simulate(lama(MakeADFun(nll_pool, par_true, silent = TRUE)), seed = 123)
head(sim)
#>   trackID state         x
#> 1       1     2 120.94832
#> 2       1     2  32.25963
#> 3       1     2 120.33054
#> 4       1     2  58.43673
#> 5       1     2  83.11626
#> 6       1     2  65.91779
dat$x = sim$x

Estimating the model

par = list(
  eta = rep(-2, 2),          # initial t.p.m. parameters (logit scale)
  log_mu = log(c(20, 50)),   # initial state-dependent means (log scale)
  log_sigma = log(c(10, 30)) # initial state-dependent sds (log scale)
)

obj_pool = MakeADFun(nll_pool, par, silent = TRUE)
system.time(
  opt_pool <- nlminb(obj_pool$par, obj_pool$fn, obj_pool$gr)
)
#>    user  system elapsed 
#>   0.204   0.000   0.204
mod_pool = lama(obj_pool)

The estimated parameters should be close to the true values used in the simulation:

mod_pool$mu     # true: 15, 60
#> [1] 15.04398 59.97678
mod_pool$sigma  # true: 10, 40
#> [1]  9.994448 40.204912
mod_pool$Gamma
#>            S1        S2
#> S1 0.89197388 0.1080261
#> S2 0.04650452 0.9534955
tpm(par_true$eta) # true t.p.m.
#>            S1        S2
#> S1 0.88079708 0.1192029
#> S2 0.04742587 0.9525741

Partial pooling

If some parameters are individual-specific while the rest are shared, we speak of partial pooling. We demonstrate this for K = 5 individuals whose transition probability matrices depend on an individual-level covariate (e.g. age), while the state-dependent parameters remain shared. We estimate the effect of this covariate on the transition probabilities.

Writing the negative log-likelihood function

Each individual has its own constant transition probability matrix, determined by its covariate value. tpm() computes all K individual-specific matrices at once, and forward() uses matrix k for track k. With delta = "stationary", each track starts in the stationary distribution of its own matrix.

nll_partial = function(par) {
  getAll(par, dat)
  Gamma = tpm(beta, matrix(z, ncol = 1)); REPORT(beta) # array of K individual-specific tpms
  mu = exp(log_mu); REPORT(mu)
  sigma = exp(log_sigma); REPORT(sigma)
  allprobs = statedens(x, dnorm, mu, sigma)
  -forward("stationary", Gamma, allprobs, trackID)
}

Simulating data

K = 5   # number of individuals
n = 200 # observations per individual

set.seed(123)
z = rnorm(K) # individual-level covariate

par_true = list(
  beta = matrix(c(-2, -2, 1, -1), nrow = 2), # intercepts and slopes of the tpm regression
  log_mu = log(c(15, 60)),
  log_sigma = log(c(10, 40))
)
dat = list(
  x = rep(NA, n * K),
  z = z,
  trackID = rep(1:K, each = n)
)

sim = simulate(lama(MakeADFun(nll_partial, par_true, silent = TRUE)), seed = 123)
dat$x = sim$x

Estimating the model

par = list(
  beta = matrix(c(-2, -2, 0, 0), nrow = 2), # regression coefficients for tpm
  log_mu = log(c(20, 50)),                  # initial state-dependent means
  log_sigma = log(c(10, 30))                # initial state-dependent sds
)

obj_partial = MakeADFun(nll_partial, par, silent = TRUE)
system.time(
  opt_partial <- nlminb(obj_partial$par, obj_partial$fn, obj_partial$gr)
)
#>    user  system elapsed 
#>   0.017   0.000   0.017
mod_partial = lama(obj_partial)

The estimated regression coefficients for the transition probabilities and the shared state-dependent parameters:

mod_partial$beta   # true: c(-2, -2, 1, -1)
#>           [,1]       [,2]
#> [1,] -1.875389  0.6667121
#> [2,] -1.823894 -0.9405922
mod_partial$mu     # true: 15, 60
#> [1] 15.17265 61.42479
mod_partial$sigma  # true: 10, 40
#> [1] 10.29714 38.78404

Random effects

In the partial pooling example above, we estimated a separate transition probability matrix for each individual based on a fixed covariate. An alternative and often more parsimonious approach is to treat the individual-level variation in transition probabilities as arising from Gaussian random effects. Instead of estimating K independent transition structures, we model each individual’s linear predictor as a deviation from a shared population-level intercept: \text{logit}(\gamma_{ij}^{(k)}) = \beta_0^{(ij)} + u_{ij}^{(k)}, \quad u_{ij}^{(k)} \overset{\text{iid}}{\sim} N\bigl(0, (\sigma_u^{(ij)})^2\bigr), for each off-diagonal element (i \neq j) and each individual k = 1, \dots, K. The population-level intercept \beta_0^{(ij)} captures the average switching tendency, while \sigma_u^{(ij)} quantifies how much individuals vary around it. The number of state-process parameters is 2 \times N(N-1) — two intercepts and two standard deviations for a two-state model — regardless of how many individuals are in the data.

To fit this model we need to integrate out the random effects from the joint likelihood. RTMB handles this automatically when the random effects are declared via the random argument of MakeADFun().

Writing the joint negative log-likelihood

The joint negative log-likelihood includes both the data likelihood (via forward()) and the Gaussian density of the random effects. The random effects are stored in a matrix u with one row per individual, so u + rep(beta_0, each = K) is the matrix of linear predictors, from which tpm(Eta = ...) computes one transition probability matrix per row.

jnll_re = function(par) {
  getAll(par, dat)
  # random effects: u[,ij] ~ N(0, sigma_u[ij]^2)
  sigma_u = exp(log_sigma_u); REPORT(sigma_u); ADREPORT(sigma_u)
  jnll = -sum(dnorm(u, 0, rep(sigma_u, each = K), log = TRUE))
  # individual-specific tpms from population intercept + random effect
  Gamma = tpm(Eta = u + rep(beta_0, each = K))
  # observation model
  mu = exp(log_mu); REPORT(mu)
  sigma = exp(log_sigma); REPORT(sigma)
  allprobs = statedens(x, dnorm, mu, sigma)
  jnll - forward("stationary", Gamma, allprobs, trackID)
}

Simulating data

For models with random effects, simulate() can work in two ways:

By default, it simulates conditional on the reported random effects, i.e. for the same individuals. For a fitted model, these are the estimated random effects. Here, they would be the zeros in par_true, so all individuals would have the same transition probabilities.
With random = TRUE, it first draws new random effects from their distribution and then the states and observations given these draws, i.e. for new individuals. This requires that the density of u is evaluated in the likelihood before u is used, as in jnll_re.

To generate data from the random effects model, we hence need random = TRUE:

K = 50  # number of individuals
n = 100 # observations per individual

par_true = list(
  beta_0 = c(-2, -2),         # population-level intercepts
  log_sigma_u = log(c(1, 1)), # log-SDs of random effects
  u = matrix(0, K, 2),        # random effects, drawn by simulate()
  log_mu = log(c(15, 60)),
  log_sigma = log(c(10, 40))
)
dat = list(
  x = rep(NA, n * K),
  trackID = rep(1:K, each = n),
  K = K
)

obj_true = MakeADFun(jnll_re, par_true, random = "u", silent = TRUE)
sim = simulate(lama(obj_true), seed = 123, random = TRUE)
dat$x = sim$x

Estimating the model

The random effects matrix u is declared via random = "u", telling RTMB to integrate it out. The optimiser only sees the marginal likelihood as a function of the fixed-effect parameters beta_0, log_sigma_u, log_mu, and log_sigma.

par = list(
  beta_0 = rep(-2, 2),        # population-level intercepts
  log_sigma_u = log(c(1, 1)), # log-SDs of random effects
  u = matrix(0, K, 2),        # random effects (integrated out)
  log_mu = log(c(20, 50)),    # initial state-dependent means
  log_sigma = log(c(10, 30))  # initial state-dependent sds
)

obj_re = MakeADFun(jnll_re, par, random = "u", silent = TRUE)
system.time(
  opt_re <- nlminb(obj_re$par, obj_re$fn, obj_re$gr)
)
#>    user  system elapsed 
#>   3.866   0.072   3.939
mod_re = lama(obj_re)

The estimated standard deviations of the random effects tell us how much the transition probabilities vary across individuals. We can check the key population-level estimates against the true values used in the simulation:

mod_re$sigma_u  # true: 1, 1
#> [1] 1.0921681 0.8144595
mod_re$mu       # true: 15, 60
#> [1] 14.72053 60.42473
mod_re$sigma    # true: 10, 40
#> [1]  9.798952 39.633995

Standard errors for all fixed effects and ADREPORT()ed quantities are available via sdreport():

sdr_re = sdreport(obj_re)
summary(sdr_re, "report") # sigma_u with SEs
#>          Estimate Std. Error
#> sigma_u 1.0921681  0.1592198
#> sigma_u 0.8144595  0.1399563

For the fitted model, the two ways of simulating answer different questions. simulate(mod_re) simulates new data for the 50 individuals in the data, using their estimated random effects, while simulate(mod_re, random = TRUE) simulates data for 50 new individuals from the estimated population. Comparing the empirical switching rates per individual shows the difference:

switch_rate = function(sim) tapply(sim$state, sim$trackID, function(s) mean(diff(s) != 0))

rates_same = switch_rate(simulate(mod_re, seed = 1))
rates_new = switch_rate(simulate(mod_re, seed = 1, random = TRUE))
cor(rates_same, switch_rate(sim)) # same individuals: related to the original rates
#> [1] 0.8934546
cor(rates_new, switch_rate(sim))  # new individuals: unrelated
#> [1] -0.02842753

Continue reading with Penalised splines.

Jan-Ole Fischer

Complete pooling

Writing the negative log-likelihood function

Simulating data

Estimating the model

Partial pooling

Writing the negative log-likelihood function

Simulating data

Estimating the model

Random effects

Writing the joint negative log-likelihood

Simulating data

Estimating the model