--- title: "Generalized linear mixed models with sommer" author: "sommer development team" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Generalized linear mixed models with sommer} %\VignetteEngine{knitr::knitr} %\VignetteEncoding{UTF-8} --- ```{r setup, message=FALSE} library(sommer) ``` # Overview `mmes()` fits generalized linear mixed models (GLMMs) when given a non-Gaussian `stats::family()` object. The model has conditional mean $$ \operatorname{E}(y_i \mid u) = \mu_i, \qquad \eta_i = g(\mu_i) = o_i + x_i^\mathsf{T}\beta + z_i^\mathsf{T}u, $$ where $g$ is the link function, $o_i$ is an optional offset, and the random effects retain sommer's usual structured Gaussian covariance model, $u \sim N(0, G)$. Thus, `random`, `rcov`, `vsm()`, relationship precision matrices, and covariance structures are specified exactly as for a Gaussian mixed model. For the detailed statistical development—including the hierarchical model, the Laplace approximation, covariance parameterization, uncertainty, prediction targets, and validation strategy—see the companion [GLMM theory vignette](sommer.glmm.theory.Rmd). The default `family = gaussian()` with its identity link takes the ordinary linear mixed-model path. Supply another family explicitly for a GLMM: ```{r basic-call, eval=FALSE} fit <- mmes(fixed = outcome ~ treatment, random = ~ subject, rcov = ~ units, data = dat, family = binomial(), control = mmes_control(update = list(blendWeight = 0))) ``` PQL uses a single numeric response, with grouped-binomial proportions carrying their trials through `binm()`. Native Laplace additionally accepts two-column binomial successes/failures and the mixed-family trials-column interface described below. # Estimation methods PQL uses Henderson working fits by default. To fit the same PQL working problems in observation space instead, select `formulation="direct"`: ```{r direct-pql-call, eval=FALSE} fit_direct <- mmes(outcome ~ treatment, random = ~ subject, data = dat, family = poisson(), control = mmes_control(formulation = "direct", working = list(maxit = 30, tol = 1e-6))) summary(fit_direct) fitted(fit_direct, type = "link") ``` Direct PQL supports ML or REML **working-model** estimation, structured random and residual covariance, grouped-binomial trials, scalar `familym()` models, supplied `W`, `weights=~block` metadata, `vcc`, and covariance warm starts. Observation filtering and offsets follow the same contract as Henderson PQL. The sparse predictor design is retained without fitting Henderson equations. `summary()` and Wald tests use the direct working-model fixed-effect covariance. As with Gaussian direct fits, `output$inverse=2` requests random-effect PEV. Prefer the default `output$inverse=0` and compute uncertainty after fitting: ```{r, eval=FALSE} predict(fit_direct, D="subject", VarU=FALSE, sed=TRUE, pairwise=TRUE) fit_direct <- postPEV(fit_direct, mode=1) # diagonal coefficient error covariance # postPEV(fit_direct, mode=2) returns full joint coefficient error covariance ``` `predict()` supports classify terms, explicit coefficient contrasts, `Dtable` rules and level selections without a prior `postPEV()` call. Direct post-fit covariance uses the retained sparse designs, relationship precisions, residual mappings and final IRLS weights in observation space, not Henderson `C`/`Cscale`. Keep `output$coefficientMatrix=TRUE` (default) to retain these inputs. `output$panelSize` bounds diagonal covariance panels; `execution$denseMemoryMB` guards post-fit workspace (or prediction's `maxMemoryMB` override). Full coefficient covariance requires quadratic storage, but diagonal PEV does not. No native factor pointers are stored in the fit. PQL prediction covariance is a link-scale working-model approximation conditional on fitted covariance parameters; coefficient contrasts do not automatically add formula offsets. Gaussian direct predictions also support `VarU` and `postVarU()`; these remain unavailable for PQL. `fitted()` and `residuals()` remain available. Direct formulation does not turn PQL into a marginal-likelihood method. PQL remains the default, and its structured covariance and relationship-matrix support are unchanged. An opt-in native Laplace method is available for binomial (logit, probit or complementary log-log), Poisson-log, zero-inflated Poisson, NB2-log, Gaussian-identity and Gamma-log models, including long-format `familym()`: ```{r laplace-call, eval=FALSE} fit_laplace <- mmes( outcome ~ treatment, random = ~ subject, data = dat, family = binomial(), method = mmes_method(objective = "ml", integration = "laplace") ) summary(fit_laplace) fit_laplace$covParNative ``` The likelihood and conditional random-effect mode are evaluated by sommer's native C++ module. It consumes the fixed/random designs and covariance precision structures prepared by sommer, including supported `CovarianceFactor` random covariance shapes and relationship precision matrices. Covariance parameters and fixed effects are optimized against the Laplace likelihood using `nlminb`, with native first derivatives of the full marginal Laplace objective (including mode dependence and the random Hessian determinant). `optimizer$maxit` controls this outer optimization and `modeMaxit`/`modeTol` control the inner random-effect mode solve. `ai_mme_sp2()` and `ai_reml_direct_sp2()` retain their Gaussian objective and are not called as if they optimized the GLMM likelihood. This creates a shared native design/covariance seam without modifying their Gaussian score equations. Preparation now stops after constructing designs and descriptors, retaining the actual response family without fitting a Gaussian working model. Reusable native family and precision modules supply stable log densities, scores, curvature derivatives and descriptor-native precision derivatives. Relationship factorizations are computed once per fitting context; covariance states and CHOLMOD symbolic analysis are reused when unchanged. Exact first-score traces use bounded dense solve panels rather than a full `H_uu^-1` or `Z H_uu^-1 Z'`. Sparse Kronecker precision blocks are still assembled: this is not a fully matrix-free implementation or a total-memory guarantee. Set the panel budget with `execution$denseMemoryMB`; other non-default Gaussian resource controls remain unsupported on this route. The fitting-local context is not stored in the fitted object. The native likelihood supports binary and grouped-binomial responses, formula offsets, and random covariance structures accepted by the Henderson setup. It requires `method$objective = "ml"`. Structured residual covariance, `W` and `weights` metadata are supported for Gaussian observation blocks (see below), not for discrete observation errors. Direct formulation and known-covariance solve mode are not yet supported. Native Laplace fitting requires `output$inverse=0`; PEV is calculated only after fitting with `postPEV(fit, mode = 1L or 2L)`. New-data predictions support conditional effects for known levels, zero-random-effect means, and response-scale population means for supported structured Gaussian random effects. Use `predict(fit, newdata, se.fit=TRUE)` for conditional or zero-random-effect standard errors. `cov.fit=TRUE` additionally returns the full prediction covariance when it fits within `execution$denseMemoryMB`. These are local observed-information approximations conditional on fitted covariance/family parameters. Population-mean SEs propagate estimated outer parameters after fitting; nominal integration has the bounds described below. No other mixed-model package is used as a fitting backend. ## Expanded response families and mixed traits ### Choosing Laplace linear algebra The same Laplace likelihood can be evaluated from two sides: - `backend="random"` factors equations indexed by hidden random coefficients. This is often best for sparse designs and relatively few random effects. - `backend="observation"` factors a dense system indexed by predictor rows, using sparse prior-precision solves. It can help when many random coefficients induce dense coupling but there are relatively few records. For nominal responses, each non-reference logit contributes a predictor row. ```{r laplace-backend, eval=FALSE} control <- mmes_control( likelihood = list(backend = "auto"), execution = list(denseMemoryMB = 1000)) fit <- mmes(count ~ treatment, random = ~ plant, data = dat, family = poisson(), control = control, method = mmes_method(objective = "ml", integration = "laplace")) fit$executionPlan$glmm$backend ``` `"auto"` is the default. It selects observation space only when the number of random coefficients is at least twice the number of predictor rows, the random design suggests dense coefficient coupling, and the estimated workspace fits the memory budget. Otherwise it retains random-space factorization. The coupling estimate is prior-precision nonzeros plus the sum of squared row support sizes, compared with one quarter of the squared random dimension. This conservative heuristic is not a measured cost model or a promise of faster fitting. Either backend can be requested explicitly; memory or numerical failures are reported, not silently switched away. The decision is fixed when an objective context is constructed. Both routes use the same family likelihoods, covariance descriptors, mode criterion, outer optimizer, determinant correction and constraints, including signed ZIP curvature and nominal/Gaussian blocks. This is **not** direct PQL and does not improve the statistical accuracy of the Laplace approximation. Keep `formulation="henderson"` at the public Laplace preparation interface; select its algebra with `likelihood$backend`, not Gaussian `solve$backend`. The observation route forms the dense covariance among predictor rows, but never a full random covariance or random-effect Hessian. Its covariance is cached while the random covariance parameters are unchanged; changing those parameters rebuilds it. Dense covariance, eigen/root and factor workspace and solve panels are checked against `execution$denseMemoryMB`. Sparse precision storage/factor fill is not a strict total-memory guarantee. Large numbers of predictor rows can make this route inappropriate. SE and PEV remain strictly post-fit. Saved observation fits contain ordinary designs/metadata rather than native pointers or a random-effect Hessian. Inference reconstructs operation-local mode solves; joint coefficient solves use a fixed-effect Schur complement. Full requested covariance output can still be large and is memory-guarded. Prediction, serialization, bootstrap and covariance constraints preserve the chosen algebra. For reproducible timing comparisons, load the working-tree package with `pkgload::load_all()` and source `inst/benchmarks/laplace-backends.R`. It checks numerical parity before measuring warm objective-plus-gradient calls with one BLAS thread, covariance-changing calls, and the cold context construction/first evaluation. The warm timings exclude preparation and the cold covariance build, and are not whole-fit speed guarantees. A source-tree CPU/OpenBLAS run of that script gave the following median seconds per **covariance-changing objective plus analytic gradient**: | Random design | Predictor rows | Random coefficients | Random backend | Observation backend | Auto | |---|---:|---:|---:|---:|---| | Dense | 32 | 256 | 0.283 | 0.026 | observation | | Dense | 64 | 512 | 1.615 | 0.160 | observation | | Sparse | 64 | 512 | 0.013 | 0.031 | random | | Dense, relatively few coefficients | 128 | 32 | 0.019 | 0.018 | random | The two dense many-coefficient cases improved roughly 10--11-fold, whereas observation space was slower for the sparse case. Backend likelihood and gradient differences were below `2e-12` in these fixtures. These small synthetic benchmarks justify the targeted route and conservative selection, not a universal timing claim; rerun them on the intended data and hardware. ### Exact fitting accelerations The outer optimizer remains `nlminb`, supplied with first derivatives; the conditional mode still uses Newton steps and the existing line search. No outer AI algorithm, relaxed inner accuracy, whitening or approximate integration is introduced. Controls retain distinct mathematical owners: | Control | Default | Alternative and meaning | |---|---|---| | `mode$warmStart` | `"auto"` | `"off"` restores zero starts; `"on"` explicitly requests successful-mode reuse. | | `optimizer$gradient` | `"analytic"` | `"central"` uses central differences of the same marginal likelihood. | | `likelihood$trace` | `"panels"` | `"selected"` contracts an exact Cholesky inverse subset. | | `likelihood$blocks` | `"off"` | `"on"` factors exact independent mode-Hessian components. | Warm starts reuse the last **successful** mode for nearby changed outer proposals, separately from the exact-parameter likelihood/score cache. Failed warm solves retry zero; convergence tolerances and line searches do not change. All current families except ZIP have concave conditional likelihoods. ZIP can have several stationary points: automatic reuse is disabled, and explicit `"on"` always cold-verifies subsequent proposals and returns the historical zero-start result. It offers no guaranteed ZIP speedup, but does not change which local root defines the likelihood. Mixed nominal analytic scores contract each observation's full logit covariance with the categorical third cumulant, then use the shared adjoint mode correction. This includes fixed, covariance, ordinal-cutpoint and other nuisance derivatives, Gaussian residual blocks and reduced `vcc` coordinates. No finite-difference fitting gradient is needed by default. Selected traces reuse an LL' Takahashi inverse subset. A required entry outside the filled Cholesky pattern triggers an **exact bounded solve**; it is never assumed zero. Block solves use the full numeric Hessian graph, including signed and zero curvature; graph membership is rebuilt when it changes. Both choices require `backend="random"`; explicit incompatible requests or insufficient workspace error. They preserve the same likelihood, but overhead and fill can make them slower. ```{r laplace-acceleration, eval=FALSE} control <- mmes_control( optimizer = list(gradient = "analytic"), mode = list(warmStart = "auto"), likelihood = list(backend = "random", trace = "selected", blocks = "on")) fit$laplaceDiagnostics$nativeWork ``` Native work records warm attempts, cold retries/checks, mode iterations, selected-subset builds, exact contraction panels and block count. The choices persist through prediction, uncertainty, bootstrap and saved fits. SE and PEV remain strictly post-fit: inverse-subset use during gradient evaluation is not public uncertainty estimation. Source `inst/benchmarks/laplace-acceleration.R` after loading the working-tree package. It uses one BLAS thread, median timings over three batches and checks every proposal/fit against the reference. Whole-fit comparisons include preparation and outer optimization. Warm mode reuse and analytic mixed scores improved the tested workloads; selected and block paths were not uniformly better and remain opt-in. Timing results below are fixture-specific rather than guarantees for arbitrary models. The final source-tree CPU/OpenBLAS run measured: | Workload | Reference seconds | Accelerated seconds | Result | |---|---:|---:|---| | Complete Poisson fit, warm reuse | 0.520 | 0.397 | 1.31x faster; mode iterations 114 to 66 | | Complete three-trait mixed fit, analytic versus numeric scores | 18.345 | 1.768 | 10.38x faster | | Mixed objective and score, analytic versus numeric | 0.645 | 0.043 | 14.96x faster | | Independent sparse predictor system, selected versus panels (1536 rows, 256 coefficients) | 0.0559 | 0.0405 | 1.38x faster | | Chain precision, selected versus panels (1536 rows, 256 coefficients) | 0.0542 | 0.0456 | 1.19x faster | | Dense design, selected versus panels (96 rows, 48 coefficients) | 0.0257 | 0.0743 | 2.89x slower | | Independent sparse system, explicit blocks versus global factorization | 0.0559 | 0.0962 | 1.72x slower | Native-system rows measure covariance-changing objective-plus-score calls, not complete fits. Timed proposals retained likelihood differences below `1e-8` and score differences below `1e-7`. Sparse native warm runs used 273 mode iterations versus 540 cold iterations across the same 90 timed proposals. Selected contractions removed all 5040 reference contraction panels in those sparse runs, but that reduction alone did not predict wall time on dense models. Consequently warm reuse and analytic scores are defaults; selected and block paths remain explicit workload-dependent choices, with global sparse Cholesky already exploiting independence well. The benchmark fixes warm starts off when comparing trace/block/gradient techniques, so each comparison measures that individual technique. ### Equality and ratio constraints in native Laplace models Native Laplace supports `vcc` for random covariance parameters in scalar, ordinal, nominal and mixed `familym()` models. Constraints are exact in the shared working-coordinate map, rather than a penalty or a post-fit repair. Inspect the parameter table using the same Laplace control as the fit: ```{r laplace-vcc, eval=FALSE} laplace_control <- mmes_control( diagnostics = list(verbose = FALSE, dateWarning = FALSE)) setup_control <- laplace_control setup_control$output$setup <- TRUE setup <- mmes(y ~ x, random = ~ group + block, data = dat, family = poisson(), control = setup_control, method = mmes_method(objective = "ml", integration = "laplace")) setup$vcParams fit <- mmes(y ~ x, random = ~ group + block, data = dat, family = poisson(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control, vcc = data.frame(parameter = c(1, 2), group = 1, scale = c(2, 1))) ``` Here the group variance is twice the block variance. Fully qualified parameter names from `setup$vcParams` can replace indices. Each constraint group needs at least two members with compatible transforms. A fixed member fixes its whole group; conflicting fixed members are rejected. This table excludes the synthetic residual scale used in Gaussian preparation. If actual structured Gaussian residual blocks are requested, their parameters are also reported under `"residual"`, but are not accepted as `vcc` members. Family dispersion, NB2 shape, ordinal cutpoints and zero-inflation coefficients are **not** `vcc` parameters. AIC/BIC count independent coordinates, not tied members. Covariance SEs, prediction uncertainty and PEV propagate the constraint map after fitting; fitting itself still computes no Hessian or PEV. `anova()` comparisons retain the nested-fixed-effects contract and require the same covariance constraints in both fits. ### Structured Gaussian observation-residual blocks Native Laplace can combine correlated Gaussian observations with the existing scalar, ordinal and nominal likelihoods. For Gaussian rows, $y_G\mid u\sim N(X_G\beta+Z_Gu,R)$; discrete observations retain their existing conditional likelihoods. With only Gaussian responses, this integration is the exact Gaussian ML likelihood, not an approximation. ```{r gaussian-residual-blocks, eval=FALSE} fit <- mmes(y ~ time + x, random = ~ group, rcov = ~ vsm(ar1m(time), ism(record)), data = dat, family = gaussian(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # Covariance between Gaussian traits paired within a record: mixed <- mmes(y ~ trait + trait:x - 1, random = ~ group, rcov = ~ vsm(usm(trait), ism(record)), data = dat_long, family = familym("trait", height = gaussian(), weight = gaussian(), count = poisson()), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) ``` `rcov` is evaluated on the Gaussian subset only. Use explicit pairing keys, especially for incomplete or unbalanced repeated measurements. The residual descriptor owns Gaussian variances; do not additionally configure Gaussian `familyParameters`. Omitting structured `rcov`, `W` and `weights` preserves the previous independent family-dispersion model. `W` uses the existing sommer residual transformation. It may cover the original or retained data, or the original or retained Gaussian subset. A full-data matrix must be identity on non-Gaussian observations and have zero cross-family entries. `weights = ~trial` declares independent blocks of `W`; it does not generate numerical weights. Missing residual coordinates are filtered with `naMethodR`, and matrices/mappings follow the retained-row order. If matrix dimensions ambiguously match different subsets, supply `W` for all original observations rather than relying on an implicit row alignment. Analytic scores include the residual determinant, mode response and Laplace determinant dependence. Numerical-Hessian SEs, prediction uncertainty and PEV remain post-fit; simulations and bootstrap refits preserve the joint Gaussian residual law. Gaussian Pearson/quantile residuals use marginal variances; deviance residuals are Cholesky-whitened and depend on observation order. The implementation reuses the residual covariance assembler with dense Gaussian-subset matrices and guards workspace through `execution$denseMemoryMB`. It does not provide Gaussian-discrete residual correlations or an arbitrary correlated discrete likelihood. Use `binm(successes, trials=trials)` or `cbind(successes, failures)` for grouped binomial responses. Positive integer trial totals and integer counts are required. The marginal likelihood retains combinatorial constants. Fitted values and response/working residuals are proportions, and Pearson variance is `p*(1-p)/trials`; trial counts are stored in `priorWeights`. Binomial logit, probit and complementary log-log links are supported. `nbinom2(theta=2)` is sommer's native NB2 family constructor: `Var(Y|u) = mu + mu^2/theta`, with a log link. Shape is jointly estimated by default. To fix it, use `familyParameters=list(theta=2, estimate=FALSE)`. `zipoisson(pi=0.1)` adds a structural-zero mixture to a Poisson count component. It estimates a separate intercept-only logit probability by default; fix it with `familyParameters=list(zeroInflation=0.2, estimate=FALSE)`. The count linear predictor models the underlying Poisson mean, while fitted and response-scale predictions report `(1-pi)*mu`. A fixed-effects zero formula is available with `zipoisson(zero=~covariates)` as described below; zero-component random effects remain unsupported. This is not a hurdle model or zero-inflated NB2. ```{r zero-inflated-poisson-call, eval=FALSE} fit_zip <- mmes( count ~ treatment, random = ~ subject, data = dat, family = zipoisson(pi = 0.1), method = mmes_method(objective = "ml", integration = "laplace"), familyParameters = list(zeroInflation = 0.1, estimate = TRUE) ) fit_zip$familyParameters ``` For a covariate-dependent structural-zero probability, use a separate fixed-effects logit formula: ```{r zero-formula-call, eval=FALSE} fit_zip_formula <- mmes( count ~ treatment, random = ~ subject, data = dat, family = zipoisson(pi=0.1, zero=~treatment), method = mmes_method(objective = "ml", integration = "laplace") ) fit_zip_formula$zeroCoefficients summary(fit_zip_formula)$zeroCoefficients ``` The zero formula has its own design and optional offset; it does not add terms to the count predictor. An intercept starts at `qlogis(pi)` and other zero coefficients at zero. A formula without an intercept starts all its coefficients at zero. Such a formula requires `estimate=TRUE`; a fixed scalar probability still uses `zero=~1` and `familyParameters`. Formula coefficients are reported on the logit scale in `zeroCoefficients`, not as a single structural-zero probability in `familyParameters`. All fitted and response predictions, simulations, and residuals use the row-specific probabilities. Missing zero covariates follow `naMethodX`, only on their declared ZIP trait's rows in a mixed-family model. Prediction preserves training factor coding and transformed terms and rejects missing covariates or unseen factor levels. Fitted ZIP models include a low-count screening report in `fit_zip_formula$laplaceDiagnostics$zip`. It summarizes retained zeros and positive counts, zero-design rank, fitted structural-zero probabilities, and the estimated split of observed zeros between the structural and Poisson components. Low-information and near-boundary flags are heuristic diagnostics, not identification tests or inferential guarantees. They use accepted-fit values only; standard errors and PEV remain post-fit. The formula extension has a bounded calibration envelope, not a universal identifiability guarantee. In development calibration with 300 datasets per scenario, the three moderate-rate settings (rate 2, 60 groups, six records per group, zero slopes 0/0.3/0.6) passed all registered gates. At rate 0.5 with zero slope 0.6, convergence was 97.7%, joint SE availability 81.3%, and both zero coefficients failed the bias tolerance. Sparse and exact-Poisson stress settings were less reliable. Use particular caution with low counts, nearly separated zero designs and boundary estimates; do not interpret an available Wald interval as proof of regular inference. With a nontrivial zero formula, `predict(se.fit=TRUE)` uses post-fit joint outer-parameter curvature, propagating zero coefficients and their cross-covariances as well as count/covariance parameters. Conditional known-level predictions also propagate mode dependence and add the local conditional latent-effect covariance. These local delta approximations differ from the parameter-fixed coefficient approximation used by the legacy intercept-only prediction path. Neither includes future-response sampling noise. Boundary or invalid joint curvature makes uncertainty unavailable; sparse-data separation can also produce a numerical-boundary error rather than a usable fit. In `familym()`, use `zipoisson()` only for a trait that needs a zero-inflated Poisson likelihood; put its `familyParameters` settings in that trait's named nested list. Gaussian-identity and Gamma-log families likewise estimate a positive dispersion by default; configure them with `list(dispersion=0.7, estimate=FALSE)` or an estimated starting value. Gaussian dispersion is residual variance. Gamma dispersion satisfies `Var(Y|u)=dispersion*mu^2`, so Gamma shape is its reciprocal. Gaussian integration is exact, but ordinary Gaussian fitting remains on its existing engine unless Laplace is explicitly selected. ```{r native-mixed-call, eval=FALSE} fit_mixed <- mmes( outcome ~ trait - 1 + trait:covariate, random = ~ vsm(usm(trait), ism(subject)), data = long_data, family = familym( "trait", binary = binomial("probit"), counts = nbinom2(theta=2), positive = Gamma("log"), continuous = gaussian() ), method = mmes_method(objective="ml", integration="laplace"), familyParameters=list( counts=list(theta=2, estimate=TRUE), positive=list(dispersion=0.3, estimate=TRUE), continuous=list(dispersion=1, estimate=FALSE) ) ) fit_mixed$familyParameters fit_mixed$dispersion ``` Every retained row has a declared observation family. Random covariance can connect traits; observation densities are conditionally independent given those random effects. Unlike PQL, this does not require a trait-structured residual term or placing a fixed-dispersion trait first. Omit `rcov` or use `~units`; structured residual correlation is still unsupported. For grouped-binomial rows among other response types, use `familym("trait", ..., trials="trial_count")`. The response contains integer successes on binomial rows and ordinary response values on other rows. The trial column holds positive integer totals for binomial rows and one for all other rows. Missing trait/trial values join the fixed-covariate filtering policy; unknown non-missing traits error. Do not additionally declare trials in a response attribute or two-column response. Native Taylor derivative kernels provide likelihood scores, observed curvature and its first/second predictor derivatives, plus log-shape/ dispersion mixed derivatives. Nuisance scores include both conditional-mode and Laplace determinant dependence. The fitted `familyParameters` reports estimated or fixed natural-scale values, including the zero-inflation probability; `summary(fit)$familyParameters` adds approximate delta-method uncertainty after fitting when joint curvature is valid. Estimated parameters count in AIC/BIC. For single-family NB2 fits, `theta.nb` is available on the fitted object; use the summary table for its post-fit standard error. The numerical Hessian and existing boundary restrictions still apply. These families do not yet imply beta-binomial, zero-inflated NB2, hurdle, other links, or arbitrary user-defined native likelihood support. ## Ordered categorical native Laplace models Single-family `ordinal("logit")` and `ordinal("probit")` models use `P(Y <= k | u) = F(cutpoint[k] - eta)`: positive coefficients favor higher categories. This is experimental cumulative-link support, not an unordered multinomial likelihood. PQL remains unsupported; categorical `familym()` support is described below. Use an explicitly ordered factor or contiguous numeric category codes `1,...,K`; labels and declared factor levels are preserved. ```{r ordinal-interface, eval=FALSE} ctrl <- mmes_control() fit_ord <- mmes(rating ~ treatment + offset(exposure), random = ~ group, data = dat, family = ordinal(), control = ctrl, method = mmes_method(objective = "ml", integration = "laplace")) fitted(fit_ord) # n by K labelled probability matrix predict(fit_ord, newdata = newdat) predict(fit_ord, re.form = NA) # zero random effects predict(fit_ord, re.form = "population") predict(fit_ord, type = "link") # latent predictor, not category codes predict(fit_ord, type = "mean") # coding-dependent expected category summary(fit_ord)$thresholds # cutpoint uncertainty, post-fit only predict(fit_ord, se.fit = TRUE, cov.fit = TRUE) simulate(fit_ord, nsim = 2, seed = 12) # ordered factors with original levels ``` By default all cutpoints are estimated, starting from retained cumulative frequencies. Preparation removes only the existing intercept column, preserving factor contrasts, and rejects any remaining location confounding. All categories must remain observed after filtering. Set `familyParameters=list(thresholds=c(-1,1), estimate=FALSE)` to fix cutpoints instead; the formula intercept is then retained, and declared ordered-factor levels without retained observations are allowed. `type="response"` and `"probabilities"` return the same probability matrix. Response, Pearson and working residuals use the numeric category coding; deviance residuals and randomized quantile residuals use category probabilities. Population probabilities use exact normal convolution for probit and 24-node Gaussian quadrature for logit under the existing restricted random-covariance eligibility. ## Categorical traits in `familym()` Native Laplace ML can combine ordinal, nominal and scalar traits in one long-format numeric response column. For each categorical trait use codes `1,...,K`; returned probability columns carry those labels. Nominal responses contribute one joint softmax likelihood per observation, including all cross-logit curvature, **not** independent binary likelihoods. ```{r categorical-familym-interface, eval=FALSE} fit_mix <- mmes(response ~ trait + trait:x - 1, random = ~ vsm(ism(group), sigma2 = 0.2), data = longdata, family = familym("trait", rating = ordinal(), choice = multinomial(reference = "3"), count = poisson()), control = ctrl, method=mmes_method(objective="ml", integration="laplace")) fitted(fit_mix) # named list: probability matrices/vectors predict(fit_mix, newdata = newlongdata, type = "response") predict(fit_mix, re.form = NA) # same list, zero random effects predict(fit_mix, se.fit = TRUE, cov.fit = TRUE, uncertainty = "full") summary(fit_mix)$thresholds # trait-labelled cutpoint SE, post-fit only postPEV(fit_mix, mode = 1) simulate(fit_mix, nsim = 2, seed = 12) # original numeric codes and row order boot_mmes(fit_mix, nsim = 199, seed = 12) ``` Nominal fixed coefficients are category-specific. Estimated ordinal locations are represented by cutpoints rather than aliased fixed intercept columns; settings such as `familyParameters=list(rating=list(thresholds=c(-1,1), estimate=FALSE))` retain fixed cutpoints and identified formula intercepts. The random formula is not automatically expanded: `~group` shares its existing random coefficient across each nominal comparison. Request separate coordinates explicitly with `~vsm(ism(.familyCoordinate), ism(group))`, or specify their covariance using the existing descriptors. Levels are `trait:location` for scalar/ordinal traits and `trait:category` for nominal non-reference categories. `.nominalCategory` belongs to single-family nominal models and errors here rather than guessing an independent expansion. Mixed residuals are trait-named, using each family's existing definitions. Nominal working/quantile residuals remain undefined. Conditional and full post-fit prediction uncertainty return matching `fit`/`se.fit` lists; `cov.fit` is the joint covariance in trait/category/within-trait-row order. Population probabilities/means use the same structured latent covariance as fitting, including nonidentity relationship precision and non-unit loadings. Population `se.fit=TRUE` propagates all estimated outer parameters and returns `uncertainty="full"`; it describes uncertainty in the integrated mean, not a future response. Nominal integration supports up to four categories with the deterministic accuracy check below. Classify aggregation remains unsupported. Mixed nominal outer gradients are native analytic, including mode and determinant dependence. `optimizer$gradient="numeric"` retains the central-difference reference alternative. Scalar/ordinal gradients remain native analytic. Gaussian and PQL routes are unchanged. All SE, parameter covariance, prediction uncertainty and PEV remain post-fit, subject to `execution$denseMemoryMB`; fitting still requires `output$inverse=0`. No numerical Hessian or PEV is calculated during fitting. Post-fit prediction uncertainty propagates the joint cutpoint, fixed-effect and covariance-parameter covariance, mode dependence, and (for conditional prediction) local latent variance. Probability SEs match the probability matrix shape; full covariance is ordered by category, with all prediction rows for each category, and labels `category:row`. `postPEV()` also propagates joint parameter uncertainty through the modes rather than treating estimated cutpoints as known. `vcov()` continues to return only the fixed-effect covariance block. `boot_mmes()` supports the same category experiment and nested slope comparisons, retaining cutpoint-location/log-gap coordinates and explicit absent-category refit failures. Implementation checks include independent `MASS::polr()` likelihood limits and probability/uncertainty references, but do not establish coverage calibration for sparse categories, separation or weakly identified random covariance. A subsequent external targeted study ran 3,000 paired datasets and 6,900 fits, including `ordinal::clmm()` higher-order quadrature and alternate-start checks. All 1,200 regular-scenario fits converged, with 99.7--100% joint SE availability and passing parameter-recovery gates. However, **18 of 712 registered regular gates failed**, mainly probability coverage/RMSE; the full targeted inference envelope is not qualified. Regular pointwise probability coverage ranged from 90.3% to 98.0%, with some Wilson intervals excluding 95%. Sparse-category fit availability reflected intentional rejection of absent declared categories; weak-variance scenarios often lacked usable uncertainty despite complete optimizer convergence. These results prioritize post-fit probability-interval calibration, not expanding the model class or relaxing the current failure guards. Random-effect PEV coverage and fixed-cutpoint inference were not assessed by this study. A separate held-out interval study used 1,200 fresh datasets and paired Wald, logit, probit and arcsine delta intervals, without changing fitting or SE calculation. All fits converged and SE availability was 99.7--100%. Of 96 registered pointwise cells per method, 8/3/4/5 respectively failed; none passed the full envelope. Logit transformation reduced the number of failed gates but is not a general coverage correction. No new interval default or public option was selected after observing these results. Refit-bootstrap calibration remains pending; these were not bootstrap intervals or joint probability regions. ## Native Laplace validation and diagnostics Use `optimizer$maxit`, `mode$maxit` and `mode$tol` to control this route. Non-default Gaussian solver, convergence, proposal, update, output and resource controls are rejected, not treated as Laplace controls. The derivative panel memory budget, BLAS thread selection and ordinary verbose/date-warning settings remain accepted. Covariance coordinates use their descriptor transforms and positive-definiteness checks, not a common statistical bound of `[-25, 25]`. Exponential coordinates outside a machine-derived numerical range are rejected before evaluation. `fit$laplaceDiagnostics` separates `outerConverged`, `modeConverged`, `curvatureValid`, `uncertaintyValid`, `boundary`, `boundaryParameters`, and `rejectedEvaluations`. Rejected proposals retain reasons; a mode iteration limit at the starting or final parameters is an error. Covariance callback errors propagate instead of being converted into an apparently valid fit. Distinct outer proposals start from the same zero random-effect vector. Only identical parameters and mode controls reuse a converged mode, so scores and numerical Hessians do not depend on the order of earlier proposals. `gradient` records the marginal log-likelihood score at the reported parameters. `nativeWork` records objective calls, mode solves, symbolic analyses, relationship factorizations and descriptor updates/reuses through the final optimizer score evaluation. Uncertainty Hessians are not calculated during fitting. A successful optimizer return is not evidence that standard errors are valid. `vcov(fit)` and `summary(fit)` calculate uncertainty after fitting, only for a converged, interior fit with a positive-definite, numerically nonsingular joint parameter Hessian, which is evaluated numerically on demand. Otherwise sommer warns and reports `NA`, rather than clamping negative inverse-Hessian diagonal entries to zero. Covariance-parameter standard errors are transformed to the natural scale by the delta method and appear in `summary(fit)$varcomp` only when the same uncertainty gate is valid. Fixed covariance parameters have no estimated standard error. These local curvature approximations are not calibrated confidence intervals. Boundary warnings use covariance eigenvalue ratios and descriptor-coordinate screening (positive coordinates relative to their starting values); they are numerical diagnostics, not calibrated hypothesis tests. Check sensitivity to covariance starts. For a valid interior fit, `vcov(fit)` returns the fixed-effect block of the joint numerical-Hessian covariance approximation. `anova(fit)` reports term-wise Wald chi-square tests using that covariance. A two-model `anova(smaller, larger)` comparison is allowed only for converged, strictly nested fixed-effect models with the same observations, response family and trials, offsets, random formula, and covariance structure. Comparisons are refused if either fit is flagged at a boundary. This does not implement variance-component likelihood-ratio tests: their null distributions may be nonstandard and require parametric bootstrap or other calibration. `wald_mmes(fit, ssType="conditional", denDF="none", terms=...)` selects conditional asymptotic term tests; small-sample denominator-df approximations are not available for the native Laplace likelihood. For an explicitly post-fit population parametric bootstrap, use `boot_mmes(fit, nsim=999, seed=17)`. It draws new random effects and responses, refits the stored native likelihood, and reports optimizer-scale parameter draws and descriptive percentile intervals without computing refit Hessians or PEV. `boot_mmes(smaller, larger, nsim=999, seed=17)` calibrates a nested count-fixed-effect likelihood-ratio statistic under the smaller model. This is not a variance-component or zero-inflation selection test. Inspect `diagnostics`: every failed or nonconvergent refit is retained. If any fail, the ordinary bootstrap p-value is unavailable and `p.value.bounds` reports sensitivity to their unknown outcomes. Boundary replicates remain flagged, not discarded. These tools still require application-specific calibration; bootstrap does not fix an unidentified generating model. `residuals(fit, type="quantile")` computes randomized probability-integral- transform residuals for discrete responses and ordinary fitted-CDF residuals for Gaussian/Gamma responses. `simulate(fit, nsim=...)` generates conditional response replicates while holding fitted latent modes fixed. `simulate(fit, re.form="population")` draws new Gaussian random effects from the fitted covariance model using sparse relationship-precision factors before sampling each response. The estimated workspace is guarded by `maxMemoryMB`, which defaults to `execution$denseMemoryMB`, and the returned simulation matrix is included in that budget. Simulations do not include parameter-estimation uncertainty. These checks do not establish population- level predictive calibration or by themselves test residual dependence. The fit stores `obsInfo` so retained rows can be traced to the original data. Training factor levels, contrast coding and transformed terms are retained for new-data predictions. The default `re.form=NULL` includes fitted conditional random-effect modes when predicting for known levels. Conditional predictions for unseen levels fail with the variable and level named; use `re.form=NA` to set their random effects to zero. On a nonlinear link this zero-random-effect response mean is **not** the population mean integrated over random effects. `re.form="population", type="response"` integrates over the fitted Gaussian latent distribution, reusing covariance descriptors and supplied relationship precision. Structured effects and random slopes contribute through their actual prediction loadings, not an intercept-only shortcut. For conditional known-level or zero-random-effect predictions, `se.fit=TRUE` returns standard errors and `cov.fit=TRUE` additionally returns the full covariance, subject to the configured dense-workspace budget. Response-scale standard errors use the delta method. They condition on fitted covariance and family parameters; population-marginal prediction uncertainty is propagated through fixed-effect, covariance and family-parameter curvature. Nominal means use each observation's joint logit covariance, including cross-category correlations, with at most three non-reference logits. The 24- and 32-node tensor rules must agree within `1e-6` absolute probability; this is an accuracy diagnostic, not a guaranteed integration error bound. Aggregation and PEV-based prediction uncertainty flags are rejected; explicit `PEV=FALSE` and `VarU=FALSE` are accepted. Working residuals are `(y - mu) / family$mu.eta(eta)`, not the working response. ## Correlated hidden processes: time, space and shared responses Repeated counts or disease observations can remain connected because an underlying condition changes gradually over time or space. Model that condition as a structured **random effect**, not as a discrete residual covariance. No new family or solver option is needed: ```{r latent-process-recipes, eval=FALSE} laplace_control <- mmes_control( diagnostics = list(verbose = FALSE, dateWarning = FALSE)) # Each plant has its own temporally correlated hidden process. # Use a factor with explicitly ordered, equally spaced time levels. dat$week <- factor(dat$week, levels = sort(unique(dat$week))) temporal <- mmes(count ~ treatment, random = ~ vsm(ar1m(week), ism(plant)), data = dat, family = poisson(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # Independent fields, with a separable hidden spatial process in each field. spatial <- mmes(diseased ~ treatment, random = ~ vsm(ar1m(row), ar1m(column), ism(field)), data = plots, family = binomial(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # In long data, this one process is shared across response types. shared <- mmes(y ~ trait + trait:treatment - 1, random = ~ vsm(ar1m(week), ism(plant)), data = long, family = familym("trait", count = poisson(), diseased = binomial()), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # These are different scientific targets: predict(temporal, newdata = visits, type = "response") # fitted plant process predict(temporal, newdata = visits, re.form = NA, type = "response") # process set to zero predict(temporal, newdata = visits, re.form = "population", type = "response", se.fit = TRUE, cov.fit = TRUE) # average over fitted process law postPEV(temporal) # uncertainty about fitted latent coefficients, post-fit only ``` In the temporal recipe, nearby weeks are linked *within* a plant; different plants are independent unless a relationship matrix is supplied. In the spatial recipe, row and column covariance shapes are multiplied. In the shared-response recipe, one latent coefficient affects each trait's link predictor. It is not automatically a separate process for each trait. Use explicit trait or `.familyCoordinate` covariance terms when distinct, correlated processes are intended. ### Ordering, missing observations and new visits `ar1m()` uses the **order of factor levels**, not elapsed numeric distances. Use it for equally spaced time/grid positions; do not code weeks 1, 2 and 8 as three consecutive positions if the intended model uses elapsed weeks. Declare the intended ordered coordinate basis when fitting and inspect the prepared random design. Unused factor levels alone do not add coordinates: the constructors drop levels with no records. A supplied relationship matrix can explicitly add unobserved individuals to the latent basis; include those individuals in the grouping factor's declared levels as well to make them available to conditional prediction. Repeated observations identify a process separately from the family variance more reliably than one binary observation per latent coefficient. Missing responses do not turn a latent process into independent effects: its covariance still connects retained observations and represented latent coordinates. Conditional prediction can use a represented plant/time combination with no response, borrowing information through the fitted process. It cannot invent a new covariance coordinate absent from the fitted basis. New independent groups can be integrated for population means when their random terms have identity relationship precision; unrelated terms may still use nonidentity relationship matrices. New time/site coordinates or unknown relationship levels are rejected rather than silently replaced by an old location. This interface is not an automatic continuous-space interpolator or an arbitrary out-of-basis forecasting API. ### Interpretation and practical checks - Dependence belongs to the hidden **link-scale process**. Discrete responses remain conditionally independent given it. This is not a discrete `rcov`. - Population prediction averages over the *prior* fitted process law, even at observed plants; it does not condition that process on their responses. - Use `vcc` and the existing descriptor controls for random covariance parameters. The Gaussian residual-block implementation is separate. - Avoid redundant intercepts, processes whose correlation approaches one, and overdispersion terms that compete with the same latent variance. Check `convergence`, boundary diagnostics and post-fit SE availability. - Population simulation draws the whole correlated process jointly. Pointwise population means need only marginal predictor variances (or a nominal observation's logit covariance), but uncertainty across predictions also reflects their shared estimated parameters. - Latent covariance prediction uses bounded dense solve panels and sparse relationship-precision factors, under `execution$denseMemoryMB`. It does not construct a dense covariance for all observations. Sparse factor fill is not covered by a strict total-memory guarantee. # How the fit works ## PQL and IRLS Sommer uses penalized quasi-likelihood (PQL). At outer iteration $t$, it linearizes the GLMM at the current link-scale predictor $\eta^{(t)}$. With $$ \mu^{(t)} = g^{-1}(\eta^{(t)}), \qquad h_i^{(t)} = \frac{d\mu_i}{d\eta_i}, \qquad V_i^{(t)} = \operatorname{Var}(y_i \mid u), $$ the working response and diagonal IRLS precision are $$ z_i^{(t)} = \eta_i^{(t)} + \frac{y_i - \mu_i^{(t)}}{h_i^{(t)}} - o_i, \qquad D_{ii}^{(t)} = \frac{\left(h_i^{(t)}\right)^2}{V_i^{(t)}}. $$ Sommer then calls its existing weighted Gaussian Henderson or direct solver on $z^{(t)}$. This updates fixed effects, BLUPs, and covariance parameters; the new predictor is $\eta^{(t+1)} = o + X\hat\beta + Z\hat u$. The outer loop stops when the deviance change is sufficiently small. Direct PQL additionally requires the final inner covariance fit to have certified convergence. This approach deliberately reuses `ai_mme_sp2()` or `ai_reml_direct_sp2()` as the weighted Gaussian inner optimizer. It does not maximize an exact marginal GLMM likelihood. PQL is practical for structured and large random-effect models, but can be biased for binary or sparse count responses, particularly with large random-effect variances. Treat standard errors and variance components as quasi-likelihood approximations in those settings. ## Dispersion and covariance structures For `binomial()` and `poisson()`, the working residual dispersion is fixed to one. For other families, sommer estimates the Gaussian working-scale residual covariance using the supplied `rcov` structure. All existing random-effect structures remain available. For example: ```{r structured-random, eval=FALSE} fit <- mmes(y ~ environment, random = ~ vsm(usm(environment), ism(genotype)), rcov = ~ units, data = dat, family = poisson(), control = mmes_control(update = list(blendWeight = 0))) ``` ## Offsets and precision weights Use `offset()` in the fixed formula as in `glm()`. A common Poisson rate model uses the logarithm of an exposure: ```{r offset, eval=FALSE} fit <- mmes(events ~ treatment + offset(log(exposure)), random = ~ site, rcov = ~ units, data = dat, family = poisson(), control = mmes_control(update = list(blendWeight = 0))) ``` `W` is an optional symmetric positive-definite observation precision matrix. It may be sparse and non-diagonal. If $W = U^\mathsf{T}U$, sommer combines it at each PQL iteration with IRLS precision as $$ W_*^{(t)} = U^\mathsf{T}D^{(t)}U. $$ This preserves both the user-supplied correlation/precision structure and the family-dependent working precision. The sparse Cholesky factor $U$ is reused across outer iterations. Missing observations are filtered before this factorization, so `W` may be supplied either for all original rows or for the retained rows. # Families and examples The interface uses standard `stats` family objects. The following examples use a random intercept to show the shared GLMM syntax. They are illustrative and therefore not run while building this vignette. ```{r family-data, eval=FALSE} set.seed(2026) n_group <- 30 n_per_group <- 8 dat <- data.frame( group = factor(rep(seq_len(n_group), each = n_per_group)), x = rep(c(0, 1), length.out = n_group * n_per_group), exposure = runif(n_group * n_per_group, 0.5, 2) ) ``` ## Gaussian The Gaussian identity model is the ordinary `mmes()` linear mixed model. ```{r gaussian-example, eval=FALSE} dat$y_gaussian <- 2 + 0.7 * dat$x + rnorm(nrow(dat), sd = 1) fit_gaussian <- mmes( y_gaussian ~ x, random = ~ group, rcov = ~ units, data = dat, family = gaussian() ) ``` A non-identity Gaussian link can also be requested, for example `gaussian(link = "log")`, provided its domain is appropriate for the response. ## Binomial ```{r binomial-example, eval=FALSE} probability <- plogis(-0.7 + 1.1 * dat$x) dat$y_binomial <- rbinom(nrow(dat), size = 1, prob = probability) fit_binomial <- mmes(y_binomial ~ x, random = ~ group, rcov = ~ units, data = dat, family = binomial(), control = mmes_control(update = list(blendWeight = 0))) ``` ## Poisson ```{r poisson-example, eval=FALSE} rate <- dat$exposure * exp(0.2 + 0.4 * dat$x) dat$y_poisson <- rpois(nrow(dat), lambda = rate) fit_poisson <- mmes(y_poisson ~ x + offset(log(exposure)), random = ~ group, rcov = ~ units, data = dat, family = poisson(), control = mmes_control(update = list(blendWeight = 0))) ``` ## Gamma The Gamma family requires a strictly positive response. ```{r gamma-example, eval=FALSE} mean_gamma <- exp(0.3 + 0.25 * dat$x) dat$y_gamma <- rgamma(nrow(dat), shape = 4, scale = mean_gamma / 4) fit_gamma <- mmes(y_gamma ~ x, random = ~ group, rcov = ~ units, data = dat, family = Gamma(link = "log"), control = mmes_control(update = list(blendWeight = 0))) ``` ## Inverse Gaussian The inverse-Gaussian family also requires a strictly positive response. The response below is positive synthetic data for demonstrating the interface. ```{r inverse-gaussian-example, eval=FALSE} dat$y_inverse_gaussian <- exp(0.2 + 0.3 * dat$x + rnorm(nrow(dat), sd = 0.2)) fit_inverse_gaussian <- mmes(y_inverse_gaussian ~ x, random = ~ group, rcov = ~ units, data = dat, family = inverse.gaussian(link = "log"), control = mmes_control(update = list(blendWeight = 0))) ``` ## Quasi families Quasi families use the same mean/link and variance definitions as their corresponding GLM families, but do not define a likelihood. Consequently, devience is useful for PQL convergence, while likelihood-based comparisons such as AIC or likelihood-ratio tests are not appropriate. ```{r quasi-example, eval=FALSE} dat$y_quasi <- 1 + 0.5 * dat$x + rnorm(nrow(dat), sd = 0.5) fit_quasi <- mmes(y_quasi ~ x, random = ~ group, rcov = ~ units, data = dat, family = quasi(link = "identity", variance = "constant"), control = mmes_control(update = list(blendWeight = 0))) fit_quasibinomial <- mmes(y_binomial ~ x, random = ~ group, rcov = ~ units, data = dat, family = quasibinomial(), control = mmes_control(update = list(blendWeight = 0))) fit_quasipoisson <- mmes(y_poisson ~ x + offset(log(exposure)), random = ~ group, rcov = ~ units, data = dat, family = quasipoisson(), control = mmes_control(update = list(blendWeight = 0))) ``` # Extracting results and controlling PQL For a GLMM, `fitted()` returns conditional fitted means on the response scale. Use `type = "link"` for $\eta$, including any offset. Residuals can be requested on response, signed-deviance, or final working scales. ```{r extract, eval=FALSE} fitted(fit_poisson) fitted(fit_poisson, type = "link") residuals(fit_poisson, type = "deviance") residuals(fit_poisson, type = "working") fit_poisson$pqlMonitor fit_poisson$pqlConverged fit_poisson$family ``` Use `working` to set outer PQL iterations and relative deviance tolerance. Set `update$blendWeight = 0` to retain pure-AI covariance updates in PQL; positive EM-direction blending is not supported. `convergence$maxIters`, `convergence$likelihoodTol`, and `proposal$stepWeight` control the Gaussian variance-component fit inside each outer iteration. `solve$method` selects the Henderson backend only; leave it at its default for direct PQL. ```{r pql-control, eval=FALSE} fit <- mmes(y_poisson ~ x + offset(log(exposure)), random = ~ group, rcov = ~ units, data = dat, family = poisson(), control = mmes_control( update = list(blendWeight = 0), working = list(maxit = 30, tol = 1e-6), convergence = list(maxIters = 20), solve = list(method = "auto") )) ``` # Practical guidance - Begin with a simple random-intercept GLMM and inspect `pqlMonitor` before adding complex covariance structures. - Check `pqlConverged`; `FALSE` means the outer PQL criterion was not met within `working$maxit` iterations. Convergence on the last allowed iteration is still reported as converged. For direct PQL, a failed final inner fit also prevents convergence; inspect `pqlMonitor[,"innerConverged"]` and `engineDiagnostics$direct` for the working-model stop reason. - Keep binomial responses away from complete separation where possible; very small derivatives can create unstable IRLS weights. - For binary outcomes with few observations per random-effect level or large random-effect variances, compare conclusions with a method based on a Laplace approximation or Bayesian posterior simulation when feasible. - Use response-scale fitted values for prediction and link-scale fitted values when interpreting additive fixed and random effects.