The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

EM + Laplace for latent-variable models

library(tulpa)

What the driver is for

Some models have a latent variable that, once known, turns the fit into an ordinary weighted GLM: occupancy (was the site occupied?), N-mixture (what was the true count?), zero-inflation (is this a structural zero?). tulpa_em_laplace() is the generic engine for that pattern. It does not know your model – you supply two callbacks and it runs the EM loop, solving each M-step submodel with tulpa_laplace() and checking convergence.

This is deliberately an engine door: model packages such as tulpaObs build their occupancy and N-mixture fitters on top of it. The callbacks below are the whole contract.

The two callbacks

e_step(fits, ...) receives the current per-submodel fits and returns the latent-variable posterior as weights:

e_step <- function(fits, ...) {
  # Use the current fits to compute the responsibility of each observation
  # (e.g. the posterior probability that a zero is a Poisson zero, not a
  # structural one). Return them as a `weights` element.
  list(weights = responsibilities)
}

m_step_encode(weights, ...) turns those weights into one or more weighted GLM submodels. Each block is a list(y, n_trials, X, family, offset); the driver fits each with tulpa_laplace() and threads the per-block family and offset through automatically:

m_step_encode <- function(weights, ...) {
  list(
    lambda = list(y = counts, n_trials = 1L, X = X_abund,
                  family = "poisson", offset = log(weights)),
    pi     = list(y = z,      n_trials = 1L, X = X_zero,
                  family = "binomial", offset = NULL)
  )
}

Running it, and the bias correction

fit <- tulpa_em_laplace(
  e_step        = e_step,
  m_step_encode = m_step_encode,
  max_iter      = 30L,
  tol           = 1e-4
)
fit$fits        # the converged per-submodel Laplace fits
fit$n_iter
fit$converged

The EM point estimate is fast but its standard errors ignore the uncertainty in the imputed latent variable. Two post-EM corrections restore it, pooled by Rubin’s rules:

# Multiple imputation: refit on several latent draws and pool.
fit_mi <- tulpa_em_laplace(e_step, m_step_encode, correction = "mi")

# Warm-started Gibbs from the EM mode, then pool.
fit_gibbs <- tulpa_em_laplace(e_step, m_step_encode, correction = "gibbs")

An optional beta_prior = list(mean, sd) threads a Gaussian fixed-effect prior into every M-step block and into the correction refits, and m_step_extra(fits, weights, ...) updates non-linear-predictor parameters (dispersions, mixing weights) between the M- and E-steps. See ?tulpa_em_laplace. For a Monte-Carlo E-step, use tulpa_em_mc().

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.