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/AI-REML 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.

This approach deliberately reuses ai_mme_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:

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:

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.

Gaussian

The Gaussian identity model is the ordinary mmes() linear mixed model.

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

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

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.

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.

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.

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)))