LaMa (version 3.1.0)

lama: Build a LaMaModel from an optimised RTMB object

Description

Having fitted a latent Markov model using automatic differentiation via RTMB, this function calls RTMB's obj$report() and does some additional processing. This then yields estimated parameters on their natural scale, allows for convenient calculation of AIC and BIC, state-decoding, and residual calculation. print() gives a short overview of the model, summary() the estimates. report() is the deprecated old name of lama().

Usage

lama(obj)
report(obj)

Value

A model object of class "LaMaModel" containing a list with the reported quantities from the RTMB object, along estimated parameters (par), the log-likelihood (llk), number of parameters (df) and observations (nobs), and the maximum gradient component (mgc) at the estimate. It also contains obj itself (not a copy), so that functions like MCreport or scores can be called on the model object.

Arguments

obj

Optimised RTMB object

Details

The number of observations (used by BIC) is the number of rows of the reported allprobs that are not all 1, i.e. missing observations are not counted. This assumes that forward is called only once in the likelihood function.

Examples

data <- trex[1:200,]

# initial parameters and observations
par <- list(
  log_mu = log(c(0.3, 1)),      # initial means for step length (log-transformed)
  log_sigma = log(c(0.2, 0.7)), # initial sds for step length (log-transformed)
  eta = rep(-2, 2)              # initial t.p.m. parameters (on logit scale)
)    
dat <- list(
  step = data$step,             # hourly step lengths
  nStates = 2                   # number of hidden states
)

# likelihood function
nll <- function(par) {
  getAll(par, dat)
  Gamma <- tpm(eta)
  delta <- stationary(Gamma)
  mu <- exp(log_mu); REPORT(mu)
  sigma <- exp(log_sigma); REPORT(sigma)
  allprobs <- matrix(1, length(step), nStates)
  ind <- which(!is.na(step))
  for(j in 1:nStates) {
    allprobs[ind,j] <- dgamma2(step[ind], mu[j], sigma[j]) 
  } 
  -forward(delta, Gamma, allprobs)
}

# automatic differentiation and optimisation
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)

### reporting ###
mod <- lama(obj)
mod          # short overview
summary(mod) # estimated state process and reported quantities

# estimated parameters
mod$par

# estimated quantities on natural scale
mod$mu
mod$sigma
mod$Gamma

# information criteria
AIC(mod)
BIC(mod)

# state decoding
states <- viterbi(mod = mod)   # global decoding
probs <- stateprobs(mod = mod) # local decoding

# residual calculation
pres <- pseudo_res(data$step, # observation sequence
                   "gamma2",  # distribution family
                   list(mean = mod$mu, sd = mod$sigma), # parameters for that family
                   mod = mod) # model object

# or from the model object alone, and plotted
plot(mod)