# Introduction to MEMWAS

###

Mixed Effects Models With Autocorrelation Structures

#### Enoch Kang

#### Aug. 04, 2026

# Installation

Install version 0.9.3 from the local source archive. A working C++ compiler supported by the installed version of ‘R’ is required because this release uses its registered native C++ backend for fitting.

    install.packages(
      "MEMWAS_0.9.3.tar.gz",
      repos = NULL,
      type = "source"
    )
    #library(MEMWAS)

Or install the current development version from GitHub:

    remotes::install_github("EnochKang/MEMWAS")
    library(MEMWAS)

The capability registry is the most direct way to inspect the fitting contract implemented by the installed version:

    MEMWAS_capabilities()

> Version scope. This vignette documents the public interface and numerical paths implemented in MEMWAS 0.9.3. This version fits models directly through `fit_MEMWAS()`. It does not implement the former `set_MEMWAS()` settings-object workflow, a selectable R-versus-C++ engine, sandwich standard errors, a separate profile-likelihood standard-error method, or a requestable skew-corrected-Laplace approximation.

# Purpose

MEMWAS is intended for longitudinal and clustered studies in which repeated observations may share several distinct sources of dependence. Examples include students followed across semesters, patients followed across visits, and participants completing repeated behavioral assessments. Stable differences in level or trajectory can be represented by Gaussian random effects. Time-ordered persistence that remains after those effects have been included can be represented by one or more latent serial coefficient processes.

The package supports Gaussian, binomial, Poisson, negative-binomial, Gamma, and exponential conditional response models. It can combine fixed effects, offsets, observation weights, clustered, crossed, or nested Gaussian random effects, fixed-effect penalties, and multiple independent named serial processes in one fitted model. A serial process may load every observation equally or may multiply one numeric predictor or one-column design, allowing time-varying coefficients with their own serial covariance.

MEMWAS also integrates optional nonlinear screening, assumption-oriented diagnostic screens, grouped tuning, autocorrelation-structure ranking, cluster bootstrap inference, directional tests, and several prediction targets. These operations are reported explicitly so that the requested approximation, the kernel actually used, convergence status, numerical limitations, unavailable tests, and post-selection steps remain visible.

# Why random effects and serial correlation are different?

Random effects and serial processes both induce dependence, but they enter different parts of the latent model and encode different scientific patterns.

For a random intercept with variance \(\tau_0^2\),

\[ \operatorname{Cov}(b_{0i},b_{0i})=\tau_0^2 \]

is contributed to every pair of observations sharing that random intercept, regardless of their temporal separation. With a random intercept and random time slope,

\[ b_i= \begin{pmatrix} b_{0i}\\ b_{1i} \end{pmatrix}, \qquad D= \begin{pmatrix} \tau_0^2 & \tau_{01}\\ \tau_{01} & \tau_1^2 \end{pmatrix}, \]

the random-effect covariance between times \(t_{ij}\) and \(t_{ik}\) is

\[ \operatorname{Cov} \left( b_{0i}+t_{ij}b_{1i}, b_{0i}+t_{ik}b_{1i} \right) = \tau_0^2+(t_{ij}+t_{ik})\tau_{01}+t_{ij}t_{ik}\tau_1^2. \]

By contrast, an ordinary AR(1) serial process with marginal variance \(\sigma_u^2\) has

\[ \operatorname{Cov}(u_{ij},u_{ik}) = \sigma_u^2\rho^{d_{ijk}}, \]

where MEMWAS defines the scaled distance as

\[ d_{ijk}=\frac{|t_{ij}-t_{ik}|}{s_t}, \]

and \(s_t>0\) is the component’s `serial_time_scale`. With the default non-negative parameterization, \(0<\rho<1\), so correlation decays as the scaled time gap increases. If `allow_negative_autocor = TRUE`, MEMWAS uses \(-1<\rho<1\), but then the scaled lags must form an integer grid so that negative-base powers are well defined.

Scientifically, random effects describe stable heterogeneity in levels or trajectories. Serial processes describe locally persistent latent deviations or persistent predictor coefficients after the fixed and random effects have been included. For example, a random intercept can represent a patient’s usual symptom level, while an AR(1) outcome-loaded process can represent short-lived departures around that level. A predictor-loaded serial process can instead represent an exposure effect that changes persistently over time.

The two sources can compete statistically. A nearly persistent AR(1) process can resemble a random intercept; a random slope can resemble smooth serial dependence; and several serial components can become weakly separated when their loading columns are similar. Boundary estimates, large uncertainty, fixed- or random-effect rank deficiencies, covariance-Jacobian rank warnings, sensitivity to starting values, and sensitivity to plausible alternative structures should therefore be examined.

For a linear Gaussian identity-link model with ordinary linear serial components, the exact marginal covariance is a sum of random-effect, serial-process, and independent observation-error contributions. For nonlinear links or for the nonlinear `expOU` loading, this simple observed-scale covariance identity no longer applies.

# What does the MEMWAS work on?

MEMWAS 0.9.3 fits longitudinal mixed-effects models with:

- Gaussian, binomial, Poisson, negative-binomial, Gamma, or exponential conditional responses;
- supported family-specific links;
- one or more clustered, crossed, or nested random-intercept or random-slope terms;
- diagonal or unstructured covariance for each random-effect term;
- offsets and non-negative observation/frequency weights;
- one or more ordinary or predictor-loaded serial components;
- L1, L2, or elastic-net fixed-effect penalization; and
- optional regularization of serial correlation parameters.

The public serial structures are `"NONE"`, `"AR(1)"`, `"OU"`, `"expOU"`, `"AR(p)"`, `"ARMA(1,1)"`, `"CS"`, `"TOEP"`, and `"UN"`. Multiple serial components are independent a priori, but they can be integrated jointly when they are connected through common observations or random-effect terms.

## Statistical framework

### Notation

Let \(r=1,\ldots,N\) index analysis rows and let \(i(r)\) identify the primary independent cluster supplied through `id`. Let:

- \(Y_r\) be the observed response;
- \(x_r\in\mathbb{R}^{p}\) be the fixed-effect design row;
- \(\beta\in\mathbb{R}^{p}\) be the fixed-effect vector;
- \(o_r\) be the total additive offset on the link scale;
- \(h=1,\ldots,H\) index random-effect terms;
- \(z_{rh}\) be the design row for random-effect term \(h\);
- \(g_h(r)\) be the grouping level used by term \(h\);
- \(b_{h,g}\sim N(0,D_h)\) be the Gaussian vector for level \(g\) of term \(h\);
- \(c=1,\ldots,C\) index serial components;
- \(s_{rc}\) be the one-column loading for serial component \(c\); and
- \(u_{i(r)c,r}\) be the serial Gaussian driver at the row’s primary cluster and time.

Random-effect vectors are independent across grouping levels and terms under the implemented prior. Serial components are independent of the random effects and of one another under the prior. Dependence across observations arises because rows share latent variables and because each serial driver has a structured covariance over time.

The implemented latent linear predictor is

\[ \eta_r = o_r+x_r^\top\beta + \sum_{h=1}^{H}z_{rh}^\top b_{h,g_h(r)} + \sum_{c=1}^{C}s_{rc}\,T_c(u_{i(r)c,r}). \]

For all ordinary serial structures,

\[ T_c(u)=u. \]

For an `expOU` component, MEMWAS uses the centered lognormal transformation

\[ T_c(u) = \exp\!\left(u-\frac{\sigma_{c}^{2}}{2}\right)-1, \]

where \(u\) is an OU Gaussian driver with variance \(\sigma_c^2\). Consequently,

\[ E\{T_c(u)\}=0, \qquad \operatorname{Var}\{T_c(u)\}=\exp(\sigma_c^2)-1, \]

and for two driver values with covariance \(K_{jk}\),

\[ \operatorname{Cov}\{T_c(u_j),T_c(u_k)\} = \exp(K_{jk})-1. \]

Thus `expOU` is not an additive Gaussian serial effect. It is a mean-centered, positive-skewed serial coefficient transformation driven by an OU process.

### Gaussian outcomes

For the identity link with only ordinary linear latent terms,

\[ Y_r=\eta_r+\varepsilon_r, \qquad \varepsilon_r\stackrel{\mathrm{ind}}{\sim}N(0,\sigma^2). \]

Stacking observations gives

\[ Y=X\beta+o+Zb+\sum_{c=1}^{C}S_cu_c+\varepsilon, \]

where \(S_c=\operatorname{diag}(s_{1c},\ldots,s_{Nc})\), with the relevant cluster-specific blocks understood. The exact marginal covariance is

\[ V = \sum_{h=1}^{H}Z_hG_hZ_h^\top + \sum_{c=1}^{C}S_cK_cS_c^\top + \sigma^2I. \]

Here \(G_h\) is the block-diagonal covariance assembled from \(D_h\), and \(K_c\) is the block-diagonal serial covariance assembled over primary clusters. This is the model handled by MEMWAS’s exact Gaussian kernel when `approximation = "laplace"` or `"auto"` dispatches a linear Gaussian identity model to that kernel.

A Gaussian model with `link = "log"`, or a Gaussian identity model containing `expOU`, is nonlinear in the Gaussian latent variables. It is therefore handled through the requested latent approximation rather than the exact linear-Gaussian marginal formula.

Observation weights multiply conditional log-likelihood contributions. They should not automatically be interpreted as known inverse residual variances.

### Non-Gaussian outcomes

For non-Gaussian responses,

\[ Y_r\mid a\sim\mathcal{F}(\mu_r,\varphi), \qquad g(\mu_r)=\eta_r, \]

where \(a\) collects every random-effect and serial latent variable connected to the row. Conditional on \(a\), rows contribute independently to the implemented conditional likelihood, subject to the selected family and observation weights.

The marginal likelihood has the generic form

\[ L(\psi) = \int \left\{ \prod_{r=1}^{N}f(Y_r\mid a;\psi)^{w_r} \right\} \phi(a;0,C_\psi)\,da, \]

where \(w_r\ge 0\) is the observation/frequency weight and \(\psi\) collects fixed effects, covariance parameters, and family-specific dispersion parameters. MEMWAS evaluates or approximates this integral by connected latent component, rather than treating each latent term as separate when observations connect them.

For nonlinear links,

\[ E(Y_r\mid x_r) = E_a\left[g^{-1}\{\eta_r(a)\}\right] \neq g^{-1}(x_r^\top\beta+o_r) \]

in general. Fixed-effect coefficients are therefore conditional or latent-specific link-scale parameters, not population-averaged response-scale effects.

The observed-scale covariance follows the law of total covariance,

\[ \operatorname{Cov}(Y) = E\{\operatorname{Cov}(Y\mid a)\} + \operatorname{Cov}\{E(Y\mid a)\}, \]

and does not generally reduce to a sum of Gaussian latent covariance matrices.

## Family-specific conditional models

### Overview on family-specific models

Family | Implemented links | Conditional mean and variance
--- | --- | ---
Gaussian | identity, log | \(Y_r\mid a\sim N(\mu_r,\sigma^2)\), so \(\operatorname{Var}(Y_r\mid a)=\sigma^2\).
Binomial | logit, probit, cloglog, loglog, cauchit | \(Y_r\mid a\sim\operatorname{Binomial}(N_r,p_r)\), \(\mu_r=N_rp_r\). Matrix responses and proportion-plus-trial-weight formulations are supported by response preparation.
Poisson | log | \(Y_r\mid a\sim\operatorname{Poisson}(\mu_r)\), so \(\operatorname{Var}(Y_r\mid a)=\mu_r\).
Negative binomial | log | NB2 parameterization with size \(\theta>0\): \(\operatorname{Var}(Y_r\mid a)=\mu_r+\mu_r^2/\theta\).
Gamma | log | Shape \(\kappa>0\), rate \(\kappa/\mu_r\): \(\operatorname{Var}(Y_r\mid a)=\mu_r^2/\kappa\).
Exponential | log | Rate \(1/\mu_r\): \(\operatorname{Var}(Y_r\mid a)=\mu_r^2\).

For negative-binomial models, `theta = NULL` requests estimation of \(\theta\); a supplied positive value fixes it. For Gamma models, `shape = NULL` requests estimation of \(\kappa\); a supplied positive value fixes it. The exponential family is the Gamma shape-one special case, but it is implemented as its own family choice.

### Choosing families and interpreting coefficients

For an identity-link Gaussian model, \(\beta_k\) is a conditional mean difference in the original response units. For a logit-binomial model,

\[ \exp(\beta_k) \]

is a conditional odds ratio for a one-unit increase in predictor \(k\), holding all other fixed predictors and latent variables constant. For Poisson, negative-binomial, Gamma, and exponential models with a log link,

\[ \exp(\beta_k) \]

is a conditional multiplicative change in the mean.

Alternative binomial links do not have an odds-ratio interpretation. Their coefficients remain conditional effects on the selected link scale.

A time-varying predictor coefficient may mix between-cluster and within-cluster associations. When the scientific question distinguishes them, decompose

\[ x_{ij}=\bar{x}_i+(x_{ij}-\bar{x}_i) \]

and include both terms. A predictor-specific serial component can then be reserved for residual temporal variation in the coefficient after those systematic components are modeled.

## Likelihood and approximation methods

### Gaussian models

For a linear Gaussian identity-link model, MEMWAS uses a direct exact-Gaussian native kernel. It forms the marginal contribution through conditional linear solves over connected latent components rather than requiring a single dense \(N\times N\) covariance matrix. Structure-specific serial precision/factorization routines and cached time-pattern factorizations are reused where possible.

Version 0.9.3 stores the global random-effects design canonically as row pointers, zero-based column indices, and nonzero loadings. The dense \(N\times q\) random-effects matrix is materialized only for an explicit diagnostic or compatibility fallback. The persistent native context copies immutable model arrays once, precomputes row-to-latent mappings, and reuses allocated linear-algebra workspaces across objective and derivative evaluations.

For covariance parameter \(\theta_j\), define \(V_j=\partial V/\partial\theta_j\). The exact REML score used by the native optimizer is

\[ \frac{\partial\ell_R}{\partial\theta_j} =\frac{1}{2}\left \{Y^\top P V_j P Y-\operatorname{tr}(P V_j)\right \}, \qquad P=V^{-1}-V^{-1}X(X^\top V^{-1}X)^{-1}X^\top V^{-1}. \]

The ML score uses the corresponding exact profiled-likelihood derivative. Random-effect Cholesky parameters, residual scale, and all ordinary serial covariance structures are differentiated analytically. The observed Hessian differentiates this exact score with centered steps and retains the former objective-stencil method only as a parameter-boundary fallback. This changes the derivative implementation, not the likelihood, REML correction, parameter constraints, or reported inferential target.

Serial precision matrices with exact band structure are factored by an exact banded Cholesky kernel. Entries are classified as structural zeros only when they are exactly zero; no numerical truncation is used. If a component is not banded or the banded factorization fails, the existing dense route is used.

With no L1 penalty, fixed effects can be profiled through the penalized generalized least-squares system

\[ \left(X^\top V^{-1}X+\lambda_2D_\beta\right)\widehat\beta = X^\top V^{-1}(Y-o), \]

where \(D_\beta\) excludes the intercept. With an L1 penalty, the profiled fixed-effect subproblem is solved by coordinate descent and soft thresholding.

Maximum likelihood and restricted maximum likelihood have distinct roles:

- `method = "ML"` is available for all supported model classes;
- `method = "REML"` is restricted to unpenalized Gaussian identity-link models handled through the Laplace/exact-Gaussian dispatch; and
- REML objective values should not be used to compare models with different fixed-effect design matrices.

### Non-Gaussian models

Let

\[ h(a;\psi) = \sum_{r=1}^{N}w_r\log f(Y_r\mid a;\psi) + \log\phi(a;0,C_\psi). \]

The exact marginal likelihood is

\[ L(\psi)=\int\exp\{h(a;\psi)\}\,da, \]

which normally lacks a closed form. MEMWAS provides the following public approximation choices: `"auto"`, `"laplace"`, `"saddlepoint"`, `"adaptive_gaussian_quadrature"`, `"variational_inference"`, and `"pql"`.

`approximation = "auto"` dispatches a model with no latent variables to direct likelihood, a linear Gaussian identity model to the exact-Gaussian kernel, and other latent models to Laplace integration.

#### Laplace approximation

Let \(\widehat a\) maximize \(h(a;\psi)\), and define

\[ H = -\left. \frac{\partial^2h(a;\psi)}{\partial a\partial a^\top} \right|_{a=\widehat a}. \]

For a latent component of dimension \(d\), the first-order Laplace approximation is

\[ \log L(\psi) \approx h(\widehat a;\psi) + \frac{d}{2}\log(2\pi) - \frac{1}{2}\log|H|. \]

The native implementation finds conditional modes component by component, uses structure-aware prior operations, and supplies an implicit outer score and matrix-free Hessian-vector products. It includes warm-mode retries, cold-start verification, multiple structural starts, and optimizer fallback paths. For the linear Gaussian identity model, the Laplace request is dispatched to an exact one-solve Gaussian calculation; the result is not merely a first-order approximation.

#### Variational inference

MEMWAS uses a full-covariance Gaussian variational distribution within each connected latent component,

\[ q(a)=N(m,S), \]

and maximizes the evidence lower bound

\[ \mathcal{L}(q,\psi) = E_q\{\log p(Y,a\mid\psi)\} - E_q\{\log q(a)\}. \]

The native updates use expected scores and information matrices. For `expOU`, the implementation uses exact first and second moments of the transformed loading under \(q\), followed by moment-matched Gauss–Hermite integration of the resulting scalar predictor contribution.

The variational objective is an ELBO, not a marginal log-likelihood. Accordingly, `logLik`, AIC, and BIC are not reported as if they were available for variational fits.

#### Adaptive Gauss–Hermite quadrature

Adaptive Gaussian quadrature centers and scales tensor-product Gauss–Hermite nodes around the native conditional mode. With \(K\) nodes per dimension and a connected latent component of dimension \(d\), the direct tensor rule requires

\[ K^d \]

node evaluations. `quadrature_points` supplies \(K\). The implementation enforces configurable limits on the maximum component dimension and total nodes, so quadrature is primarily practical for small connected latent components.

#### Saddlepoint approximation

The public `"saddlepoint"` path uses a response-level saddlepoint likelihood and then performs native Laplace integration over Gaussian latent effects. Binomial, Poisson, and negative-binomial contributions use cumulant-generating-function saddlepoint mass approximations with boundary extensions. Gamma and exponential contributions use a Daniels-type saddlepoint density. Numerical response derivatives use five-point finite differences, with fallback to exact conditional derivatives when the saddlepoint calculation becomes singular.

A requestable `"skew_corrected_laplace"` method is not part of the 0.9.3 public approximation registry.

#### Penalized quasi-likelihood

The PQL path iteratively builds a working response and Fisher-weighted Gaussian mixed-model problem, solves it through the native Gaussian machinery, and applies damping or line search to update the latent mode. Convergence is assessed from changes in the mode and working objective.

PQL is an approximation-specific criterion rather than a marginal likelihood. PQL values must not be compared to Laplace or quadrature log-likelihoods as though they had the same definition, and AIC/BIC are unavailable for PQL fits.

#### Fixed-effect and serial penalties

MEMWAS applies the fixed-effect penalty

\[ P(\beta) = \lambda_1\sum_{k\notin\mathcal I}|\beta_k| + \frac{\lambda_2}{2}\sum_{k\notin\mathcal I}\beta_k^2, \]

where \(\mathcal I\) contains the intercept column. The optimization minimizes the negative selected objective plus this penalty. The fitted object retains:

- adjusted coefficients from the requested penalized fit, used for fitted values and predictions; and
- original coefficients from an unpenalized refit of the selected model, when that refit is available.

`control$autocor_regularization` can add L1, L2, or elastic-net regularization to serial correlation/shape parameters. It does not penalize the serial variance itself. For an unstructured serial covariance, the regularization targets off-diagonal Cholesky parameters.

# How does the MEMWAS work?

## Data requirements and identifiability

Before fitting, verify that:

1. the primary `id` identifies independent top-level clusters;
2. time values are finite, meaningful, and correctly ordered within primary clusters;
3. duplicate times have a defensible interpretation for the requested serial structure;
4. all random-effect grouping variables and design columns are correctly specified;
5. the fixed-effect and random-effect design matrices have adequate rank;
6. the panel contains enough primary clusters and repeated measurements for the proposed covariance complexity;
7. predictor-specific serial loading columns contain useful variation;
8. discrete-grid structures have compatible scaled integer lags where required;
9. an unstructured serial covariance has a sufficiently small and well-supported common set of time levels; and
10. the missing-data assumptions required by the fitted likelihood are scientifically plausible.

A correlation estimate near a boundary can reflect genuine persistence, but it can also indicate inadequate detrending, an omitted random effect, competing latent components, sparse time patterns, or weak identification. Numerical convergence and statistical identification are related but not identical.

## Workflows

A typical MEMWAS workflow is:

1. inspect `MEMWAS_capabilities()` for the installed fitting contract;
2. define the fixed, random, and serial parts of the model;
3. fit the model directly with `fit_MEMWAS()`;
4. inspect `summary()` and `diagnose_approximation()`;
5. compare plausible approximation or covariance choices when scientifically justified;
6. optionally run nonlinear screening, assumption-oriented screens, grouped tuning, autocorrelation ranking, or cluster bootstrap inference; and
7. state the prediction target explicitly when calling `predict()`.

For Gaussian identity models, the term serial residual process can be interpreted as an additive latent Gaussian component in the response model. For non-Gaussian models, the serial process is on the link scale and must not be described as an ordinary additive observed-scale residual.

## Core functions

Function | Implemented role | Typical result
--- | --- | ---
`MEMWAS_capabilities()` | Report the registered families, links, approximations, kernels, covariance structures, and inference contract | Capability table
`serial_component()` | Define one named outcome-loaded or predictor-loaded serial component | `MEMWAS_serial_component`
`fit_MEMWAS()` | Validate, assemble, optimize, diagnose, and optionally screen or bootstrap one model | `MEMWAS_fit`
`diagnose_approximation()` | Expose requested/used approximation, kernel, objective, convergence, Hessian, boundary, and latent-dimension diagnostics | Diagnostic list/table
`compare_approximations()` | Refit one model specification under several active approximation choices | Comparison object and component fits
`screen_MEMWAS_nonlinearity()` | Screen eligible numeric predictors using restricted cubic splines and one of two implemented procedures | Screening object
`check_MEMWAS_assumptions()` | Run selected package-specific diagnostic screens and report unavailable procedures explicitly | Diagnostic-screen object
`rank_MEMWAS_autocorrelation()` | Rank candidate serial structures by grouped prediction or compatible likelihood criteria | Ranking object
`rank_temporal_autocorrelation()` | Compatibility entry point for temporal autocorrelation ranking | Ranking object
`predict()` | Produce conditional, zero-latent, population-marginal, or new-cluster predictions | Numeric vector or interval table

## Key parameters

The main direct interface is:

    fit <- fit_MEMWAS(
      formula,
      family = "gaussian",
      data,
      id,
      time,
      random = ~ 1,
      autocor = "AR(1)",
      serial = NULL,
      predictor_autocor = NULL,
      autocor_predictors = NULL,
      L1_penalty = 0,
      L2_penalty = 0,
      control = list(),
      method = "ML",
      random_cov = "unstructured",
      approximation = "laplace",
      init_approximation = "variational_inference",
      quadrature_points = 7L,
      se_method = "hessian",
      link = NULL,
      offset = NULL,
      weights = NULL,
      theta = NULL,
      shape = NULL,
      screen_nonlinear = FALSE,
      check_assumptions = FALSE,
      bootstrap_inference = FALSE,
      prediction_inference = FALSE,
      verbose = TRUE
    )

Parameter | Implemented meaning and guidance
--- | ---
`formula` | Two-sided fixed-effect formula. Standard formula transformations and factor coding are used to form \(X\beta\).
`family`, `link` | Conditional distribution and supported link. Unsupported family-link combinations are rejected.
`data` | Long-format analysis data. Each row supplies one response and the required predictor, grouping, time, offset, and weight values.
`id` | Name of the primary independent-cluster column. Serial processes are constructed within this clustering level.
`time` | Name of the measurement-time column used for distances or discrete lags.
`random` | One-sided formula, structured term, list of terms, or `NULL`. `~ 1` means a random intercept at the primary `id`; bar syntax supports crossed or nested random terms.
`random_cov` | `"diagonal"` or `"unstructured"`, globally or by term. The executable function default in 0.9.3 is `"unstructured"`.
`autocor` | Convenient scalar or named serial specification. The default is one outcome-loaded `"AR(1)"` process. Use `NULL` when `serial` supplies all components.
`serial` | One `serial_component()` or a list of components. Each component has one loading column and its own covariance parameters.
`predictor_autocor`, `autocor_predictors` | Compatibility interface for adding independent predictor-loaded serial processes.
`L1_penalty`, `L2_penalty` | Non-negative penalties on non-intercept fixed effects.
`control` | Numerical and structural controls, including optimizer limits/tolerances, inner-mode controls, number of starts, Hessian strategy, quadrature limits, VI/PQL limits, serial regularization, time scaling, and permission for negative AR(1).
`method` | `"ML"`, or `"REML"` for an unpenalized Gaussian identity model in the supported exact-Gaussian path.
`approximation` | `"auto"`, `"laplace"`, `"saddlepoint"`, `"adaptive_gaussian_quadrature"`, `"variational_inference"`, or `"pql"`.
`init_approximation` | Supported initialization route. It may differ from the final approximation and does not redefine the final reported objective.
`quadrature_points` | Positive Gauss–Hermite order used by adaptive quadrature and relevant predictive-density integrations.
`se_method` | `"hessian"` or `"none"`. Subject-cluster case bootstrap is activated separately through `bootstrap_inference = TRUE`.
`weights` | Non-negative log-likelihood weights. For binomial proportions, positive integer values are interpreted as trial totals.
`theta`, `shape` | Optional fixed NB2 size or Gamma shape. `NULL` requests estimation when applicable.
`screen_nonlinear` | Run restricted-cubic-spline screening before the final fit. Disabled by default.
`check_assumptions` | Run selected package-specific diagnostic screens after fitting. Disabled by default.
`bootstrap_inference` | Resample primary IDs as whole clusters and refit. This is the implemented dependence-robust resampling route.
`post_selection_inference` | `"none"` or `"bootstrap_after_selection"`. The latter repeats L1 active-set selection in cluster bootstrap refits but remains conditional on the screened formula.
`prediction_inference` | Calculate fitted-row uncertainty after fitting using the requested interval settings.
`...` | No additional unnamed interface is implemented; unused arguments are rejected. In particular, `engine` and settings-object arguments are not accepted.

A fully explicit multi-component serial model can be declared as follows:

    serial_spec <- list(
      outcome = serial_component(
        structure = "OU",
        name = "outcome_persistence"
      ),
      exposure = serial_component(
        structure = "expOU",
        predictor = "exposure",
        name = "exposure_persistence"
      )
    )

    fit_multi <- fit_MEMWAS(
      outcome ~ time_c + exposure + treatment,
      family = "gaussian",
      data = dat,
      id = "subject_id",
      time = "visit_day",
      random = ~(1 + time_c | subject_id),
      autocor = NULL,
      serial = serial_spec
    )

## Approximation diagnostics and uncertainty

After a latent non-Gaussian or nonlinear fit, use:

    diagnose_approximation(fit)

The returned diagnostics distinguish:

- the requested and actually used approximation;
- the native kernel and objective type;
- whether a true marginal log-likelihood is available;
- whether the exact-Gaussian route was used;
- whether fixed effects were analytically profiled;
- the number and maximum dimension of connected latent components;
- whether a dense serial fallback was required;
- optimizer, gradient, step, objective-change, and KKT criteria;
- inner-mode and outer-parameter convergence;
- Hessian and rank diagnostics;
- boundary flags and iteration-limit status; and
- strict versus operational convergence.

The implemented uncertainty routes have different meanings:

Route | Interpretation
--- | ---
`se_method = "hessian"` | Model-based covariance from the observed native Hessian of the fitted objective. Its interpretation depends on the selected likelihood or approximation and model specification.
`se_method = "none"` | Skip covariance and standard-error calculation.
`bootstrap_inference = TRUE` | Primary-ID cluster case-resampling. It propagates refitting variability under resampling of independent clusters and can better display nonlinearity or boundary behavior, subject to bootstrap success requirements.
`post_selection_inference = "bootstrap_after_selection"` | Repeat L1 active-set selection inside cluster bootstrap refits. This addresses active-set variability, not every source of earlier model-selection uncertainty.

A sandwich covariance estimator is deliberately outside this version’s registered inference contract.

## Within-subject correlation structures

### Overview on autocorrelation structures in MEMWAS

Each ordinary serial component has a Gaussian driver covariance \(K_c\). Its diagonal scale parameter is represented through a log standard deviation, ensuring positivity. The remaining raw parameters are transformed to satisfy the structure-specific constraints.

Structure | Implemented covariance or recursion | Native computational path and conditions
--- | --- | ---
`"NONE"` | No serial latent variables | No serial factorization
`"AR(1)"` | \(K_{jk}=\sigma_u^2\rho^{d_{jk}}\) | Exact Markov precision. By default \(\rho=\operatorname{logit}^{-1}(\alpha)\in(0,1)\). With negative autocorrelation allowed, \(\rho=\tanh(\alpha)\) and scaled lags must be integer.
`"OU"` | \(K_{jk}=\sigma_u^2\exp(-d_{jk}/r)\), \(r>0\) | Exact OU Markov precision for ordered distinct times; \(r=\exp(\alpha_r)\). Suitable for irregular elapsed times.
`"expOU"` | OU covariance for Gaussian driver; loading is \(\exp(u-\sigma_u^2/2)-1\) | OU driver Markov precision plus nonlinear transformed loading in likelihood/approximation calculations.
`"AR(p)"` | \(u_t=\sum_{r=1}^{p}\phi_ru_{t-r}+a_t\) | Raw parameters are transformed to reflection/partial-autocorrelation coefficients by `tanh`, then recursively converted to a stationary AR coefficient vector. State-space precision on a regular grid; exact dense fallback when needed.
`"ARMA(1,1)"` | \(u_t=\phi u_{t-1}+a_t+\vartheta a_{t-1}\) | \(\phi=\tanh(\alpha_\phi)\), \(\vartheta=\tanh(\alpha_\vartheta)\). Durbin–Levinson innovations on a regular complete grid; dense exact fallback otherwise.
`"CS"` | Diagonal \(\sigma_u^2\), off-diagonal \(\sigma_u^2\rho\) | Analytical precision/log determinant. The transformed \(\rho\) is constrained above \(-1/(m_{\max}-1)\) and below 1 for the largest relevant cluster.
`"TOEP"` | \(K_{jk}=\sigma_u^2r_{|j-k|}\) | One reflection coefficient per observed lag, transformed by `tanh`; Schur recursion maps these to a positive-definite ACF. Durbin–Levinson innovations on regular grids, dense fallback otherwise.
`"UN"` | Free positive-definite covariance over global fitted time levels | Lower-triangular Cholesky factor with exponentiated diagonal. Requires unique within-cluster times and is limited by `max_unstructured_times` (default 12).

For the ARMA(1,1) implementation, the autocorrelation at lag one is

\[ \rho_1 = \frac{(\phi+\vartheta)(1+\phi\vartheta)} {1+\vartheta^2+2\phi\vartheta}, \]

and for \(k\ge 2\),

\[ \rho_k=\rho_1\phi^{k-1}. \]

For `TOEP`, unconstrained raw parameters \(\alpha_1,\ldots,\alpha_p\) are converted to reflection coefficients

\[ \kappa_j=\tanh(\alpha_j), \]

and the Schur recursion constructs an admissible autocorrelation sequence. This is not a direct unconstrained estimation of separate raw lag correlations.

### Unequally spaced measurements

`OU` is the explicit continuous-time structure:

\[ \operatorname{Corr}(u_{ij},u_{ik}) = \exp\left[-\frac{|t_{ij}-t_{ik}|}{s_t r}\right], \]

where `serial_time_scale` is \(s_t\) and \(r\) is the estimated range.

The AR(1) implementation also uses scaled time distance,

\[ \operatorname{Corr}(u_{ij},u_{ik}) = \rho^{|t_{ij}-t_{ik}|/s_t}, \]

when \(\rho>0\). Thus it can represent fractional positive powers. AR(p), ARMA(1,1), and Toeplitz structures require scaled times to lie on a common integer lag grid. The time unit and `serial_time_scale` must therefore be reported because they determine the interpretation of every correlation or range parameter.

## Model comparison and nonlinear screening

Likelihood-based comparison requires a common response set and compatible objective definitions. AIC and BIC are reported only when the fitted objective exposes a marginal log-likelihood. They are not available for VI or PQL.

For a regular nested likelihood comparison,

\[ \Lambda = 2\{\ell(\widehat\psi_{\mathrm{full}}) - \ell(\widehat\psi_{\mathrm{reduced}})\}. \]

The usual chi-square calibration can fail when the null places a variance or correlation parameter on a boundary. A cluster or parametric bootstrap calibration should be considered when that boundary is central to the scientific conclusion.

The package uses

\[ \operatorname{AIC}=-2\ell(\widehat\psi)+2k \]

and

\[ \operatorname{BIC}=-2\ell(\widehat\psi)+k\log(N_{\mathrm{eff}}), \]

subject to the fitted object’s definition of effective sample size. In clustered data, the analyst should state the independence unit underlying the interpretation of BIC.

`screen_MEMWAS_nonlinearity()` builds restricted cubic spline alternatives. Two procedures are implemented:

1. Nuisance-adjusted score screening. It fits an unpenalized ML-Laplace null, forms nuisance-adjusted conditional score blocks, and uses one shared primary-cluster Rademacher multiplier bootstrap to obtain maxT familywise-adjusted p values. This route requires all random grouping factors to be strictly nested within the primary `id` and uses at least 99 bootstrap replicates.
2. Likelihood-ratio screening. It refits each spline alternative by ML-Laplace and applies the selected p-value adjustment (`holm`, `hochberg`, `hommel`, `bonferroni`, `BH`, `BY`, or `none`). Boundary and repeated-selection caveats still apply.

Because screening is data-dependent model selection, ordinary Hessian standard errors after screening are conditional on the selected formula. They do not automatically include selection uncertainty.

`rank_MEMWAS_autocorrelation()` can compare candidate serial structures by grouped cross-validation or by AIC, BIC, or log-likelihood when those criteria are available and comparable. Grouped folds again preserve primary IDs as indivisible units.

## Predictions and random-effect estimates

`predict.MEMWAS_fit()` supports four explicit targets:

`mode` | Target
--- | ---
`"fitted_cluster_conditional"` | Uses fitted random-effect and serial conditional modes for retained clusters.
`"zero_random_effect"` | Sets random and serial contributions to zero and returns the fixed-plus-offset prediction.
`"population_marginal_mean"` | Integrates over the fitted zero-mean latent distributions and returns a response-scale marginal mean.
`"new_cluster_predictive_distribution"` | Integrates over new latent effects and conditional response variation to represent a new-cluster response distribution.

For a linear Gaussian identity model, a fitted-cluster conditional value is

\[ \widehat y_r^{\mathrm{cond}} = o_r+x_r^\top\widehat\beta + \sum_h z_{rh}^\top\widehat b_{h,g_h(r)} + \sum_c s_{rc}\widehat u_{i(r)c,r}. \]

The Gaussian conditional modes are BLUP-type estimates under the fitted model. In nonlinear models, fitted random and serial values are empirical-Bayes/conditional-mode estimates under the chosen approximation.

For nonlinear links, the population-marginal response mean is

\[ \widehat\mu_r^{\mathrm{marg}} = E_{a\mid\widehat\psi} \left[ g^{-1}\{\eta_r(a)\} \right]. \]

The prediction implementation accumulates the relevant latent variance, including \(\exp(\sigma_c^2)-1\) for an `expOU` transformed coefficient, and uses native Gauss–Hermite moment integration for nonlinear response means. New-cluster predictive bounds are moment-matched bounds, clipped to the response support where necessary; they should not be described as exact marginal quantiles.

For `newdata`, fitted serial states are not available. Unseen random grouping levels are assigned zero fitted-mode contribution when `allow_new_levels = TRUE`; population or new-cluster modes should be used when integration over latent variation is the intended target.

# Examples

## Minimal educational example

The following simulation has a student-specific random intercept, an AR(1) latent serial process, and independent Gaussian observation error.

    set.seed(2026)

    m <- 80L
    n_time <- 4L
    rho <- 0.45
    sigma_u <- 3
    sigma_eps <- 2
    sigma_b <- 5

    student_dat <- data.frame(
      student_id = rep(seq_len(m), each = n_time),
      semester = rep(seq_len(n_time), times = m)
    )

    student_dat$study_minutes <- rnorm(m * n_time, 90, 25)
    student_dat$attendance_rate <- runif(m * n_time, 0.70, 1.00)
    b0 <- rnorm(m, 0, sigma_b)

    simulate_ar1 <- function(n, rho, sigma) {
      u <- numeric(n)
      u[1L] <- rnorm(1L, sd = sigma)
      innovation_sd <- sigma * sqrt(1 - rho^2)
      if (n > 1L) {
        for (tt in 2L:n) {
          u[tt] <- rho * u[tt - 1L] + rnorm(1L, sd = innovation_sd)
        }
      }
      u
    }

    u <- unlist(
      lapply(seq_len(m), function(i) simulate_ar1(n_time, rho, sigma_u)),
      use.names = FALSE
    )

    eps <- rnorm(m * n_time, 0, sigma_eps)

    student_dat$reading_score <-
      55 +
      0.08 * student_dat$study_minutes +
      18 * student_dat$attendance_rate +
      b0[student_dat$student_id] +
      u + eps

    fit_student <- fit_MEMWAS(
      reading_score ~ study_minutes + attendance_rate,
      family = "gaussian",
      data = student_dat,
      id = "student_id",
      time = "semester",
      random = ~ 1,
      autocor = "AR(1)",
      method = "REML",
      approximation = "laplace",
      verbose = FALSE
    )

    summary(fit_student)
    diagnose_approximation(fit_student)

Because this is a linear Gaussian identity-link model, the Laplace request is dispatched to the exact-Gaussian kernel. The random intercept represents persistent between-student level differences. The AR(1) process represents serial variation remaining around those student-specific levels, and the Gaussian family variance represents independent observation error.

## Minimal healthcare example

For irregular follow-up times, use `OU` rather than a nonexistent `continuous_time` switch.

    set.seed(2026)

    m <- 120L
    visit_times <- c(0, 2, 5, 9, 14, 20)
    n_time <- length(visit_times)

    clinic_dat <- data.frame(
      patient_id = rep(seq_len(m), each = n_time),
      visit_week = rep(visit_times, times = m)
    )

    clinic_dat$visit_week_c <-
      clinic_dat$visit_week - mean(visit_times)
    clinic_dat$therapy_sessions <- rpois(m * n_time, lambda = 2)
    clinic_dat$sleep_hours <- rnorm(m * n_time, mean = 6.5, sd = 1.1)

    # Illustrative response only; a realistic simulation would also generate
    # the stated random-slope and OU latent processes.
    clinic_dat$depression_score <-
      14 -
      0.12 * clinic_dat$visit_week_c -
      0.30 * clinic_dat$therapy_sessions -
      0.60 * clinic_dat$sleep_hours +
      rnorm(m * n_time, sd = 4)

    fit_clinic <- fit_MEMWAS(
      depression_score ~ visit_week_c + therapy_sessions + sleep_hours,
      family = "gaussian",
      data = clinic_dat,
      id = "patient_id",
      time = "visit_week",
      random = ~(1 + visit_week_c | patient_id),
      random_cov = "unstructured",
      autocor = "OU",
      control = list(serial_time_scale = 1),
      method = "REML",
      verbose = FALSE
    )

    summary(fit_clinic)

The random intercept and slope represent persistent differences in patient trajectories. The OU process represents remaining serial dependence as a function of elapsed weeks. Changing `serial_time_scale` changes the time unit used to interpret the fitted OU range.

## Non-Gaussian example and approximation check

Adaptive tensor quadrature is feasible only for small connected latent dimensions. The following negative-binomial example uses a random intercept and no serial process so that each latent component is one-dimensional.

    set.seed(2026)

    m <- 100L
    n_time <- 4L
    count_dat <- data.frame(
      patient_id = rep(seq_len(m), each = n_time),
      visit = rep(seq_len(n_time), times = m)
    )

    count_dat$visit_c <- count_dat$visit - mean(seq_len(n_time))
    count_dat$therapy_sessions <- rpois(m * n_time, 2)
    count_dat$sleep_hours <- rnorm(m * n_time, 6.5, 1.1)
    b0 <- rnorm(m, 0, 0.45)

    eta <-
      -0.4 -
      0.08 * count_dat$visit_c +
      0.10 * count_dat$therapy_sessions -
      0.06 * count_dat$sleep_hours +
      b0[count_dat$patient_id]

    count_dat$emergency_visits <- rnbinom(
      nrow(count_dat),
      mu = exp(eta),
      size = 2.5
    )

    fit_count <- fit_MEMWAS(
      emergency_visits ~ visit_c + therapy_sessions + sleep_hours,
      family = "negative_binomial",
      data = count_dat,
      id = "patient_id",
      time = "visit",
      random = ~ 1,
      autocor = "NONE",
      approximation = "laplace",
      init_approximation = "variational_inference",
      se_method = "hessian",
      verbose = FALSE
    )

    summary(fit_count)
    diagnose_approximation(fit_count)

    approx_comparison <- compare_approximations(
      approximations = c(
        "laplace",
        "adaptive_gaussian_quadrature",
        "variational_inference",
        "pql"
      ),
      formula = emergency_visits ~ visit_c + therapy_sessions + sleep_hours,
      family = "negative_binomial",
      data = count_dat,
      id = "patient_id",
      time = "visit",
      random = ~ 1,
      autocor = "NONE",
      quadrature_points = 7L
    )

    print(approx_comparison)

For this log-link model, exponentiated fixed effects are conditional mean ratios. Agreement of coefficients across methods supports numerical sensitivity assessment, but objective values should be compared only within compatible objective types. In particular, VI reports an ELBO and PQL reports a quasi-likelihood criterion, so their objective values, AIC, and BIC must not be ranked against marginal-likelihood rows.

## Predictor-specific serial coefficient example

A named serial component can model a time-varying exposure coefficient:

    fit_exposure <- fit_MEMWAS(
      outcome ~ time_c + exposure + treatment,
      family = "gaussian",
      data = dat,
      id = "subject_id",
      time = "visit_day",
      random = ~(1 | subject_id),
      autocor = NULL,
      serial = list(
        outcome_persistence = serial_component(
          structure = "OU",
          name = "outcome_persistence"
        ),
        exposure_effect = serial_component(
          structure = "expOU",
          predictor = "exposure",
          name = "exposure_effect"
        )
      )
    )

The explicit `serial` list contains an outcome-loaded OU process and a separate predictor-loaded `expOU` process. In version 0.9.3, supply either a non-`NULL` `autocor` specification or `serial`, not both; therefore `autocor = NULL` is required here.

# Reporting recommendations

A reproducible MEMWAS analysis should report:

- the MEMWAS version and the capability registry version;
- the response family, link, response encoding, weights, and offsets;
- the fixed-effect formula and categorical-variable coding;
- every random-effect term, grouping factor, and term-specific covariance parameterization;
- the primary independent `id`;
- the time variable, its unit, and each serial time scale;
- every serial component, its loading variable/design, and its covariance structure;
- whether a component is ordinary Gaussian or `expOU` transformed;
- ML or REML and any L1, L2, or serial regularization;
- the requested and used approximation and native kernel;
- the initialization route and quadrature order, where relevant;
- the standard-error or cluster-bootstrap route;
- operational and strict convergence, gradient/KKT status, Hessian status, rank diagnostics, boundary flags, and dense fallback use;
- whether nonlinear screening, tuning, covariance ranking, directional tests, assumption screens, or post-selection bootstrap was used;
- the prediction mode, response/link scale, included latent contributions, and treatment of new clusters;
- sensitivity analyses across scientifically plausible covariance or approximation choices; and
- the number of primary clusters, observations per cluster, time pattern, and missing-data pattern.

Diagnostic screens should be reported by their MEMWAS method names. Several are intentionally package-specific screens rather than canonical textbook tests. For example, `WithinClusterDifferenceRatioScreen` is not the Durbin–Watson test, `PooledWithinClusterPortmanteauScreen` is not the canonical Ljung–Box test, `PearsonResidualJarqueBeraScreen` is not a randomized-quantile-residual test, and `SquaredResidualFittedLinearScreen` is not the full Breusch–Pagan procedure. A non-significant screen does not prove an assumption.

# References

Breslow, N. E., & Clayton, D. G. (1993). Approximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88(421), 9–25.

Daniels, H. E. (1954). Saddlepoint approximations in statistics. The Annals of Mathematical Statistics, 25(4), 631–650.

Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963–974.

Littell, R. C., Pendergast, J., & Natarajan, R. (2000). Modelling covariance structure in the analysis of repeated measures data. Statistics in Medicine, 19(13), 1793–1819.

Liu, Q., & Pierce, D. A. (1994). A note on Gauss–Hermite quadrature. Biometrika, 81(3), 624–629.

Ormerod, J. T., & Wand, M. P. (2010). Explaining variational approximations. The American Statistician, 64(2), 140–153.

Tierney, L., & Kadane, J. B. (1986). Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association, 81(393), 82–86.