LaMa (version 2.1.2)

report: Get reported quantities from and RTMB object and return a LaMaModel

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.

Usage

report(obj)

Value

A model object of class "LaMaModel" containing a list with the reported quantities from the RTMB object, along estimated parameters and other quantities.

Arguments

obj

Optimised RTMB object

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 <- report(obj)

# 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