--- title: "Estimation and inference in lmmix" output: rmarkdown::html_vignette bibliography: references.bib link-citations: true vignette: > %\VignetteIndexEntry{Estimation and inference in lmmix} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` ## Gaussian mixed-model representation For a response vector (Y), `lmmix` uses \[ Y = X\beta + Zu + \varepsilon, \qquad u \sim N(0, G), \qquad \varepsilon \sim N(0, R), \qquad u \perp \varepsilon. \] The marginal model is \[ Y \sim N(X\beta, V), \qquad V(\eta) = ZG(\eta)Z^\mathsf{T} + R(\eta), \] where (eta) is an unconstrained numerical parameter vector. This is the standard Gaussian mixed-model decomposition described by @harville1977 and @pinheiro2000. ## Covariance parameterization The optimizer works on an unconstrained scale. The package maps those values to admissible covariance matrices as follows. | Component | Parameterization | |:--|:--| | Random-effect covariance `G` | Lower-triangular Cholesky factor with log diagonal | | Independent residuals `id` | One variance on the log scale | | Compound symmetry `cs` | Log variance and logistic correlation constrained to \((-1/(q-1), 1)\) | | Autoregressive `ar1` | Log variance and `tanh` correlation in \((-1, 1)\) | | Toeplitz `toep` | Log variance and partial autocorrelations mapped to a stationary Toeplitz correlation | | Banded Toeplitz `toep(k)` | First `k` covariance bands estimated; longer-lag covariances fixed to zero | | Unstructured `un` | Lower-triangular Cholesky factor with log diagonal | Here (q) is the number of distinct repeated-measurement occasions. Cholesky parameterizations guarantee positive definiteness for `G` and `un`. The partial-autocorrelation construction guarantees a positive-definite full Toeplitz correlation matrix. For `toep(k)`, candidate banded blocks are checked by Cholesky factorization during likelihood evaluation, and non-positive- definite candidates are rejected. ## Profiled ML and REML criteria For a fixed covariance matrix (V), the generalized least-squares estimator is \[ \widehat\beta(\eta) = (X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}Y. \] Let (r = Y-X\widehat\beta). The profiled negative ML log-likelihood is \[ \ell_{\mathrm{ML}}^-(\eta) = \frac{1}{2}\left[ n\log(2\pi) + \log|V| + r^\mathsf{T}V^{-1}r \right]. \] The profiled negative REML log-likelihood is \[ \ell_{\mathrm{REML}}^-(\eta) = \frac{1}{2}\left[ (n-p)\log(2\pi) + \log|V| + \log|X^\mathsf{T}V^{-1}X| + r^\mathsf{T}V^{-1}r \right]. \] REML originates with @patterson1971. General derivations are given by @harville1977 and @lamotte2007. The implementation uses Cholesky solves for (V^{-1}X), (V^{-1}Y), log determinants and the GLS normal equations. It does not form an explicit matrix inverse in the likelihood calculation. ## Optimization and covariance-parameter uncertainty `lmm_control()` can select `nlminb()` or `optim()` directly. Its default `optimizer = "auto"` starts with `nlminb()` and uses deterministic perturbed starts and alternative optimizers only when an earlier attempt does not converge. Every attempt is retained in `fit$convergence$attempts`. After optimization, `numDeriv` computes the observed likelihood Hessian on the unconstrained (eta) scale. The inverse Hessian estimates \[ \widehat\Sigma_\eta = \operatorname{Var}(\widehat\eta). \] The package checks the Hessian eigenvalues. A non-positive-definite Hessian triggers a warning because covariance-parameter standard errors and Satterthwaite calculations may then be unavailable. Natural-scale covariance standard errors are obtained by applying the delta method to the transformation from (eta) to variances and correlations. ## Satterthwaite denominator degrees of freedom For a one-row contrast (L), define \[ q(\eta) = L\operatorname{Var}(\widehat\beta)L^\mathsf{T}. \] If \[ g = \frac{\partial q(\eta)}{\partial\eta}, \] the implemented Satterthwaite approximation is \[ \nu = \frac{2q(\widehat\eta)^2} {g^\mathsf{T}\widehat\Sigma_\eta g}. \] This follows @satterthwaite1946. The gradient of (q) is evaluated numerically on the same unconstrained covariance scale used by the optimizer. With `ddf = "residual"`, the package instead uses (n-p). ## Kenward-Roger adjustment For `ddf = "kenward-roger"`, the package applies the small-sample adjustment of @kenward1997. This method is available only with REML. Let \[ \Phi = (X^\mathsf{T}V^{-1}X)^{-1} \] and let (W) denote the covariance matrix of the estimated unconstrained covariance parameters. First and second numerical derivatives of (V) are used to construct the (P_h), (Q_{hj}), and (R_{hj}) matrices in the Kenward-Roger expansion. The adjusted fixed-effect covariance is \[ \widehat\Phi_A = \widehat\Phi + 2\widehat\Phi\left\{\sum_{h,j}W_{hj} \left(Q_{hj} - P_h\widehat\Phi P_j - \frac{1}{4}R_{hj}\right) \right\}\widehat\Phi. \] Coefficient tests use this adjusted covariance. Multi-degree-of-freedom tests also use the Kenward-Roger scale factor and denominator degrees of freedom. `vcov(fit)` returns the adjusted covariance, while `vcov(fit, adjusted = FALSE)` returns the unadjusted GLS covariance. ## Type III tests Type III hypotheses are constructed with sum contrasts, then transformed back to the fitted coefficient parameterization. For a multi-row hypothesis matrix (L), the covariance of (L\widehat\beta) is eigendecomposed into orthogonal one-dimensional components. Their Satterthwaite degrees of freedom are combined with the Fai-Cornelius moment-matching formula [@fai1996]. For one model, `anova()` returns type III tests. With two or more nested models, it performs likelihood-ratio comparisons. Models that differ in fixed effects are compared under ML; REML fits are refitted automatically by default. Covariance-model comparisons must already use ML. Their chi-squared reference distribution can be approximate when the null value lies on the boundary of the parameter space [@self1987]. Setting `test = "parametric.bootstrap"` simulates the reference distribution under each smaller fitted model, refits both models by ML, and uses the finite- simulation corrected tail proportion [@davison1997]. Covariance-parameter confidence intervals are available with `confint(fit, parm = "theta_")`. Variance intervals use a log-scale delta method and correlation intervals use a Fisher-z delta method, so their limits respect the natural parameter bounds. They remain local Wald intervals based on the observed Hessian, not profile-likelihood intervals. ## Estimated marginal means `lsmeans()` constructs a fixed-effects reference grid from the model terms. Requested factors vary over their levels, nuisance factors receive equal weights, and unrequested numeric covariates are held at their observed means. The corresponding rows of the fixed-effects model matrix are averaged to form the contrast matrix. This targets population marginal means as described by @searle1980. Pairwise contrasts are differences between automatically generated marginal mean rows. Multiplicity correction is applied to pairwise p-values with `stats::p.adjust()`. When p-values are adjusted, `conf_adjust = "auto"` produces simultaneous Bonferroni confidence intervals. Pointwise intervals remain available with `conf_adjust = "none"`. ## Empirical BLUPs and prediction For each independent random-effect term (k), covariance estimation is followed by the empirical BLUP \[ \widehat u_k = G_kZ_k^\mathsf{T}V^{-1}(Y-X\widehat\beta). \] Conditional fitted values use (X\widehat\beta + Z\widehat u); marginal fitted values use (X\widehat\beta). For new data, known grouping levels use their empirical random effects. New levels cause an error unless `allow.new.levels = TRUE`, which requests their population-level prediction with a zero random contribution. ## Computational scope Version `0.4.0` forms the connected components induced jointly by repeated groups and random-effect groups. If those components are independent, the ML or REML objective is evaluated from block-specific Cholesky factorizations and aggregated GLS quantities. The computational cost of that stage is therefore driven by the component sizes rather than necessarily by one factorization of the complete (n \times n) covariance. This is not a fully sparse implementation. The fitted object retains the full dense marginal covariance for method compatibility, finite-difference Kenward-Roger calculations remain dense, and crossed random effects may connect all observations into one component. Sparse Woodbury-based calculations, analytic covariance derivatives, and generalized responses are outside the current implementation. Independent crossed and nested random-effect terms are supported, with an unstructured covariance within each term.