--- title: "Theory and practice of generalized linear mixed models in sommer" author: "sommer development team" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Theory and practice of generalized linear mixed models in sommer} %\VignetteEngine{knitr::knitr} %\VignetteEncoding{UTF-8} --- ```{r setup, message=FALSE} library(sommer) ``` # A field guide to GLMMs in sommer A generalized linear mixed model (GLMM) is a way to study outcomes such as yes/no results, counts, or proportions when observations are grouped or related. Plants may share a field, patients may share a hospital, and repeated measurements may come from the same person. Such observations are not as independent as a pile of unrelated data points. Here is the big idea, before the symbols: 1. Choose a probability model that makes sense for the outcome (for example, Bernoulli for yes/no or Poisson for counts). 2. Connect the outcome's average to explanatory variables using a link function. 3. Add random effects to represent shared group-to-group differences. Random effects are a bit like group-specific adjustments. Two plants in the same field share the field adjustment, so their outcomes become related. The adjustment is called "random" because the model treats group adjustments as draws from a probability distribution, rather than estimating a separate, unrelated fixed number for every group. You do not need advanced calculus to follow the main ideas. When calculus does help, we translate the derivative into the practical question it answers. A quick algebra reminder: a superscript $\mathsf{T}$ means "transpose" (turn a column into a row), and variance describes how spread out values are. The notation $E(Y)$ means the average we would expect over many repetitions. This is the theoretical companion to the shorter [GLMM user guide](sommer.glmm.Rmd). It explains what the formulas mean, why the two fitting routes work, and what their limitations imply in practice. Not every family, covariance structure, and fitting route can be combined; the guide documents the combinations currently supported by `sommer`. The central challenge is that group adjustments are hidden: we do not observe them directly, so the model must average over all plausible values. That average is usually too complicated to calculate exactly. `sommer` therefore offers penalized quasi-likelihood (PQL) and an opt-in native Laplace likelihood. Both account for random effects, but they approximate the calculation differently. They are not interchangeable versions of the same likelihood analysis. | Word | Everyday meaning | | --- | --- | | Mean | The average outcome the model expects in repeated similar situations | | Link | A conversion that lets us model a probability or positive mean using an additive formula | | Fixed effect | An overall predictor effect, such as the treatment coefficient on the model's link scale | | Random effect | A group-specific adjustment, such as one field being better than another | | Variance | How much values tend to spread out around their average | | Covariance | Whether two quantities tend to move up and down together | | Likelihood | How well a proposed set of model settings explains the data we observed | # 1. The model ## Traditional nonlinear contributions and post-fit uncertainty `nonlinear_curve()` adds a shared exponential, logistic or saturation contribution to the link predictor. Existing scalar, ordinal, nominal, zero-inflated and mixed-family likelihoods are reused. Positive curve parameters and their random deviations use log coordinates. Exact conditional curvature includes likelihood score times second curve derivatives, not just a working weighted design crossproduct. Nonlinear uncertainty is requested only after fitting with `post_nonlinear()`, `postPEV()`, `vcov()` or prediction uncertainty arguments. For all free marginal parameters $a$, let $V=(-\partial^2\ell/\partial a\partial a^\mathsf{T})^{-1}$, $S$ be inverse exact conditional random-mode curvature, and $D=\partial\hat u/\partial a$. The local empirical Bayes approximation is $S+DVD^\mathsf{T}$; it differs from $S$ and is not an exact unconditional frequentist PEV. Family nuisance parameters and moving reduced-rank loading maps are included in these derivatives. For conditional response prediction $g(a,u)$, the total parameter derivative is $g_a+g_uD$. Its covariance approximation is $(g_a+g_uD)V(g_a+g_uD)^\mathsf{T}+g_uSg_u^\mathsf{T}$. For zero-deviation and prior-population means only marginal parameter uncertainty is propagated. Population derivatives differentiate the actual joint-prior integration, rather than substituting conditional mode covariance for prior variability. Future-response normal intervals also include conditional observation variability and, for population prediction, the variance across the joint prior. Probability confidence intervals retain categorical shapes; numeric future-category intervals are rejected. Checked numerical choices use `inference` in `mmes_control()`. Invalid curvature, weak identification, nonstationarity and estimated boundaries reject inference explicitly. Ordinary chi-square references are not supplied for random-variance boundaries. See `?post_nonlinear` for assumptions, scales and complete output semantics. ## 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. ::: {.link-reminder style="border: 1px solid #999; border-radius: 4px; padding: 1em; margin: 1em 0;"} **Reminder: the logit link in a mixed model** **1. Separate the observed outcome from its probability.** Suppose $y_i=1$ means that seedling $i$ survives, and $y_i=0$ means that it does not. The observed outcome $y_i$ is either 0 or 1, but its conditional mean $p_i$ is a probability: for example, $p_i=0.20$ means a 20% chance of survival for seedlings with those explanatory variables and field effect. It does **not** mean that an individual seedling survives "20%." We explain survival using treatment $t_i$ (0 for control, 1 for treated), initial height $h_i$ in cm, and a random intercept $u_{j(i)}$ for the field containing seedling $i$. Here $j(i)$ simply tells us which field seedling $i$ belongs to; seedlings in the same field share its adjustment. **2. Use the link to say what the additive model describes.** An additive formula can produce any real number, including negative numbers or numbers above 1. It therefore cannot, without restrictions, describe a probability directly. Instead, we model the **log-odds of survival**. The logit link first converts a probability into odds, then takes their natural logarithm: $$ g(p)=\operatorname{logit}(p)=\log\!\left(\frac{p}{1-p}\right). $$ For example, a probability of 0.20 becomes odds of $0.20/0.80=0.25$ (one survival for every four failures, on average), then log-odds of $\log(0.25)\approx -1.3863$. A negative log-odds value is perfectly valid; it means the probability is below 0.50, not that the probability is negative. Probabilities strictly between 0 and 1 map to log-odds anywhere on the real line. **This is where the link enters the model:** we set the transformed conditional mean equal to the additive predictor: $$ y_i\mid u \sim \operatorname{Bernoulli}(p_i),\qquad \operatorname{logit}(p_i)=\eta_i =\beta_0+\beta_1t_i+\beta_2h_i+u_{j(i)},\qquad u_j\sim N(0,\sigma_u^2). $$ Read this as: "the log-odds of survival equal an intercept, plus treatment and height contributions, plus the field adjustment." The link is not another term added to the predictor; it defines the **scale on which all those contributions add together**. There is no offset in this example. The fixed-effect part $x_i^\mathsf{T}\beta$ is $\beta_0+\beta_1t_i+\beta_2h_i$, while $z_i^\mathsf{T}u$ selects the seedling's field effect $u_{j(i)}$. **3. Fit the model using the observed zeros and ones.** Before fitting, the coefficients and field effects are unknown. For proposed values, the model computes each seedling's predictor and uses the inverse link to obtain its survival probability. A surviving seedling contributes $p_i$ to the conditional likelihood; a seedling that dies contributes $1-p_i$. For instance, if $p_i=0.80$, survival contributes 0.80 and death contributes 0.20: the model regards survival as more likely. Fitting seeks settings that explain the observed pattern, while accounting for the Gaussian distribution of field effects. In a GLMM, the likelihood also averages over the unobserved random effects; Laplace approximates that calculation, while PQL uses a working-response approximation instead of maximizing that same likelihood. These are fitting details, but both routes use the chosen link to connect means to the predictor. Crucially, we **do not take the logit of the observed zeros and ones**: $\operatorname{logit}(0)$ and $\operatorname{logit}(1)$ are not finite. The link applies to the unknown conditional mean $p_i$, not directly to the observed outcome $y_i$. Choosing the link specifies this relationship; it need not appear as a separate forward-transformation step during prediction. **4. After fitting, predict by using the inverse link.** Suppose the fitted coefficients are $\beta_0=-2$, $\beta_1=1$, and $\beta_2=0.1$, and the estimated random effect for one field is $u_j=0.5$. For a **treated seedling of height 5 cm in that field**, insert $t_i=1$, $h_i=5$, and $u_j=0.5$: $$ \eta_i = -2 + 1\times 1 + 0.1\times 5 + 0.5 = 0. $$ The **inverse link** converts this sum into a survival probability: $$ p_i = \frac{1}{1+\exp(-\eta_i)} = \frac{1}{1+\exp(0)} = 0.50. $$ Going the other way, the **link** gives $\operatorname{logit}(0.50)=\log\{0.50/(1-0.50)\}=\log(1)=0$, recovering the same linear predictor. Thus the two directions are **probability $\longrightarrow$ log-odds** (link) and **log-odds $\longrightarrow$ probability** (inverse link). The equation defining the model uses the first direction; calculating a probability from fitted coefficients uses the second. **5. Change an explanatory variable or the field effect.** For a **control seedling of the same height in the same field**, only $t_i$ changes to 0: $\eta_i=-2+0+0.5+0.5=-1$, so $p_i=1/\{1+\exp(1)\}\approx 0.269$. Thus the treatment coefficient of 1 adds 1 to the **log-odds**, not to the probability: predicted survival rises from about 27% to 50%. Likewise, one additional cm of height adds 0.1 to the log-odds, holding treatment and field fixed. For a treated 6 cm seedling in this field, $\eta_i=-2+1+0.1\times6+0.5=0.1$, giving $p_i\approx 0.525$. The height coefficient is not a fixed 10-percentage-point increase in survival probability. The random effect matters too. For the treated 5 cm seedling, setting $u_j=0$ instead gives $\eta_i=-0.5$ and $p_i\approx 0.378$, rather than 0.50. These are **conditional predictions** at specified field effects; setting $u_j=0$ is not the same as averaging probabilities over all fields. ::: ## 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. Covariance structures and random effects ## 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 also accepts structured Gaussian observation-residual blocks and compatible supplied `W` and block metadata. Non-Gaussian observation errors remain conditionally independent (apart from each nominal response's category block). For a structured Gaussian subset, its scalar factors are replaced by the joint density $$ \ell_G=-\frac12\left[n_G\log(2\pi)+\log|R| +(y_G-\eta_G)^\top R^{-1}(y_G-\eta_G)\right]. $$ The predictor score is $R^{-1}(y_G-\eta_G)$ and its negative Hessian is the full $R^{-1}$, not just its diagonal. The remaining observation factors and the random-effect prior are unchanged. ## 2.3 Formulae and covariance specification The syntax for fixed and random effects follows the usual `mmes()` conventions. For example: ```{r random-intercept-example, eval=FALSE} fit <- mmes( count ~ treatment + offset(log(exposure)), random = ~ site, data = dat, family = poisson(), method = mmes_method(objective = "ml", integration = "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. Response families and parameterizations The table gives the main conditional mean/variance parameterizations relevant to the native Laplace route. PQL delegates family calculations to family objects, and the set of supported combinations differs. | Family | Conditional mean | Conditional variance | Native Laplace link(s) | | --- | --- | --- | --- | | Binomial | $p$ | $p(1-p)/m$ for $m$ trials | logit, probit, complementary log-log | | Poisson | $\lambda$ | $\lambda$ | log | | Negative binomial NB2 | $\mu$ | $\mu+\mu^2/\theta$ | log | | Gaussian | $\mu$ | $\phi$ | identity | | Gamma | $\mu$ | $\phi\mu^2$ | log | | Zero-inflated Poisson | $(1-\pi)\lambda$ | $(1-\pi)\lambda+(1-\pi)\pi\lambda^2$ | log for $\lambda$ | The mean column says what the model expects on average; the variance column says how much individual outcomes naturally vary around that average. In the binomial row, $m$ is the number of trials. In the negative-binomial row, $\theta$ controls extra spread beyond a Poisson model. In the Gamma row, $\phi$ controls spread relative to the squared mean. For grouped binomial data, $m$ is the number of trials. `binm(successes, trials=...)` represents a proportion with its trial count attached, and the native Laplace method also accepts a two-column successes/failures response. With a mixed-family model, trial counts can be provided through the documented `familym()` interface. Supplying proportions without their denominator loses information about their precision. Why does the denominator matter? A proportion of 1/2 could mean 1 success in 2 tries or 500 successes in 1,000 tries. Both proportions are the same, but the second gives much stronger evidence about the underlying probability. The trial count tells the model how much information each proportion carries. ### Where these mean and variance formulas come from The response family describes how an outcome varies around its conditional mean. The link connects that mean, or another family parameter, to the predictor $\eta$. For example, with a logit link, $p=\operatorname{logit}^{-1}(\eta)=1/(1+\exp(-\eta))$. At $\eta=\log(3)$, $p=0.75$. The family and link answer different questions: the family describes the outcome's randomness, while the link translates predictor values into the distribution's parameters. For a binomial response, imagine $m$ independent trials, each succeeding with probability $p$. Each trial has mean $p$ and variance $p(1-p)$. Adding the trials gives a success count with mean $mp$ and variance $mp(1-p)$. Dividing that count by $m$ gives a proportion with mean $p$ and variance $p(1-p)/m$, as shown in the table. For example, if $p=0.3$ and $m=10$, the expected count is 3, its variance is 2.1, and the sample proportion has variance $0.021$. A Poisson count can be understood as many small, independent chances for an event to occur. In the limit where each chance is tiny but their total expected count stays at $\lambda$, both the count's mean and variance are $\lambda$. For $\lambda=2.5$, the model expects 2.5 events and has variance 2.5; it assigns probability $\exp(-2.5)\approx0.082$ to observing no events. The NB2 distribution allows the Poisson rate itself to vary between observations. Let the rate $\Lambda$ have a Gamma distribution with mean $\mu$ and variance $\mu^2/\theta$, then draw $Y\mid\Lambda$ from a Poisson distribution with rate $\Lambda$. Averaging over rates gives $E(Y)=E(\Lambda)=\mu$. The law of total variance gives $\operatorname{Var}(Y)=E\{\operatorname{Var}(Y\mid\Lambda)\} +\operatorname{Var}\{E(Y\mid\Lambda)\}=\mu+\mu^2/\theta$. For example, $\mu=4$ and $\theta=2$ give variance $4+16/2=12$, three times the Poisson variance of 4. Smaller $\theta$ means more variation in the hidden rates and therefore more overdispersion. The Gaussian row starts with an additive-error model, $Y=\mu+\epsilon$, where $\epsilon$ has mean zero and variance $\phi$. Therefore $Y$ has mean $\mu$ and variance $\phi$. If $\mu=10$ and $\phi=4$, outcomes have standard deviation 2 around the mean. The Gamma row is useful for positive, right-skewed outcomes. A Gamma distribution with shape $k$ and scale $s$ has mean $ks$ and variance $ks^2$. Choosing $k=1/\phi$ and $s=\phi\mu$ makes its mean $\mu$ and variance $\phi\mu^2$. For example, $\mu=10$ and $\phi=0.25$ give shape 4, scale 2.5, and variance 25. Its standard deviation is 5, so its variability grows with the size of its mean rather than staying constant. ## 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 `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. $$ For a numerical example, take $\pi=0.25$ and $\lambda_i=2$. The structural zero process contributes 25% zeros; the remaining 75% of outcomes come from a Poisson count. The observed mean is $0.75(2)=1.5$, and the variance is $0.75(2)+0.75(0.25)(2^2)=2.25$. The probability of a zero is $0.25+0.75\exp(-2)\approx0.3515$: it includes both structural zeros and ordinary Poisson zeros. The probability of exactly two events is $0.75\exp(-2)2^2/2!\approx0.2030$, because a positive count can only come from the Poisson part. 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 `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. ```{r mixed-family-example, eval=FALSE} 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() ), method = mmes_method(objective = "ml", integration = "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. Fitting route one: PQL and IRLS ## 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. ### One numerical PQL/IRLS iteration Consider four binary observations: two successes in group A and two failures in group B. Use a logit link, start with fixed intercept $\beta=0$ and group effects $u_A=u_B=0$, and hold the random-intercept variance at 1 for this illustration. The starting predictor is zero for every observation, so each fitted probability is $p=0.5$. For a Bernoulli response, the derivative of the mean with respect to the logit predictor is $h=p(1-p)=0.25$, and the variance is $V=p(1-p)=0.25$. Thus each working precision is $d=h^2/V=0.25$. The working response is $z=\eta+(y-p)/h$. Each success gives $z=0+(1-0.5)/0.25=2$; each failure gives $z=0+(0-0.5)/0.25=-2$. The inner weighted Gaussian mixed-model fit now uses two copies of $z=2$ for group A, two copies of $z=-2$ for group B, precision $0.25$ per row, and a unit-variance Gaussian prior for each group effect. Symmetry gives $\widehat\beta=0$. For group A, its random-effect equation is $0.25\{(u_A-2)+(u_A-2)\}+u_A=0$, so $\widehat u_A=2/3$; group B gives $\widehat u_B=-2/3$. After this one iteration, the predictors are $2/3$ in group A and $-2/3$ in group B. Their updated probabilities are approximately $0.661$ and $0.339$. PQL would recompute the working responses and precisions at these new predictors, then repeat. Holding the variance fixed keeps this example focused on one IRLS update; a real fit can also update covariance parameters inside the Gaussian step. 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. `working$maxit`, `working$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. ```{r pql-example, eval=FALSE} fit_pql <- mmes( response ~ treatment, random = ~ subject, rcov = ~ units, data = dat, family = binomial(), control = mmes_control( working = 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. Fitting route two: native Laplace likelihood ## 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. ### One numerical Laplace calculation Take two binary successes from one group, fix the population intercept at $\beta=0$, and give the group's random intercept $u$ a standard normal prior. Each success has probability $p(u)=\operatorname{logit}^{-1}(u)$. Including the normal-prior normalizing constant, the log joint density is $$ h(u)=2\log\{p(u)\}-\frac{u^2}{2}-\frac12\log(2\pi). $$ At $u=0$, $p=0.5$. The score is $h'(u)=2\{1-p(u)\}-u$, so $h'(0)=1$. The negative second derivative is $H(u)=1+2p(u)\{1-p(u)\}$, giving $H(0)=1.5$. One Newton step to maximize $h$ is $u_{\mathrm{new}}=u+H^{-1}h'(u)=0+1/1.5=0.667$. Repeating that update converges to $\widehat u=0.6748$, where $p(\widehat u)=0.6626$, $h(\widehat u)=-1.9699$, and $H(\widehat u)=1.4471$. There is one random effect, so the Laplace log likelihood is $h(\widehat u)+\tfrac12\log(2\pi)-\tfrac12\log H(\widehat u) \approx-1.2357$. The first two terms involving $\log(2\pi)$ cancel here because the prior is standard normal. Direct numerical integration gives log likelihood $-1.2263$ for this example, close to the Laplace result. The calculation shows both parts of the approximation: evaluate the joint density at its best random-effect value, then adjust for the local width using the curvature $H$. ## 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 interchangeable Newton linear algebra 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. ### Equivalent random-effect and observation-space algebra Write the prior precision as $P=G^{-1}$, the random design as $Z$, and the negative observation log-likelihood Hessian as $C$. Here $C$ may have nominal or Gaussian blocks, zero eigenvalues, or negative curvature (for example ZIP). The mode Newton system is $$H=P+Z^\top C Z.$$ The random backend factors $H$ using sparse Cholesky. The observation backend applies $P^{-1}$ in bounded panels to construct $K=ZP^{-1}Z^\top$, then obtains a square root $K=LL^\top$ by symmetric eigendecomposition. Only roundoff-sized negative eigenvalues of this positive-semidefinite covariance are clipped; observation curvature itself is not clipped. It factors $$B=I+L^\top C L.$$ Positive definiteness of $B$ is equivalent to that of $H$, including when $K$ is singular or $C$ is signed. There is no inversion of $C$ or $K$ and no ridge used to manufacture curvature. The determinant identity is $$\log|H|-\log|P|=\log|B|.$$ For a right-hand side $r$, let $v=P^{-1}r$ and $w=CZv$. The equivalent solve is $$H^{-1}r=v-P^{-1}Z^\top\{w-CLB^{-1}L^\top w\}.$$ These are algebraic reformulations of the **same** Laplace objective, not another approximation to the integration. Floating-point conditioning can nevertheless differ. Existing mode line searches, analytic score contractions, constraints and outer optimization use the shared solve interface. Standalone nominal and mixed nominal analytic scores also use it. `likelihood$backend` is `"auto"`, `"random"` or `"observation"`. Automatic selection uses predictor/random dimensions, a structural coupling estimate and memory admission, not elapsed-time measurements or estimator fallback. Observation covariance/root construction is reused until prior precision changes. Its dense cost is quadratic in predictor rows for storage and cubic for eigendecomposition/factorization; prior-precision solves and derivative panels add work. It avoids assembling/factoring $H$, but does not eliminate the existing sparse prior precision or guarantee bounded sparse factor fill. Post-fit observation-space joint solves eliminate the random block with $H^{-1}$ and factor the fixed-effect Schur complement $X^\top C X-X^\top C ZH^{-1}Z^\top C X$. This preserves the existing conditional coefficient uncertainty target. Ordinal/full-parameter nominal and mixed uncertainty retain their existing outer-parameter propagation. No SE/PEV or outer-parameter Hessian is computed during fitting. ### Exact score contractions and safeguarded reuse For a nominal observation with non-reference category probabilities $p$, let $T_{abc}$ be the third categorical cumulant and $V_i=Z_iH^{-1}Z_i^\top$ its conditional predictor covariance. The determinant's predictor correction is $$c_{ia}=\tfrac12\sum_{b,c}T_{abc}V_{i,bc}.$$ For scalar observations this reduces to $c_i=\tfrac12 C_{\eta,i}z_i^\top H^{-1}z_i$. With $a=H^{-1}Z^\top c$, the fixed score is $X^\top(s-c+CZ a)$. For a covariance coordinate with precision derivative $P_j$, its score is $$-\tfrac12 u^\top P_j u+ \tfrac12\partial_j\log|P|-\tfrac12\operatorname{tr}(H^{-1}P_j) +a^\top P_j u.$$ These are derivatives of the full marginal Laplace objective, not frozen-mode or Gaussian AI scores. Mixed nominal rows use their actual expanded observation/logit indices; ordinal, nuisance and Gaussian residual derivatives retain their existing analytic formulas. Reduced `vcc` coordinates apply the same parameter-map chain rule. Analytic first scores are the default; explicit numeric scores remain a diagnostic central-difference reference. Outer second derivatives remain post-fit. For $H=LL^\top$ after permutation, the LL' Takahashi recurrence gives inverse entries within the filled Cholesky pattern: $$S_{ij}=-\sum_{k>i}\frac{L_{ki}}{L_{ii}}S_{kj}\quad(j>i),\qquad S_{ii}=L_{ii}^{-2}-\sum_{k>i}\frac{L_{ki}}{L_{ii}}S_{ki}.$$ Only the selected subset is stored. Sparse predictor and precision contractions use it when every needed entry is present; otherwise they use the exact bounded-panel solve reference. Missing entries are not zeros. Extraction is operation-local, workspace-guarded and invalidated when the numeric mode factor changes. Independent numeric Hessian graph components may instead be factored separately, with summed log determinants and exact gather/scatter solves. Graph connectivity is recomputed when curvature changes; no small couplings are dropped. These are opt-in random-backend choices because workload-dependent overheads can exceed savings. Successful-mode warm reuse is separate from exact-proposal caching. Concavity of the conditional likelihood plus positive-definite prior precision ensures a unique mode; the same score certificate and line search are required at every proposal. Failed warm solves retry zero. ZIP is nonconcave: automatic reuse is disabled and explicit reuse returns the authoritative cold-verified result, avoiding a path-dependent local root. No estimator, tolerance, stochastic trace, whitening or adaptive accuracy change is introduced. 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. ```{r laplace-fit, eval=FALSE} fit_laplace <- mmes( outcome ~ treatment, random = ~ subject, data = dat, family = binomial(), method = mmes_method(objective = "ml", integration = "laplace"), control = mmes_control( optimizer = list(maxit = 1000L), mode = list(maxit = 100L, tol = 1e-8)) ) logLik(fit_laplace) AIC(fit_laplace) fit_laplace$laplaceDiagnostics ``` `optimizer$maxit` controls the outer optimizer. `mode$maxit` and `mode$tol` control the inner conditional-mode solve. 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. Random covariance equality/ratio constraints (`vcc`) are supported. Native proposals expand from independent coordinates as $w=o+T\phi$, and scores are pulled back by $T^\top$. Post-fit information is computed in the reduced space; natural covariance uncertainty uses the constraint and report-transform Jacobians. AIC/BIC count independent coordinates. The Laplace `output$setup=TRUE` parameter table contains random covariance parameters and, when requested, actual Gaussian observation-residual parameters: family dispersion, NB2 shape, thresholds and zero-inflation coefficients cannot be constrained with `vcc`. Nested-fixed-effects likelihood comparisons require identical covariance constraints. Structured Gaussian residual blocks follow $y_G\mid u\sim N(X_G\beta+Z_Gu,R)$, with the full $R^{-1}$ entering observation scores and mode curvature. The residual descriptor owns Gaussian variance, replacing (not multiplying) Gaussian family dispersion. Exact residual derivatives include both the Gaussian and mode-Hessian determinant terms. The same shared covariance assembler applies `W` on the Gaussian subset and validates independent `weights` groups. Gaussian-only integration equals exact marginal Gaussian ML. Gaussian residual workspace is dense and bounded by `execution$denseMemoryMB`; SE and PEV remain post-fit. Non-Gaussian residual dependence, 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: | Route | What it approximates | Main practical point | | --- | --- | --- | | PQL | Repeatedly replaces the curved mean-response problem with a weighted Gaussian one | Reuses Gaussian mixed-model tools; can be biased for sparse or binary data | | Laplace | The average likelihood over hidden random effects, using a local bell-shaped approximation | Provides an approximate marginal ML likelihood for its supported models | 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. Estimation, uncertainty, and PEV ## 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: ```{r post-fit-pev, eval=FALSE} 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. * **Population prediction** (`re.form="population", type="response"`): integrates over supported fitted Gaussian latent 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. For any supported positive-definite random covariance descriptor, predictor variance is obtained from the actual loadings as $v_i=z_i^\top Gz_i$, with $G$ incorporating the inverse of supplied relationship precision. For nominal observations, the corresponding small joint logit covariance $Z_iGZ_i^\top$ is required instead. Correlation between observations does not enter an individual marginal mean directly, but affects estimation, joint simulation and conditional prediction. Bounded solve panels avoid a dense covariance for all observations. Sparse factor fill is not a strict total-memory guarantee. 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. Population mean uncertainty differentiates these integrated means with respect to all estimated outer parameters and propagates their post-fit covariance. For log-linked counts/Gamma, the mean is $\exp(x_i^\top\beta+v_i/2)$ (multiplied by the nonzero-component probability for ZIP). Probit ordinal probabilities use exact normal convolution; scalar logit integration uses the existing 24-node rule. Nominal integration supports up to three non-reference logits and requires agreement of 24- and 32-node tensor rules within $10^{-6}$ absolute probability. Rule agreement is not a rigorous error bound. Temporal/spatial coordinates must be represented in the fitted covariance basis. Missing responses at represented coordinates can borrow information through the conditional process. Unknown covariance/relationship coordinates are not replaced by an arbitrary old location; new independent groups are allowed only under identity relationship precision. Population means integrate the prior fitted process law, not its posterior conditional on observed responses. See the practical GLMM guide for temporal/spatial recipes, factor-level ordering and identifiability cautions. 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 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. Diagnostics and validation ## 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. ```{r diagnostics, eval=FALSE} 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<\cdots0$ 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. The checks in this section concern single-family cumulative-logit/probit support, not the nominal likelihood. Categorical `familym()` now combines ordinal, nominal and scalar observation likelihoods under native Laplace ML. Its nominal blocks retain full cross-logit curvature, and existing random effects are shared unless distinct `.familyCoordinate` effects are explicitly requested. Mixed nominal outer derivatives are analytic derivatives of the native joint Laplace likelihood; central differences remain an explicit reference option, and scalar/ordinal derivatives remain analytic. Uncertainty remains post-fit. 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. Random covariance `vcc` constraints use the same reduced-coordinate implementation as scalar and ordinal fits. Standalone nominal models do not support `W`, weights formulas or structured observation-residual covariance because they have no Gaussian observations. Mixed `familym()` models can use these features on their Gaussian subset. Non-default Gaussian solver controls remain unsupported. 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. ### Choose one category as the comparison point Suppose green is the **reference category**. The model describes red versus green and blue versus green. It does not estimate a third, redundant set of coefficients for green versus itself. For the two comparisons, $$ \eta_{\mathrm{red}}=\log(p_{\mathrm{red}}/p_{\mathrm{green}}),\qquad \eta_{\mathrm{blue}}=\log(p_{\mathrm{blue}}/p_{\mathrm{green}}). $$ A positive red coefficient increases the red-to-green odds, holding the other predictors fixed. It does not mean "a higher category", because colors have no order. The **softmax** calculation converts all comparison scores into probabilities together: $$ p_{\mathrm{red}}=\frac{\exp(\eta_{\mathrm{red}})} {1+\exp(\eta_{\mathrm{red}})+\exp(\eta_{\mathrm{blue}})}. $$ Blue has the same denominator with its own numerator; green has numerator one. The three probabilities sum to one. The comparisons are therefore coupled: fitting separate yes/no models would not give this joint model. The numerical implementation uses stable log probabilities even when a category is extremely unlikely. Factor levels retain their declared order and names. The reference defaults to the last declared level, or can be specified by its label. The reference is moved to the final *internal* position without changing the recorded original level order. Plain numeric category codes must be contiguous integers starting at one. Character responses are not silently sorted into categories, and ordered factors must be deliberately converted if their order is to be ignored. Every declared category must remain represented after missing-data filtering; none is silently merged or dropped. For example: ```r ctrl <- mmes_control() fit <- mmes(color ~ temperature, random=~vsm(ism(.nominalCategory), ism(group)), data=dat, family=multinomial(reference="green"), control=ctrl, method=mmes_method(objective="ml", integration="laplace")) ``` The fitter creates `.nominalCategory` internally; it is not a column the user adds to `dat`. `fitted(fit)` returns one probability column per original category level. The category labels remain in the factor's original order, even though the reference is stored last internally. ### Give each comparison its own coefficients For a fixed formula such as `color ~ temperature`, the model creates a separate intercept and temperature slope for red versus green and for blue versus green. An intercept is kept here: unlike an ordinal model with estimated cutpoints, there are no cutpoint locations competing with it. The original transformations, factor contrasts and training levels are recorded for later prediction work. The computer stores all red-comparison rows first, then all blue-comparison rows. This is only bookkeeping: there is still one observed color and one joint likelihood contribution per original observation, not two independent responses. A row map records the original observation and comparison label. Missing data must remove or retain all comparison rows of an observation together. A formula offset is repeated across comparisons in this initial preparation; category-specific offset formulas and alternative fixed-effect sharing patterns are not part of this slice. Random effects also need a clear meaning. The stacked design exposes `.nominalCategory`, containing only non-reference categories. It can be used in the public random formula: | Random formula | Meaning | |---|---| | `~group` | One group adjustment shared by every non-reference comparison. It is not a separate adjustment per category. | | `~vsm(ism(.nominalCategory), ism(group))` | Separate category-by-group adjustments, independent across categories, with a common variance scale. | | `~vsm(usm(.nominalCategory), ism(group))` | Category-by-group adjustments whose variances and correlations can differ. | For example, a group that favors red may also favor blue relative to green. The last model can represent that association rather than forcing it to zero. Relationship precision supplied through `Gu` retains its existing meaning and level alignment. Changing reference recodes the comparisons; to represent exactly the same random-effect model, its design and covariance must be transformed consistently too. ### Reuse the covariance machinery, not a second set of rules The optimizer does not work directly with unrestricted variance values. For example, a log-variance coordinate keeps a variance positive after it is transformed back. Existing sommer covariance descriptors already define these transformations, the starting values, and which coordinates are fixed. The nominal objective now reuses that same native machinery rather than inventing a separate covariance parameterization. It takes the fixed coefficients followed by the free covariance coordinates in descriptor order. Fixed coordinates remain unchanged. At each proposal it builds the random-effect precision (the inverse-covariance representation), finds the conditional random-effect modes, and calculates the Laplace likelihood and its analytic score. The score is a slope telling the optimizer which parameter direction improves the likelihood; it is **not** a standard error. It includes the effect of moving the modes and changing the log-determinant adjustment. Small-block contractions and bounded sparse-solve panels supply these derivatives without building a full inverse or PEV. Invalid proposals and unsuccessful modes remain explicit. Controls stay in the existing `mmes_options`/`mmes_control()` system: `likelihood` supplies mode iteration limits and tolerance, `execution$denseMemoryMB` limits admitted workspace, and `execution$blasThreads` controls evaluation threads. The fit requires Henderson Laplace ML in fit mode and `output$inverse=0`; it rejects unused Gaussian controls and unsupported family-parameter settings. Uncertainty remains strictly post-fit. The original Gaussian/LMM and scalar GLMM fitting routes are unchanged, and no package dependency is added. `fitted()` and `predict()` return category-probability matrices; `type="link"` returns the category logits. Conditional prediction uses fitted modes for known random-effect levels, and `re.form=NA` sets those effects to zero. Population-marginal prediction is not yet implemented. `simulate()` preserves the original factor labels and can draw new random effects from the fitted population covariance. Response and Pearson residuals are returned by category; deviance, working and quantile residuals are not yet defined. `summary()` and `vcov()` compute fixed-effect uncertainty after fitting from the marginal Laplace objective. `postPEV()` reconstructs the full cross-category curvature after fitting, and `predict(se.fit=TRUE)` uses the joint fixed/random conditional curvature and a delta-method probability transformation. Prediction standard errors condition on the fitted covariance parameters. `boot_mmes()` does not yet support nominal refits. The likelihood, AIC/BIC, and nested fixed-effect likelihood comparisons are available. Independent dense calculations and numerical derivatives test the likelihood, mode, marginal score, optimizer and these public outputs. Numerical agreement does not establish confidence-interval coverage or good behavior with rare categories; nominal statistical calibration remains a separate release gate. The existing Gaussian/LMM and scalar-GLMM fitting routes are unchanged, and no package dependency is added. # 8. A practical modeling workflow 1. **Define the estimand.** Decide whether the scientific question concerns conditional effects, subject-specific predictions, or a population average. 2. **Choose the response model.** Check the outcome support, trial counts, offsets, exposure, mean-variance relation, and whether extra zeros have a defensible structural interpretation. 3. **Specify random effects scientifically.** Use the covariance and relationship structures that represent the design and dependence, not simply the largest available structure. 4. **Choose the fitting approximation deliberately.** PQL is the default and reuses the established Gaussian solver; native Laplace provides an approximate marginal ML likelihood for its supported set and may be preferable when likelihood-based inference is central. 5. **Check identifiability and stability.** Inspect convergence, gradients, boundaries, uncertainty, starts, and model sensitivity. Treat numerical success separately from statistical validation. 6. **Keep uncertainty post-fit when using Laplace.** Fit with `output$inverse=0`, then request summaries, prediction uncertainty, or `postPEV()` as needed. Respect the memory guard for dense covariance outputs. 7. **Validate predictions.** Use response-scale simulations and design-aware residual checks. For population predictions, ensure the integration target and supported covariance assumptions match the question. 8. **Report the approximation.** State the family/link, covariance structures, estimation method, optimization controls, boundary behavior, and whether inference is conditional, marginal, approximate, or simulation-calibrated. # 9. Reporting template A transparent GLMM report should state: * the conditional response distribution, link, offsets, trial counts, and any estimated family parameters; * the fixed-effects formula and the random-effect covariance/relationship structures; * whether estimates used PQL or native Laplace ML, including relevant controls and software version; * convergence, mode and boundary diagnostics, and how uncertainty was calculated; * the prediction target (conditional, zero random effect, or population-integrated) and any integration restrictions; * validation, calibration, residual and simulation checks, plus known limitations. For zero-inflated Poisson models, additionally report the interpretation of the structural-zero state, whether its probability was fixed or estimated, the boundary behavior and uncertainty limitations, and the evidence supporting the mixture rather than a different count or dependence model. # 10. Further reading The following references provide broader treatments of the statistical ideas behind GLMMs, quasi-likelihood, and likelihood approximations: * Breslow, N. E. and Clayton, D. G. (1993). Approximate inference in generalized linear mixed models. *Journal of the American Statistical Association*, 88, 9--25. * McCulloch, C. E., Searle, S. R., and Neuhaus, J. M. (2008). *Generalized, Linear, and Mixed Models*, second edition. Wiley. * McCulloch, C. E. and Neuhaus, J. M. (2011). Prediction of random effects in linear and generalized linear models under model misspecification. *Biometrics*, 67, 270--279. * Pinheiro, J. C. and Bates, D. M. (2000). *Mixed-Effects Models in S and S-PLUS*. Springer. * Wolfinger, R. and O'Connell, M. (1993). Generalized linear mixed models: a pseudo-likelihood approach. *Journal of Statistical Computation and Simulation*, 48, 233--243. The method-specific restrictions and implementation details in this vignette describe `sommer`'s current GLMM interfaces; consult the package release documentation when reproducing an analysis.