TheoMichelot / hmmTMB

Fit hidden Markov models using Template Model Builder (TMB): flexible state-dependent distributions, transition probability structures, random effects, and smoothing splines.
53 stars 7 forks source link

Error in HMM$predict for tpm with covariates #8

Closed TheoMichelot closed 3 years ago

TheoMichelot commented 3 years ago

An error arises when I try to plot the effect of a covariate on the transition probabilities (with traceback):

Error in self[[comp]]()$update_coeff_fe(new_par) : 
  'coeff_fe' should be of length 4 (one parameter for each column of the design matrix) 
7. stop("'coeff_fe' should be of length ", ncol_total, " (one ", 
    "parameter for each column of the design matrix)") at markovchain.R#221
6. self[[comp]]()$update_coeff_fe(new_par) at hmm.R#1162
5. func(x, ...) 
4. jacobian.default(lfn, par, oldpar = oldpar, ind_fe = ind_fe, 
    ind_re = ind_re, comp = comp) 
3. numDeriv::jacobian(lfn, par, oldpar = oldpar, ind_fe = ind_fe, 
    ind_re = ind_re, comp = comp) at hmm.R#1177
2. self$predict(name, t = "all", newdata = newdata, level = 0.95) at hmm.R#1428
1. hmm$plot("tpm", "tod") 

The problem seems to be in the delta method in the HMM$predict function. It doesn't come up when predicting covariate effects on the state-dependent observation parameters. I include a complete reproducible example below, based on some dummy data. The same issue happens in the examples from the vignette.

library(hmmTMB)
set.seed(56834)

# Dummy data
n <- 500
data <- data.frame(z = rnorm(n),
                   tod = (1:n)%%24)

# Create MarkovChain
f <- ~ tod
hid <- MarkovChain$new(n_states = 2, structure = f, data = data)

# Create Observation model
dist <- list(z = "norm")
par0 <- list(z = list(mean = c(-0.5, 0.5), sd = c(0.6, 0.6)))
formulas <- list(z = list(mean = f, sd = ~1))
obs <- Observation$new(data = data, dists = dist, n_states = 2, 
                       par = par0, formulas = formulas)

# Fit HMM
hmm <- HMM$new(obs = obs, hidden = hid)
hmm$fit(silent = FALSE)

# This works
hmm$plot("obspar", "tod")

# This fails
hmm$plot("tpm", "tod")
r-glennie commented 3 years ago

This was caused by a naming issue in $predict where "hidden" was named as "hid". Fixed.