1.1 Conditional response distributions

Let \(y_i\) denote observation \(i\), and let \(u\) collect all random effects. Conditional on \(u\), the observations are modeled as independent, with conditional means \(\mu_i\) and a specified exponential-family or other supported response distribution:

\[ y_i \mid u \sim f_i(y_i \mid \mu_i, \phi_i), \qquad g_i(\mu_i) = \eta_i = o_i + x_i^\mathsf{T}\beta + z_i^\mathsf{T}u. \]

Here \(g_i\) is the link, \(o_i\) is an optional offset, \(x_i\) and \(z_i\) are design rows, and \(\beta\) contains fixed effects. The conditional mean and variance are functions of the family parameters:

\[ \operatorname{E}(y_i\mid u)=\mu_i,\qquad \operatorname{Var}(y_i\mid u)=\phi_i V_i(\mu_i). \]

For Bernoulli and Poisson families, dispersion is fixed by the distribution. For other families, dispersion or shape parameters may be estimated or fixed where the selected family and fitting method support them.

In everyday terms, \(x_i^\mathsf{T}\beta\) is the part of the prediction made from measured explanatory variables. The term \(z_i^\mathsf{T}u\) adds a group-specific adjustment, and \(o_i\) adds a known quantity such as log exposure. Their sum is called the linear predictor, \(\eta_i\). It is linear because the effects are added together there, even though the prediction on the outcome scale may not be a straight line.

The link \(g\) translates between the outcome’s mean and this linear predictor. For a yes/no outcome, the logit link turns a probability between 0 and 1 into log-odds, which can range over all real numbers. For counts, a log link turns a positive mean into its logarithm. This lets us use an ordinary additive predictor without accidentally predicting a probability below zero or a negative count mean.

The variance formula says that different outcomes have different natural amounts of noise. A Bernoulli outcome with probability \(p\) has variance \(p(1-p)\): it is most uncertain near \(p=0.5\) and less uncertain near 0 or 1. A Poisson count with mean \(\lambda\) has variance \(\lambda\), so larger expected counts naturally have more spread.

The random effects follow a zero-mean multivariate normal distribution,

\[ u \sim N(0,G(\theta)), \]

where \(\theta\) denotes covariance parameters. In sommer, the covariance matrix is assembled from the random-effect terms, their covariance descriptors, and any supplied relationship matrices. Thus the fixed-effect formula describes the conditional linear predictor, while random and its covariance structures describe the latent Gaussian distribution.

The normal distribution here is a convenient model for the collection of group adjustments: zero is the average adjustment, while \(G(\theta)\) says how much the adjustments vary and which ones tend to move together. “Zero mean” does not mean every group effect is zero; it means positive and negative group deviations balance around zero across the population of groups.

1.2 Conditional and marginal quantities

The model is specified conditionally on \(u\). Its inverse-link mean \(g^{-1}(x_i^\mathsf{T}\beta+z_i^\mathsf{T}u+o_i)\) is generally not the mean after integrating over the random effects:

\[ \operatorname{E}(y_i\mid x_i) = \operatorname{E}_{u}\!\left[ g^{-1}(o_i+x_i^\mathsf{T}\beta+z_i^\mathsf{T}u) \right]. \]

For a nonlinear inverse link, this integral usually differs from \(g^{-1}(o_i+x_i^\mathsf{T}\beta)\), the response mean at zero random effects. For example, with a log link and a Gaussian random intercept \(u\sim N(0,\sigma_u^2)\),

\[ \operatorname{E}_{u}\{\exp(\eta+u)\} =\exp(\eta+\sigma_u^2/2), \]

whereas setting \(u=0\) gives \(\exp(\eta)\). This distinction is essential when describing predictions: a conditional prediction, a zero-random-effect prediction, and a population-integrated prediction answer different questions.

Why do those two numbers differ? The exponential curve bends upward, so a positive deviation increases the value more than an equally sized negative deviation decreases it. Averaging the values after applying the curve is therefore not the same as applying the curve to the average deviation. In calculus language, the inverse link is nonlinear; in plain language, “take the average, then transform” and “transform, then take the average” are different instructions.

The log-link example follows a familiar normal-distribution identity. If \(u\) is normal with variance \(\sigma_u^2\), then \(E\{\exp(u)\}=\exp(\sigma_u^2/2)\). Multiplying by \(\exp(\eta)\) gives \(\exp(\eta+\sigma_u^2/2)\). The extra term comes from averaging across groups; it is not a change to the conditional model.

To see where the identity comes from, write the expectation as an integral of \(\exp(u)\) times the normal density. Completing the square in the exponent turns that product into a shifted normal density times the constant \(\exp(\sigma_u^2/2)\). The shifted normal density still integrates to one, so only that constant remains. This is the same “complete the square” algebra used when working with quadratic expressions.

1.3 A concrete binary example

For a Bernoulli-logit model with subject intercepts,

\[ y_{ij}\mid u_j\sim \operatorname{Bernoulli}(p_{ij}),\qquad \operatorname{logit}(p_{ij}) =\beta_0+\beta_1 x_{ij}+u_j,\qquad u_j\sim N(0,\sigma_u^2). \]

The coefficient \(\beta_1\) is a conditional log-odds contrast, holding the latent subject effect fixed. It is not, in general, the log odds ratio obtained after averaging predicted probabilities over subjects. The random intercept induces within-subject dependence even though observations are independent conditional on that intercept.

For example, suppose two people have the same measured characteristics and the same subject effect. The model gives them the same success probability. If one subject’s hidden adjustment is higher, all observations from that subject tend to have higher probabilities together. This shared shift is how the model represents repeated measurements that resemble one another.

2.1 Covariance of coefficient blocks

A random term can have one or more coefficients per level. For example, a random intercept-and-slope model has a two-dimensional covariance matrix \(\Sigma\) for the intercept and slope. If the levels also have a known relationship matrix \(A\), the coefficient covariance has a Kronecker-product structure, with ordering determined by the model design:

\[ \operatorname{Var}(u) \ \text{has the structure}\ \Sigma\otimes A. \]

Equivalently, its precision is assembled from the inverse descriptor and the relationship precision. An identity relationship gives independent levels; a pedigree or genomic relationship matrix couples their effects. The covariance descriptor controls covariance among coefficient dimensions, while the relationship matrix controls covariance among levels.

For a less abstract picture, imagine three sites, each with a random intercept and slope. The small matrix \(\Sigma\) describes whether sites with higher intercepts also tend to have steeper slopes. The matrix \(A\) describes whether the sites themselves are related. Their Kronecker product combines these two questions into one covariance matrix for all six coefficients. In short: \(\Sigma\) is about “which effects within a level travel together?” and \(A\) is about “which levels are related?”

This separation is one of sommer’s key modeling capabilities. Structures such as unstructured, diagonal, compound-symmetric, autoregressive, factor analytic, and reduced-rank forms constrain \(\Sigma\) in different ways. Random slopes, multi-trait effects, genotype-by-environment terms, and relationship-informed genetic effects can therefore be expressed using the same underlying idea as in Gaussian mixed models. The existence of a covariance structure in the package does not imply that every structure is eligible in every GLMM engine; consult the GLMM guide and inspect the execution plan for the selected route.

2.2 Conditional independence is not marginal independence

For the likelihood methods described below, observation densities factor conditional on \(u\):

\[ p(y\mid u,\beta,\theta)=\prod_i p(y_i\mid u,\beta). \]

After integrating over \(u\), this product no longer factors across observations that share random effects. A random intercept therefore captures dependence through the shared latent variable; it is not an extra residual correlation term added independently to each response distribution.

Conditional independence means: if we already knew each group’s hidden adjustment, observations would be independent under the model. But we do not know those adjustments. When we average over them, two observations from the same group tend to move together because both share the same unknown value. That is why “independent given the random effect” does not mean “independent overall.”

PQL’s working Gaussian problem can use rcov and other supported working covariance structures. Native Laplace instead uses the conditional observation likelihood and the random-effect covariance. It currently does not accept structured residual covariance or the supplied observation precision matrix W.

2.3 Formulae and covariance specification

The syntax for fixed and random effects follows the usual mmes() conventions. For example:

fit <- mmes(
  count ~ treatment + offset(log(exposure)),
  random = ~ site,
  data = dat,
  family = poisson(),
  control = mmes_control(
    task = list(objective = "ml",
      glmm = list(method = "laplace"))
  )
)

The offset is known and enters the linear predictor with coefficient fixed at one. In a log-rate model, offset(log(exposure)) makes the expected count proportional to exposure while the fixed and random effects describe the log-rate. For example, doubling the exposure adds \(\log(2)\) to the predictor; after the log link is undone, that doubles the expected count while leaving the modeled rate unchanged.

3.1 Negative binomial and dispersion parameters

The NB2 parameter \(\theta\) is a positive shape parameter. Large \(\theta\) approaches the Poisson variance; smaller \(\theta\) allows greater overdispersion. On the native Laplace route, it may be estimated jointly or fixed with task$glmm$familyParameters. Gaussian dispersion is residual variance \(\phi\). Gamma dispersion satisfies \(\operatorname{Var}(Y\mid u)=\phi\mu^2\), so the Gamma shape is \(1/\phi\). Family-specific parameter estimates are reported on their natural scale.

Dispersion should not be used as a substitute for omitted structure. Overdispersion may arise from latent heterogeneity, excess zeros, dependence, or an incorrect mean model, and those mechanisms imply different interpretations and predictions.

The phrase “overdispersion” just means that the data vary more than the basic model expects. It describes a symptom, not its cause. For example, mixing together high-rate and low-rate groups can create extra variation, but adding an NB2 dispersion parameter does not automatically explain why those groups differ.

3.2 Zero-inflated Poisson

The current M5 implementation is a restricted zero-inflated Poisson (ZIP) model. Let \(S_i\) indicate membership in a structural-zero state:

\[ S_i\sim\operatorname{Bernoulli}(\pi),\qquad Y_i = \begin{cases} 0, & S_i=1,\\ \operatorname{Poisson}(\lambda_i), & S_i=0, \end{cases} \qquad \log(\lambda_i)=\eta_i. \]

Consequently,

\[ \Pr(Y_i=0)=\pi+(1-\pi)e^{-\lambda_i},\quad \Pr(Y_i=y>0)=(1-\pi)e^{-\lambda_i}\frac{\lambda_i^y}{y!}, \]

and

\[ \operatorname{E}(Y_i\mid u)=(1-\pi)\lambda_i,\qquad \operatorname{Var}(Y_i\mid u) =(1-\pi)\lambda_i+(1-\pi)\pi\lambda_i^2. \]

Think of ZIP as a two-step chance experiment. First, with probability \(\pi\), the observation is in a state that can only produce zero. Otherwise, with probability \(1-\pi\), we draw an ordinary Poisson count with mean \(\lambda_i\). A zero can therefore come from either step; a positive count can only come from the Poisson step.

The formulas follow from those two possibilities. A zero occurs either in the structural-zero state or when the Poisson draw happens to be zero, so the two probabilities add. The average count is \((1-\pi)\lambda_i\) because the structural-zero state contributes zero and only the other state contributes a Poisson mean. For the variance, use the law of total variance: average the within-state variance, then add the variation caused by switching between the two states. In symbols, \(\operatorname{Var}(Y)=E\{\operatorname{Var}(Y\mid S)\} +\operatorname{Var}\{E(Y\mid S)\}\). The first part is \((1-\pi)\lambda_i\); the second is \(\pi(1-\pi)\lambda_i^2\). Adding them gives the formula above.

zipoisson() estimates one common intercept-only structural-zero probability by default. The count component’s linear predictor models \(\log(\lambda_i)\); the fitted response mean is \((1-\pi)\lambda_i\). The probability can be fixed or estimated using task$glmm$familyParameters. A separate fixed-effects formula, zipoisson(zero=~covariates), extends this to \(\operatorname{logit}(\pi_i)=a_i+w_i^\mathsf{T}\gamma\), with its own offset \(a_i\) and fixed design \(w_i\). Its logit coefficients are jointly optimized with the count coefficients and covariance parameters and appear in zeroCoefficients; summary() computes their SEs only after fitting. The mixture mean and variance above then use \(\pi_i\) row by row. The default zero=~1 preserves the scalar family-parameter interface. Zero-component random effects remain unsupported. ZIP is not a hurdle model: zeros may arise from either component, and positive counts follow a zero-truncated Poisson distribution conditional on positivity.

For a zero formula, \(\exp(\gamma_j)\) is a structural-zero odds ratio conditional on the other zero covariates, not a count rate ratio. A covariate present in both components can change the observed mean through both predictors:

\[ \frac{\partial\log\{\operatorname{E}(Y_i\mid u)\}}{\partial x} =\beta_x-\pi_i\gamma_x \]

when \(x\) enters both linear predictors without transformations or interactions. Count and zero coefficients must therefore not be interpreted as two independent regressions on observed positives and zeros.

Here is the algebra behind that derivative. The mean is \((1-\pi_i)\lambda_i\). Taking logs turns the product into a sum: \(\log(1-\pi_i)+\log(\lambda_i)\). A one-unit increase in \(x\) changes the second term by \(\beta_x\). For the logit model of \(\pi_i\), the derivative of \(\log(1-\pi_i)\) is \(-\pi_i\gamma_x\). Adding the two changes gives \(\beta_x-\pi_i\gamma_x\). So a count-rate increase can be partly offset if the same covariate also makes a structural zero more likely.

To see why, let \(\tau_i=\Pr(S_i=1\mid Y_i,u)\) under the fitted model. For a positive observation \(\tau_i=0\); for a zero,

\[ \tau_i=\frac{\pi_i}{\pi_i+(1-\pi_i)\exp(-\lambda_i)},\qquad \frac{\partial\log p(Y_i\mid u)}{\partial\gamma} =w_i(\tau_i-\pi_i). \]

The quantity \(\tau_i\) is the model’s estimated chance that a particular observed zero came from the zero-only state. It is not necessarily 0 or 1: the data usually cannot tell which process produced that zero. The score formula says that the zero-component coefficients are updated by comparing this observation-specific chance with the model’s overall probability \(\pi_i\), weighted by its predictors \(w_i\).

An observed zero does not reveal which component generated it. The native optimizer differentiates the complete marginal Laplace approximation, including the mode and determinant corrections, rather than merely summing these conditional scores. Zero coefficients enter the same joint outer parameter vector as count coefficients and covariance parameters. There is no second latent random-effect block for the zero predictor.

The null case \(\pi=0\) is a boundary of the parameter space. Standard interior-point Wald approximations are not automatically valid there, and near-zero inflation can be weakly identified from the Poisson count mean. Model convergence alone does not establish that the ZIP component is identified. Use simulation-based calibration, boundary diagnostics, and scientific justification before interpreting \(\pi\).

3.3 Mixed families and long data

familym() assigns a response family to each trait label in a long-format data set. Conditional on the shared random effects, each row contributes its own family density. Random-effect structures may link traits, but the likelihood still uses each row’s declared observation family and its appropriate trials or family parameters.

In long format, each row is one measurement, and a column says which trait was measured. For instance, disease rows can use a binomial model while count rows use a negative binomial model. Shared random effects can connect the traits, but each row still needs the probability model that matches its own outcome.

fit <- mmes(
  value ~ trait - 1 + trait:environment,
  random = ~ vsm(usm(trait), ism(genotype)),
  data = long_data,
  family = familym(
    "trait",
    disease = binomial("probit"),
    count = nbinom2(theta = 2),
    biomass = Gamma("log"),
    yield = gaussian()
  ),
  control = mmes_control(
    task = list(objective = "ml",
      glmm = list(method = "laplace",
        familyParameters = list(
          count = list(theta = 2, estimate = TRUE),
          biomass = list(dispersion = 0.4, estimate = TRUE),
          yield = list(dispersion = 1, estimate = FALSE)
        )))
  )
)

The family-specific likelihood is useful for heterogeneous traits, but it does not make the residual covariance structure of a Gaussian multivariate model available to the native Laplace likelihood. Model the cross-trait dependence through eligible shared random effects.

4.1 The working response

PQL replaces the nonlinear response model, locally at the current predictor, with a weighted Gaussian working problem. At outer iteration \(t\), define

\[ \mu_i^{(t)}=g^{-1}(\eta_i^{(t)}),\qquad h_i^{(t)}=\frac{d\mu_i}{d\eta_i}\bigg|_{\eta_i^{(t)}},\qquad V_i^{(t)}=V_i(\mu_i^{(t)}). \]

With prior precision weight \(w_i\) (for example, binomial trial count), the working response and working precision are

\[ z_i^{(t)} =\eta_i^{(t)}+\frac{y_i-\mu_i^{(t)}}{h_i^{(t)}}-o_i, \qquad d_i^{(t)}=w_i\frac{(h_i^{(t)})^2}{V_i^{(t)}}. \]

These symbols look intimidating, but they describe a familiar strategy: replace a curved problem by a straight-line approximation near the current guess. The derivative \(h_i^{(t)}\) answers, “If the predictor changed a little here, how much would the mean change?” A large \(h_i^{(t)}\) means the mean responds strongly; a small one means it barely moves.

Remember from algebra that a straight line through a point can be written as its current value plus slope times the change in input. Applying that idea to the mean gives the local approximation \(\mu_i\approx\mu_i^{(t)}+h_i^{(t)}(\eta_i-\eta_i^{(t)})\). Solving this approximation for the predictor that would match the observed response \(y_i\) gives \(\eta_i^{(t)}+(y_i-\mu_i^{(t)})/h_i^{(t)}\). This is the “working response”; the known offset is subtracted because the Gaussian working fit estimates the rest of the predictor, not the offset.

The working precision \(d_i^{(t)}\) says how strongly each working response should count in that fit. Responses whose means are locally sensitive to the predictor get more precision, while responses with naturally high variance get less. This is the same idea as weighted averages: more precise measurements should pull the answer more.

The offset is removed from \(z_i\) because it remains fixed and known in the linear predictor. sommer passes this working response and precision to its existing weighted Gaussian mixed-model machinery. That inner fit updates fixed effects, random-effect predictions, and covariance parameters. The linear predictor is refreshed and the process repeats until the PQL convergence criterion is met.

This procedure is computationally attractive because it reuses established Gaussian solvers and their covariance machinery. But iteratively solving working Gaussian models is not the same as maximizing the integrated GLMM likelihood. PQL estimates and standard errors can be biased for binary, rare-event, or sparse-count data, especially with substantial random-effect variance. Convergence of the outer iteration is numerical convergence of the chosen approximation, not proof of unbiasedness.

An analogy: PQL repeatedly draws a tangent line to a curved road and solves the easier straight-road problem, then moves to a new point and draws another tangent. This can be fast and useful, but the answer is still based on those local straight-line approximations. It does not become an exact likelihood just because the iterations stop changing.

4.2 PQL usage and controls

PQL remains the default non-Gaussian route. task$pql$maxit, task$pql$tol, and the warm-start option govern the outer procedure. Other mmes_control() settings configure the Gaussian inner fit and therefore depend on the selected Gaussian engine and covariance model.

fit_pql <- mmes(
  response ~ treatment,
  random = ~ subject,
  rcov = ~ units,
  data = dat,
  family = binomial(),
  control = mmes_control(
    task = list(pql = list(maxit = 30L, tol = 1e-5))
  )
)

The actual supported families, response representation, rcov, weights, and random structures depend on the PQL path. In particular, the native ZIP family is Laplace-only. Consult the current GLMM guide rather than assuming that a family supported by stats::glm() can be combined with every mmes() mixed-model option.

In the example, maxit is the maximum number of outer re-linearizations, and tol is how small the changes must become before the procedure stops. These settings control the stopping rule; they do not remove approximation bias.

5.1 Integrating out random effects

The likelihood is the model’s way of scoring how plausible the observed data are for a proposed set of coefficients and variance settings. Because the random effects are hidden, we cannot just plug in one guessed value and pretend it is known. We must average the data score over the possible random effects, giving more weight to values that are plausible under their normal distribution. That average is the marginal likelihood:

\[ L(\beta,\theta,\psi;y) =\int \left\{\prod_i p(y_i\mid u,\beta,\psi)\right\} \varphi(u;0,G(\theta))\,du, \]

where \(\psi\) collects any estimated family parameters. Except for special cases, this integral has no closed form. Define the log joint density \(h(u)\) as the conditional observation log likelihood plus the Gaussian random-effect log density. Its conditional mode is

\[ \widehat u=\operatorname*{arg\,max}_u h(u), \qquad H=-\frac{\partial^2h(u)}{\partial u\,\partial u^\mathsf{T}} \bigg|_{u=\widehat u}. \]

The Laplace approximation to the marginal log likelihood is

\[ \ell_{\mathrm{LA}}(\beta,\theta,\psi) \approx h(\widehat u) +\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|H|, \]

with \(q=\dim(u)\). This is a local Gaussian approximation to the integral around its mode. It is exact when the log joint density is quadratic in \(u\), as in a Gaussian response with a Gaussian random effect. For non-Gaussian families, accuracy depends on how well the local curvature represents the integral; small numbers of independent clusters, sparse data, and strongly skewed conditional posteriors warrant special care.

The challenge is calculating that average over every possible value of the hidden effects. Usually there are too many values to sum exactly. Laplace’s idea is to look near the most plausible value \(\widehat u\), called the mode, and approximate the log density by a downward-opening parabola there.

Remember the Gaussian bell curve from algebra/statistics: its total area is its height at the center times a width factor. Near the mode, the curvature matrix \(H\) tells us how narrow or wide that bell is. In one dimension the width factor is \(\sqrt{2\pi/H}\); taking its logarithm gives \(\tfrac12\log(2\pi)-\tfrac12\log(H)\). For many random effects, the determinant \(|H|\) is the corresponding joint width factor. That is where the determinant term in the Laplace formula comes from. It corrects the height at the best-fitting random effects to account for the range of nearby plausible values too.

When the log density really is a quadratic, this bell-shaped approximation is exact. With a non-Gaussian outcome it may only be approximately bell-shaped, so the approximation can be less accurate, especially with sparse data or few independent groups.

5.2 Numerical optimization in sommer

The native route obtains fixed/random designs and covariance descriptors through sommer’s model-preparation machinery. It then uses a sparse Newton-type solve to find the conditional random-effect mode and evaluate its Hessian determinant. The outer nlminb optimization jointly updates fixed effects, free covariance parameters, and any free family parameters using native first scores of the marginal Laplace objective. It is an ML method, not REML.

This architecture preserves the meaning of the random-effect covariance structures and avoids treating a Gaussian AI-REML score as if it were the gradient of a non-Gaussian likelihood. The Gaussian objective and AI-REML engines remain distinct. The optimizer’s first derivatives include dependence of both the mode and the determinant correction on the parameters. Numerical second derivatives are calculated on demand for uncertainty, after fitting.

There are two nested numerical jobs. The inner job finds the best-fitting random-effect values for the current model parameters. The outer job changes the fixed effects, covariance values, and any estimated family parameters to improve the approximate marginal likelihood. This is maximum likelihood (ML): it chooses parameters that make the observed data as probable as possible under the model. It is not REML, which adjusts likelihood calculations to account for estimating fixed effects.

fit_laplace <- mmes(
  outcome ~ treatment,
  random = ~ subject,
  data = dat,
  family = binomial(),
  control = mmes_control(
    task = list(objective = "ml",
      glmm = list(method = "laplace",
        maxit = 1000L, modeMaxit = 100L, modeTol = 1e-8))
  )
)
logLik(fit_laplace)
AIC(fit_laplace)
fit_laplace$laplaceDiagnostics

maxit controls the outer optimizer. modeMaxit and modeTol control the inner conditional-mode solve. These are GLMM-specific controls; Gaussian AI-REML iteration and proposal settings do not control the Laplace optimizer. Inspect the optimizer result, convergence, gradient, laplaceDiagnostics, and executionPlan, rather than treating a returned object as sufficient evidence of a well-behaved fit.

In other words, both jobs need to behave: the inner solver must find its random-effect mode, and the outer optimizer must settle on model parameters. The gradient is a vector of slopes of the objective; values close to zero suggest that small parameter changes no longer improve it. These checks are useful evidence, not a guarantee that the scientific model is correct.

5.3 Supported set and deliberate restrictions

Native Laplace currently supports the documented binomial, Poisson, ZIP, NB2, Gaussian, and Gamma family/link combinations, including long-format familym() models. It reuses sommer covariance and relationship precision descriptors where supported by the native setup. It requires ML, Henderson formulation, a full-rank fixed-effects design, at least one random term, and output$inverse=0 during fitting.

Structured residual covariance, W, weights formulas, vcc, direct formulation, known-covariance solve mode, arbitrary user-defined likelihoods, hurdle models, zero-inflated NB2, and random effects in the ZIP zero component are not supported on this route. Unsupported settings are rejected rather than silently interpreted as Gaussian-engine controls. These limits matter when comparing models: a familiar covariance syntax doesn’t establish that it is accepted or scientifically identical for the selected likelihood.

These restrictions are boundaries of this particular fitting engine, not claims that the listed models are impossible in statistics. For example, the native Laplace route needs an identifiable fixed-effects design: “full rank” means no column in that design is an exact copy or combination of the others. Otherwise, different coefficient values could describe the same predictions, so the data could not determine a unique answer.

Here is the short version of the two fitting routes:

Neither route turns an approximation into certainty. The right choice depends on the scientific question, the available model support, and how much approximation error is acceptable for the data at hand.

6.1 Point estimates are not uncertainty estimates

An optimizer can converge while the likelihood is nearly flat, parameters are at a boundary, or the Hessian is indefinite. In the native Laplace method, summary() and vcov() request a numerical Hessian of the joint outer objective after fitting. Local standard errors are reported only when the fit converged, is not flagged near a boundary, and the joint curvature is positive-definite and numerically nonsingular. Covariance and family parameter standard errors are transformed to the natural scale using the delta method.

These are local curvature approximations, not finite-sample calibrated confidence intervals. A Wald interval for a variance near zero or a ZIP probability near zero is especially suspect. Boundary parameters have nonstandard null distributions, so ordinary chi-square or normal reference distributions should not be assumed without justification.

Think of the fitted likelihood as a hill whose peak gives the estimate. The estimate tells us where the peak is; uncertainty tells us how broad or flat the hill is around it. A flat hill means many nearby values fit almost equally well, so the estimate is less precise. The Hessian summarizes this local curvature. If it has the wrong shape or is nearly singular, a usual standard error cannot be trusted.

The delta method is a calculus shortcut for translating uncertainty through a smooth transformation. If a computer estimates uncertainty for a log-variance but we want the standard error of the variance itself, the method uses the slope of the exponential transformation near the estimate. This works as a local approximation; it does not magically make small-sample or boundary inference exact.

6.2 Post-fit PEV and computational separation

The Laplace fit avoids computing PEV or a full inverse during optimization. This keeps potentially expensive uncertainty calculations out of the likelihood iterations. Once fitting is finished, postPEV() can calculate the joint fixed/random coefficient curvature approximation:

fit_with_pev <- postPEV(fit_laplace, mode = 1L)
head(fit_with_pev$uPevList)

# Request the dense coefficient covariance only when its size is appropriate.
fit_with_full_covariance <- postPEV(fit_laplace, mode = 2L)
dim(fit_with_full_covariance$Ci)

Mode 1 stores random-effect prediction-error variances without retaining a dense full coefficient covariance. Mode 2 also materializes the full joint coefficient covariance and is therefore subject to a dense-memory guard. Mode 0 removes stored PEV/covariance output. These post-fit coefficient curvature calculations do not rerun or alter the fitted parameter estimates. They are distinct from the numerical outer-parameter Hessian used for fixed-effect, covariance-parameter, and family-parameter standard errors.

As with any approximate latent-variable model, PEV is conditional on the fitted model and its local curvature. It should not be presented as including all sources of uncertainty or as a substitute for external validation.

PEV means prediction-error variance: how uncertain the model is about a predicted random effect, given this model and these data. It is not the same as the population variance of the random effects. The population variance describes real differences among groups; PEV describes how uncertain we are about a particular group’s estimated adjustment.

6.3 Prediction targets

For native Laplace fits, the prediction interface distinguishes:

• Conditional prediction (re.form=NULL): includes fitted conditional modes for known random-effect levels.
• Zero-random-effect prediction (re.form=NA): sets the random effects to zero. On a nonlinear link, this is generally not a population mean.
• Restricted population prediction (re.form="population", type="response"): integrates over supported random-intercept distributions for response-scale means.

These are three different questions. A conditional prediction asks, “What does the model expect for this particular group, using its estimated adjustment?” A zero-random-effect prediction asks, “What would it predict for an average group before seeing that group’s data?” A population prediction asks, “What is the average response across the distribution of groups?” Because links are nonlinear, the second is not generally the third.

The population integration is currently restricted to independent random-intercept terms with identity relationship precision. Other covariance structures fail explicitly for this prediction target. Conditional and zero-random-effect se.fit=TRUE predictions use post-fit local uncertainty; cov.fit=TRUE additionally requests a full prediction covariance subject to the dense-memory limit. Response-scale uncertainty uses a delta method. Uncertainty for general population integration is not implemented.

For ZIP fits with a nontrivial fixed-effects zero formula, link- and response-scale uncertainty propagate the full outer-parameter covariance, including zero coefficients and cross-covariances. Conditional predictions additionally propagate the fitted mode’s parameter dependence and its local conditional latent covariance. This is a post-fit local approximation, not an exact posterior or a future-response prediction interval. The legacy scalar-ZIP coefficient-based prediction approximation remains unchanged.

Other native-family prediction approximations condition on fitted family and covariance parameters except where the documented restricted population prediction propagates its available curvature information. New levels cannot be assigned fitted conditional modes; use zero-random-effect or a supported population target when that matches the question.

7.1 Numerical diagnostics

At minimum, inspect:

1. Whether the optimizer and conditional-mode iteration converged.
2. The marginal score/gradient and any rejected evaluations.
3. Boundary flags and covariance estimates relative to their scale.
4. Whether requested standard errors are available and whether their curvature is plausible.
5. Sensitivity to reasonable starting values and model simplifications.

The Laplace diagnostics distinguish outer convergence, mode convergence, curvature validity, uncertainty validity, boundaries, and rejected proposals. They answer numerical questions. They do not prove model identifiability, correct likelihood specification, calibrated uncertainty, or adequate population-level prediction.

Think of convergence as checking whether the calculator finished the requested steps, not whether the answer is true. A model can finish successfully and still be the wrong description of the data. Numerical diagnostics ask “did the algorithm behave?”; statistical checks ask “does the model describe the real process well enough for this question?”

7.2 Response diagnostics

Randomized quantile residuals transform a discrete response through its fitted conditional CDF, randomizing within the probability mass at the observed value. Under a correctly specified model they should be broadly compatible with a standard normal distribution, subject to dependence and parameter estimation. Simulation can compare observed patterns with replicates from the fitted response model:

In plain language, a residual compares what happened with what the model expected. For count or yes/no data, ordinary subtraction can be hard to interpret because the outcomes are discrete. A randomized quantile residual first asks where the observed value falls in the model’s predicted probability distribution, then converts that position to a normal-score scale. This makes different outcome types easier to compare, although grouping and parameter estimation still affect the result.

rq <- residuals(fit_laplace, type = "quantile")
replicates <- simulate(fit_laplace, nsim = 100)

Conditional simulations hold fitted latent modes fixed. Population simulations draw new random effects from the fitted covariance model before sampling responses. Neither includes parameter-estimation uncertainty. Check residual patterns against fitted values, covariates, clusters, time, exposure, and design strata; a marginal histogram alone can hide a conditional misspecification.

A conditional simulation asks, “If these were the same groups, what responses might we see again?” A population simulation asks, “If we sampled new groups as well, what responses might we see?” Those are different scientific questions, so choose the simulation to match the question.

7.3 Statistical calibration, particularly for ZIP

Simulation is especially important when parameters are weakly identified or lie on boundaries. A useful calibration design varies sample size, cluster count, mean response rate, random-effect variance, and the true structural-zero probability. It should record not just optimizer convergence, but also bias, interval coverage when intervals are justified, boundary frequency, uncertainty availability, variance-component recovery, and predictive diagnostics.

For ZIP, include the Poisson boundary \(\pi=0\) and plausible alternatives near it. The null distribution of the zero-inflation estimate is not an ordinary interior normal distribution. Also test scenarios in which the true zero probability varies with covariates while fitting the current intercept-only model. This measures the consequences of misspecification; it does not by itself establish that adding a formula is the right remedy. Compare against scientifically plausible alternatives, including a better count mean, exposure offset, NB2 overdispersion, and omitted cluster or observation structure.

Calibration means checking the method on many simulated examples where we know the truth. For example, if a method advertises 95% confidence intervals, we can generate data with a known coefficient and count how often the intervals contain it. “Coverage” is that fraction. Always report how many intervals were successfully produced as well: 95% coverage among the easiest 20% of datasets is not 95% performance across all datasets.

The implemented fixed-effects formula was subsequently evaluated on 1,800 paired datasets, with 300 replicates per scenario and criteria specified before the analysis. Three moderate-rate scenarios passed every check. A low-rate scenario did not: only 81.3% of datasets supplied joint SEs, even though zero-slope coverage among available intervals was 93.9%. Across all planned replicates, only 76.3% both supplied an interval and covered the true zero slope. This illustrates why selected-subset coverage, convergence and uncertainty availability must be reported separately. It also limits the current formula extension’s calibration envelope; it does not license zero-component random effects or additional mixture families.

For orientation, these numbers are a stress test of this implementation, not universal properties of ZIP models. The practical lesson is that “converged” and “gave a usable, accurate uncertainty estimate” are separate outcomes. That distinction matters especially when data are sparse.

Model comparisons require care. Native Laplace uses ML and can provide likelihood-based summaries, but variance-component tests at boundaries have nonstandard reference distributions. For these cases, parametric bootstrap or a dedicated calibration is preferable to an unqualified standard likelihood-ratio test.

For the currently supported nested count-fixed-effect comparison, boot_mmes(smaller, larger, nsim=B, seed=...) draws new population random effects and responses under the null fit, then refits both native Laplace models. The statistic is \(T=2(\ell_{\rm larger}-\ell_{\rm smaller})\). With all \(B\) replicates successful, the reported Monte Carlo test is \((1+\sum_b I(T_b\geq T_{\rm observed}))/(B+1)\). Failures are not counted as non-rejections: if any fail, the ordinary p-value is unavailable and the output bounds its possible value by assigning each failure either outcome. These are sensitivity bounds, not confidence limits. Boundary replicates are retained and flagged.

Single-model boot_mmes(fit, ...) also returns descriptive percentile intervals in native optimizer coordinates. Neither form computes refit Hessians or PEV, mutates the original fit, or changes the fitting path. Inference remains conditional on the fitted generating model and the Laplace approximation. This API does not yet test covariance components or zero-inflation mechanisms, nor guarantee coverage at weakly identified boundaries.

7.4 Ordered categorical likelihood and uncertainty

An ordered outcome has categories with a natural order, such as small, medium, and large. The model describes the chance of being at or below each cutpoint. For an ordered response with \(K\) categories, the cumulative-link convention is

\[ P(Y_i\leq k\mid u)=F(c_k-\eta_i),\qquad c_1<\cdots<c_{K-1}. \]

Here \(F\) is a cumulative distribution function (for example, logistic or normal), and each \(c_k\) is a dividing point between neighboring categories. The increasing cutpoints are essential: category 3 must come after category 2. The model first calculates cumulative probabilities such as \(P(Y_i\leq 1)\) and \(P(Y_i\leq 2)\). The probability of category 2 is their difference, \(P(Y_i\leq 2)-P(Y_i\leq 1)\). The first and last categories use the lower and upper tails. This is why probabilities across all categories add to one.

The native kernel evaluates log probabilities stably rather than subtracting nearly equal CDFs in a tail. Predictor derivatives through fourth order and analytic cutpoint scores feed the same native Laplace objective. The residual link scale is fixed; there is no additional ordinal dispersion parameter.

Estimated cutpoints use coordinates \(a=(c_1,\log(c_2-c_1),\ldots,\log(c_{K-1}-c_{K-2}))\). The existing fixed intercept column is removed without rebuilding contrasts, and equivalent constant-span confounding is rejected. Fixed cutpoints retain the formula intercept. Estimated cutpoints require every declared category after filtering; no category is silently merged or dropped.

Why estimate log-gaps instead of free cutpoints? Exponentials are always positive. If the model estimates \(\log(c_2-c_1)\), then \(c_2-c_1=\exp\{\log(c_2-c_1)\}>0\) automatically. This simple algebra trick prevents the cutpoints from crossing during optimization.

Let \(\theta\) contain fixed effects, free covariance coordinates and free cutpoint coordinates. Only after fitting, marginal observed information gives \(V_\theta\), and the natural cutpoint covariance is \(J_c V_\theta J_c'\) with the appropriate coordinate selection. For a conditional prediction \(g(\theta,\hat u(\theta))\), the numerical parameter Jacobian includes dependence of the fitted modes on every parameter. With conditional random-mode precision \(H\), its local covariance is

\[ J_\theta V_\theta J_\theta' + J_u H^{-1}J_u'. \]

The same decomposition supplies ordinal random-effect PEV and fixed/random cross-covariances after fitting. It avoids conditioning away cutpoint-location uncertainty or counting fixed-effect uncertainty twice. Category-probability Jacobians include all cutpoint cross-covariances; their rows sum to zero over categories, as does uncertainty in a probability sum of one. Population probabilities integrate over new random effects rather than substituting their zero mean: probit uses normal convolution and logit uses 24-node Gaussian quadrature for the existing eligible covariance structures.

The uncertainty equation has the same two-part logic as many familiar measurement problems: uncertainty in the estimated model parameters (\(V_\theta\)) contributes one part, and uncertainty in the hidden random effects contributes another. The Jacobians \(J_\theta\) and \(J_u\) are derivatives: they describe how much the prediction changes when each parameter or random effect changes a little. Multiplying by these derivatives translates uncertainty on the input side into uncertainty in the prediction. This remains a local approximation, just like the delta method described above.

Default fitted/predicted outputs are labelled category probability matrices. Expected category codes are available explicitly, but their values and response/Pearson/working residuals depend on the numeric coding. Quantile and deviance residuals use the category distribution itself. Simulations retain ordered-factor levels, whereas bootstrap refits retain numeric codes and report absent-category failures.

For example, the model’s prediction for one observation might be \((0.2, 0.5, 0.3)\): a 20% chance of “small”, 50% of “medium”, and 30% of “large”. Those probabilities are more informative than a single average category code, because the meaning of averaging codes depends on how the categories were numbered.

This is single-family cumulative-logit/probit support, not a nominal multinomial model or ordinal mixed-family extension. Independent likelihood, probability, derivative and delta-method checks establish implementation correctness only; larger recovery/coverage calibration and cross-platform qualification remain necessary before broad inferential claims.

For the detailed checks below, RMSE (root mean squared error) summarizes the typical size of an estimate’s error, while coverage checks whether intervals include the known value often enough. A Wilson interval is one way to show sampling uncertainty around an estimated coverage rate. These terms describe the test of the method, not extra model parameters.

The subsequent frozen-source targeted calibration used 3,000 paired datasets, 300 per link/design combination, with four ordered categories and random intercepts. Four regular scenarios contributed 1,200 primary fits; all converged, with 99.7–100% joint parameter SE availability. Registered fixed-effect/cutpoint and natural random-variance bias gates passed. Nevertheless 18 of 712 registered regular gates failed: 14 probability coverage/RMSE gates, three conservative parameter-coverage checks, and one external AGQ1 agreement gate. Every raw regular coverage-range gate passed, but several pointwise Wilson intervals excluded 95%. The large number of correlated pointwise gates is not a simultaneous-coverage test; the full requested envelope still cannot be claimed as qualified.

A separate scalar-Laplace calculation at each implementation’s own fitted parameters checked 120 regular sensitivity datasets. Maximum absolute likelihood discrepancies were below \(7.4\times10^{-9}\) for sommer and \(1.76\times10^{-4}\) for the AGQ1 reference, localizing the comparison failure to reference precision without revising its original gate. Higher-order quadrature remains an approximation diagnostic, not generated-parameter truth.

In sparse designs, 66/300 logit and 258/300 probit datasets lacked a declared category and were intentionally rejected. Every primary dataset with all categories observed converged, but additional sparse uncertainty failures remained. Near-zero-variance scenarios supplied joint SEs in only 52.3%/54.3% of all datasets despite complete optimizer convergence. Coverage and uncertainty availability must therefore be reported jointly, with all planned replicates retained. The next statistical task is bounded post-fit probability-interval calibration (Wald versus transformed or bootstrap intervals), not assuming nominal coverage from a converged fit. Random-effect PEV coverage, fixed-cutpoint inference, arbitrary covariance and cross-platform release qualification remain outside this study.

7.5 Unordered categories: nominal multinomial models

The public multinomial() family fits unordered categories with the native Laplace likelihood. It requires ML, at least one random-effect term, and the Henderson formulation. It does not use PQL. W, weights formulas, vcc, structured residual covariance, and non-default Gaussian solver controls are not supported. All uncertainty work is deferred until after fitting.

An unordered, or nominal, outcome has names without a natural ranking. For example, “red”, “green” and “blue” are categories, not low, medium and high. Giving them numeric codes does not make the distance from red to green meaningful. That is why an ordered model is not a substitute.

ctrl <- mmes_control(task=list(objective="ml",
  glmm=list(method="laplace")))
fit <- mmes(color ~ temperature,
  random=~vsm(ism(.nominalCategory), ism(group)),
  data=dat, family=multinomial(reference="green"), control=ctrl)