LaMa (version 2.0.5)

pseudo_res: Calculate pseudo-residuals

Description

For HMMs, pseudo-residuals are used to assess the goodness-of-fit of the model. These are based on the cumulative distribution function (CDF) $$F_{X_t}(x_t) = F(x_t \mid x_1, \dots, x_{t-1}, x_{t+1}, \dots, x_T)$$ and can be used to quantify whether an observation is extreme relative to its model-implied distribution.

This function calculates such residuals via probability integral transform, based on the local state probabilities obtained by stateprobs or stateprobs_g and the respective parametric family.

Usage

pseudo_res(
  obs,
  dist,
  par,
  stateprobs = NULL,
  mod = NULL,
  normal = TRUE,
  discrete = NULL,
  randomise = TRUE,
  seed = NULL
)

Value

vector of pseudo residuals

Arguments

obs

vector of continuous-valued observations (of length n)

dist

character string specifying which parametric CDF to use (e.g., "norm" for normal or "pois" for Poisson) or CDF function to evaluate directly. If a discrete CDF is passed, the discrete argument needs to be set to TRUE because this cannot determined automatically.

par

named parameter list for the parametric CDF

Names need to correspond to the parameter names in the specified distribution (e.g. list(mean = c(1,2), sd = c(1,1)) for a normal distribution and 2 states). This argument is as flexible as the parametric distribution allows. For example you can have a matrix of parameters with one row for each observation and one column for each state.

stateprobs

matrix of local state probabilities for each observation (of dimension c(n,N), where N is the number of states) as computed by stateprobs, stateprobs_g or stateprobs_p

mod

optional model object containing initial distribution delta, transition probability matrix Gamma, matrix of state-dependent probabilities allprobs, and potentially a trackID variable

If you are using automatic differentiation either with RTMB::MakeADFun or qreml and include forward, forward_g or forward_p in your likelihood function, the objects needed for state decoding are automatically reported after model fitting. Hence, you can pass the model object obtained from running RTMB::report() or from qreml directly to this function and avoid calculating local state proabilities manually. In this case, a call should look like pseudo_res(obs, "norm", par, mod = mod).

normal

logical, if TRUE, returns Gaussian pseudo residuals

These will be approximately standard normally distributed if the model is correct.

discrete

logical, if TRUE, computes discrete pseudo residuals (which slightly differ in their definition)

By default, will be determined using dist argument, but only works for standard discrete distributions. When used with a special discrete distribution, set to TRUE manually. See pseudo_res_discrete for details.

randomise

for discrete pseudo residuals only. Logical, if TRUE, return randomised pseudo residuals. Recommended for discrete observations.

seed

for discrete pseudo residuals only. Integer, seed for random number generation

Details

When used for discrete pseudo-residuals, this function is just a wrapper for pseudo_res_discrete.

Examples

## continuous-valued observations
obs = rnorm(100)
stateprobs = matrix(0.5, nrow = 100, ncol = 2)
par = list(mean = c(1,2), sd = c(1,1))
pres = pseudo_res(obs, "norm", par, stateprobs)

## discrete-valued observations
obs = rpois(100, lambda = 1)
par = list(lambda = c(1,2))
pres = pseudo_res(obs, "pois", par, stateprobs)

## custom CDF function
obs = rnbinom(100, size = 1, prob = 0.5)
par = list(size = c(0.5, 2), prob = c(0.4, 0.6))
pres = pseudo_res(obs, pnbinom, par, stateprobs, 
                  discrete = TRUE)
# if discrete CDF function is passed, 'discrete' needs to be set to TRUE

## full example with model object
step = trex$step[1:200]

nll = function(par){
  getAll(par)
  Gamma = tpm(logitGamma)
  delta = stationary(Gamma)
  mu = exp(logMu); REPORT(mu)
  sigma = exp(logSigma); REPORT(sigma)
  allprobs = matrix(1, length(step), 2)
  ind = which(!is.na(step))
  for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

par = list(logitGamma = c(-2,-2), 
           logMu = log(c(0.3, 2.5)), 
           logSigma = log(c(0.3, 0.5)))
           
obj = MakeADFun(nll, par)
opt = nlminb(obj$par, obj$fn, obj$gr)

mod = obj$report()

pres = pseudo_res(step, "gamma2", list(mean = mod$mu, sd = mod$sigma),
                  mod = mod)

qqnorm(pres)
abline(a = 0, b = 1)