--- title: "Introduction to MEMWAS" subtitle: |

Mixed
Effects
Models
With
Autocorrelation
Structures

author: "[Enoch Kang](https://orcid.org/0000-0002-4903-942X)" date: "Aug. 20, 2026" description: > Mixed Effects Models with Autocorrelation Structures in the current release. output: rmarkdown::html_vignette: toc: true toc_depth: 3 keep_md: true vignette: > %\VignetteIndexEntry{Introduction to MEMWAS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE ) # library(MEMWAS) ```
# Installation Install the package from its 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. ```{r install-tarball, eval=FALSE} install.packages( "MEMWAS.tar.gz", repos = NULL, type = "source" ) #library(MEMWAS) ``` The capability registry is the most direct way to inspect the fitting contract implemented by the installed version: ```{r capabilities, eval=FALSE} MEMWAS_capabilities() ``` > **Release scope.** This vignette documents the public interface and numerical paths in the installed package. Models are fitted directly through `fit_MEMWAS()`.
# 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, autocorrelation-structure ranking, automatic dependence-aware 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 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. ## Penalized smooth mean functions MEMWAS provides a dependency-free P-spline engine alongside restricted cubic splines. An explicitly specified smooth augments the conditional mean as \[ g\{E(Y_{it}\mid b_i,u_i)\} =X_{it}\beta+\sum_j B_j(x_{it})\gamma_j +Z_{it}b_i+S_{it}u_i+o_{it}. \] For term $j$, its contribution is \[ \frac{1}{2}\lambda_{s,j}\gamma_j^\mathsf{T}S_j\gamma_j, \qquad S_j=D_{d_j}^\mathsf{T}D_{d_j} \] for the default coefficient-difference P-spline penalty of order $d_j$. This penalty is distinct from elastic-net shrinkage: each smooth has its own $S_j$ and $\lambda_{s,j}$, smooth coordinates are excluded from the L1 penalty, and the smooth penalty is stored and optimized separately from `L1_penalty` and `L2_penalty`. MEMWAS diagonalizes each term penalty, partitions its penalized range from its null space, applies term-appropriate alias and identifiability constraints, and retains the identifiable smooth coordinates in the profiled fixed-coefficient block. Ordinary, cyclic, and tensor terms are centered/projected against the parametric design. For factor-by and varying-coefficient smooths, residualization diagnoses exact aliases only; the fitted forms remain $B(x)I(\mathrm{level})$ and $zB(x)$, so unrelated parametric covariates cannot create unrequested interactions. The smooth mean, random effects, and serial covariance are consequently estimated jointly rather than in a two-stage residual smoother. The native univariate P-spline is the default smooth type: ```{r pspline-fit, eval=FALSE} fit_smooth <- fit_MEMWAS( y ~ x2, data = sim_data, id = "id", time = "time", random = ~ (1 | id), autocor = "AR(1)", smooth = list( list( name = "x1_curve", type = "pspline", variable = "x1", k = 10L, degree = 3L, difference_order = 2L, lambda = "auto" ) ), smooth_control = list( optimizer = "grouped_cv", cv_folds = 5L, log_lambda_range = c(-4, 4), metric = "RMSE" ) ) summary(fit_smooth)$smooth_summary ``` `smooth` may be one term list or a list of term lists. The implemented classes are: - ordinary P-splines (`type = "pspline"`); - cyclic P-splines (`"cyclic_pspline"`), with optional `period` and periodic prediction; - factor-by/treatment-specific smooths (`"factor_by"`, with a factor named by `by`); - varying-coefficient smooths (`"varying_coefficient"`, with a numeric `by`); - tensor products (`"tensor_product"`, with two or more `variables` and margin-specific controls); and - shrinkage smooths (`type = "shrinkage"` or `select = TRUE`) and grouped-CV whole-term selection (`selection = "whole_term"`). Null-space shrinkage and whole-term removal are alternative selection strategies and cannot be requested together. `boundary_knots` controls the fitted domain; term-level `lambda_initial` and `rank_tolerance` override their `smooth_control` defaults. Noncyclic terms default to `extrapolation = "constant"` and may instead use `"error"` or `"linear"`; cyclic terms wrap periodically. Subject-specific random smooths are intentionally rejected until sparse latent integration is redesigned. Set a finite non-negative `lambda` to fix a smoothing parameter, or use `NULL`/`"auto"` to select it internally. Tensor products accept one smoothing parameter per margin. Automatic selection uses connected dependence-component folds among positive-weight rows: observations that share a primary ID or any nested or crossed grouping level stay together, while zero-weight rows cannot bridge active components. Within every fold, boundaries, bases, factor levels, penalty diagonalizations, null-space constraints, and optional nonlinear screening are derived from training rows only. Validation rows are transformed with the stored training blueprint. Every smoothing candidate refits the smooth mean and serial covariance jointly. This fold-local engine is used inside an ordinary final fit, and the augmented fixed design makes smooths available to all supported families. Nonlinearity screening remains opt-in. Choose `nonlinear_spline = "pspline"` only together with an explicit screening request such as `screen_nonlinear = TRUE` or a direct call to `screen_MEMWAS_nonlinearity()`. It is not run when `screen_nonlinear = FALSE`, and `run_checks_and_screening = FALSE` suppresses the integrated screen. Restricted cubic splines remain the default screening engine. For P-splines, screening and final summaries use smooth-specific effective and reference degrees of freedom, statistics, p-values, and inferential status rather than interpreting the raw number of basis columns as degrees of freedom. ## 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. The current implementation stores the global random-effects design canonically as row pointers, zero-based column indices, and nonzero loadings. Its pivoted, rank-revealing diagnostic operates directly on that representation and retains the pivot decision and aliased-column report. Set `control$dense_fallback = FALSE` to forbid compatibility fallback that materializes the global $N\times q$ random-effects design. `control$block_factorization = FALSE` explicitly selects the dense reference path and therefore cannot be combined with this prevention request. 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 executes 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"`, `"native_fixed_effect"`, `"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`, MEMWAS does not label a two-moment Gaussian substitution as an ELBO: the substitution is not generally the actual expected log likelihood, and its covariance update is not the derivative of that substituted objective. A requested variational initialization is therefore replaced by the native fixed-effect start, and a requested variational final fit fails explicitly with guidance to use Laplace or adaptive Gaussian quadrature. For supported variational fits, 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 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 a refit without elastic-net L1/L2 penalties, while retaining any smooth penalty, 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. When any fixed-effect, smooth, or serial-covariance penalty is active, the objective Hessian is retained only as optimization curvature; it is not treated as an ordinary sampling covariance. Wald standard errors, tests, confidence intervals, and Hessian/delta prediction inference are therefore unavailable for the penalized estimate. Dependence-preserving bootstrap refits repeat the complete fitted penalized estimation and active-set operation. Empirical standard errors and percentile intervals are reported only for invariant parametric coefficient coordinates. Primary-ID case resampling rebuilds a global smooth projection from each resample, so its nominal parametric coefficients are resample-local; the same issue arises when whole-term replay can add or remove a global projection. MEMWAS then retains refit and selection-replay diagnostics but explicitly suppresses that coefficient table in favor of common-scale function or prediction contrasts. Bootstrap output does not report sign-count or universal selection-corrected p-values.
# 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, autocorrelation ranking, or dependence-aware 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 or native P-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_autocorrelation_structures()` | Rank candidate serial structures by dependence-component grouped prediction or compatible likelihood criteria | 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: ```{r fit-interface, eval=FALSE} fit <- fit_MEMWAS( formula, family = "gaussian", data, id, time, random = ~ (1 | id), autocor = "AR(1)", serial = NULL, predictor_autocor = NULL, L1_penalty = 0, L2_penalty = 0, smooth = NULL, smooth_control = list(), 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, nonlinear_spline = "restricted_cubic", 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 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` | Interface for adding independently parameterized predictor-loaded serial processes. | | `L1_penalty`, `L2_penalty` | Non-negative penalties on non-intercept fixed effects. | | `smooth` | One smooth term list, a list of term lists, or `NULL`. Supports ordinary/cyclic P-splines, factor-by, varying-coefficient, tensor-product, shrinkage, and whole-term-selection specifications. | | `smooth_control` | Basis defaults and internal grouped-CV controls for term-specific smoothing parameters. | | `control` | Numerical and structural controls, including optimizer limits/tolerances, inner-mode controls, concurrent starts, Hessian strategy, quadrature limits, VI/PQL limits, serial regularization, time scaling, and permission for negative AR(1). `dense_fallback = FALSE` forbids compatibility materialization of the global random-effects design while retaining compressed pivoted rank diagnostics. | | `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"`. Automatic dependence-aware bootstrap inference 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`, `nonlinear_spline` | Optionally screen with `"restricted_cubic"` or `"pspline"` before the final fit. Screening is disabled by default. | | `check_assumptions` | Run selected package-specific diagnostic screens after fitting. Disabled by default. | | `bootstrap_inference` | Select a dependence-preserving bootstrap automatically. Primary-ID case resampling is used only for primary-ID-only grouping and draws only positive-weight rows from IDs having at least one such row; zero-weight rows and zero-weight-only IDs are ignored. Nested or crossed models use joint parametric random-effect, serial-process, and response simulation on the fitted incidence graph. The joint parametric route requires every positive observation/frequency weight to equal one; zero-weight rows are ignored. | | `post_selection_inference` | `"none"` or `"bootstrap_after_selection"`. The latter is a compatibility alias that turns on the same dependence-aware bootstrap as `bootstrap_inference = TRUE`; every penalized bootstrap repeats the fitted active set, automatic smoothing and grouped-CV whole-term selection are replayed from the retained preselection specification, and nonlinear-screen selection remains conditional on the screened formula. | | `prediction_inference` | Calculate fitted-row Hessian/delta uncertainty after an unpenalized fit. It fails explicitly when a fixed-effect, smooth, or serial-covariance penalty is active; use a dependence-preserving bootstrap workflow for penalized uncertainty. | | `...` | 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: ```{r multi-serial, eval=FALSE} set.seed(1L) sim_data <- MEMWAS:::.simulate_panel_data( n_id = 30L, n_time = 4L, beta = c(x1 = 0.6, x2 = -0.3, x3 = 0.2), cor_matrix = diag(3L), intercept = 0.5, sigma_eps = 0.5, sigma_b = 0.4, autocor = "NONE" ) serial_spec <- list( outcome = serial_component( structure = "OU", name = "outcome_persistence" ), exposure = serial_component( structure = "expOU", predictor = "x1", name = "exposure_persistence" ) ) fit_multi <- fit_MEMWAS( y ~ x1 + x2 + x3, family = "gaussian", data = sim_data, id = "id", time = "time", random = ~ (1 + x3 | id), autocor = NULL, serial = serial_spec ) ``` ## Approximation diagnostics and uncertainty After a latent non-Gaussian or nonlinear fit, use: ```{r diagnostics, eval=FALSE} 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` | Automatic dependence-aware bootstrap. It uses primary-ID case resampling only when every grouping partition matches the primary ID; those native draws contain only positive-weight rows from positive-weight-eligible IDs. Otherwise it simulates the fitted joint nested/crossed random, serial, and response model on the original incidence graph. On the joint parametric route, every positive observation/frequency weight must equal one; zero-weight rows remain likelihood-inactive. | | `post_selection_inference = "bootstrap_after_selection"` | Turn on the dependence-aware bootstrap as an alias for `bootstrap_inference = TRUE`. It is not a second refit algorithm: every penalized bootstrap repeats the fitted active set and replays retained automatic-smoothing/whole-term choices, while nonlinear screening remains conditional. Percentile intervals are limited to invariant parametric coordinates; a resample-local global smooth projection suppresses that coefficient table. No sign-count or universal selection-corrected p-value is reported. | 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. Except for compound symmetry, serial states are keyed by distinct `(id, time)` pairs: duplicate rows at the same time share one latent state. Compound symmetry is instead indexed by observations, so distinct equal-time rows have the fitted off-diagonal correlation rather than correlation one. | 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 on integer scaled lags. By default $\rho=\tanh(\alpha)^2\in[0,1)$, so $\alpha=0$ is independence. With negative autocorrelation allowed, $\rho=\tanh(\alpha)$. | | `"OU"` | $K_{jk}=\sigma_u^2\exp(-d_{jk}/r)$ for $r>0$ | Exact OU Markov precision for ordered distinct times; $r=\alpha_r^2$, with $r=0$ defined as independence between distinct times. Suitable for irregular elapsed times. | | `"expOU"` | OU covariance for Gaussian driver; loading is $\exp(u-\sigma_u^2/2)-1$ | The driver range is $r=\alpha_r^2$, with raw zero denoting independence. OU driver Markov precision plus nonlinear transformed loading is used in supported likelihood approximations. | | `"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. Equal-time duplicate rows share the corresponding time-level state; the number of distinct levels 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 penalized likelihood fits, both criteria additionally require effective degrees of freedom from a valid generalized-trace calculation; an active-parameter-count fallback is diagnostic only and is not substituted into either criterion. 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. When inferential bootstrap is requested, MEMWAS automatically uses primary-ID case resampling only when that partition captures every fitted grouping level and otherwise uses joint parametric simulation on the nested or crossed incidence graph. The package uses \[ \operatorname{AIC}=-2\ell(\widehat\psi)+2k \] and \[ \operatorname{BIC}=-2\ell(\widehat\psi)+k\log(N_{\mathrm{eff}}), \] where $N_{\mathrm{eff}}$ is the number of independent connected components among positive-weight rows in the joint dependence graph induced by the primary ID and every nested or crossed random-effect grouping factor. Longitudinal serial dependence is contained within the primary-ID edges. Frequency weights alter likelihood contributions but do not create additional BIC sampling units. MEMWAS requires at least two connected components for BIC; a fully connected crossed design can therefore retain an available AIC while BIC is unavailable. `screen_MEMWAS_nonlinearity()` builds restricted cubic spline alternatives by default or native univariate P-spline alternatives when `nonlinear_spline = "pspline"`. P-spline screening uses a term-specific difference penalty and reports penalty-adjusted effective degrees of freedom; selected terms undergo fold-local grouped smoothing selection in the final fit. 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`). A P-spline comparison uses its smooth effective degrees of freedom rather than its raw basis count. 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. Screening never runs solely because the P-spline engine is available; it requires a direct call or an explicit integrated screening flag. `rank_autocorrelation_structures()` can compare candidate serial structures by grouped cross-validation or by AIC, BIC, or log-likelihood when those criteria are available and comparable. For grouped CV, rows connected through the primary ID or any nested or crossed random-effect grouping level form one dependence component. Components are assigned intact to folds, preventing shared latent levels from leaking between training and validation data. Non-Gaussian targets and predictive-density/deviance metrics fail explicitly because they require integration over the actual shifted-lognormal driver; MEMWAS does not replace it with a two-moment Gaussian distribution. ## 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 integration for nonlinear response means. For Gaussian identity-link models, `joint = TRUE` additionally returns the full cross-row covariance, separated into fixed-parameter, nested or crossed random, serial, and observation components. Except for Gaussian identity-link models, conditional predictive intervals are available only on the response scale and use quadrature of the actual response-family mixture over asymptotic fixed-parameter uncertainty. This includes Gaussian models with a log link. Their `fit` column is the mixture-distribution mean; `plug_in_fit` retains the inverse-link prediction at the fitted coefficients. New-cluster intervals instead integrate the declared latent random and serial variance with the actual response family, respecting discrete or positive response support. MEMWAS does not substitute a two-moment Gaussian distribution for a non-Gaussian `expOU` predictive target; that request fails explicitly. 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. ```{r education-example, eval=FALSE} set.seed(1L) sim_data <- MEMWAS:::.simulate_panel_data( n_id = 30L, n_time = 4L, beta = c(x1 = 0.7, x2 = -0.3), cor_matrix = diag(2L), intercept = 1, sigma_eps = 0.4, sigma_b = 0.5, autocor = "AR(1)", autocor_param = list(rho = 0.45) ) fit_student <- fit_MEMWAS( y ~ x1 + x2, family = "gaussian", data = sim_data, id = "id", time = "time", random = ~ (1 | id), autocor = "AR(1)", method = "REML", approximation = "laplace", control = list(n_starts = 1L, cold_start_verification = FALSE), 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 elapsed-time dependence, request `OU` directly rather than using a nonexistent `continuous_time` switch. ```{r healthcare-example, eval=FALSE} set.seed(2L) sim_data <- MEMWAS:::.simulate_panel_data( n_id = 30L, n_time = 4L, beta = c(x1 = -0.5, x2 = -0.3, x3 = 0.2), cor_matrix = diag(3L), intercept = 1.5, sigma_eps = 0.5, sigma_b = 0.4, autocor = "NONE" ) fit_clinic <- fit_MEMWAS( y ~ x1 + x2 + x3, family = "gaussian", data = sim_data, id = "id", time = "time", random = ~ (1 + x1 | 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 Poisson example uses a random intercept and no serial process so that each latent component is one-dimensional. ```{r nongaussian-example, eval=FALSE} set.seed(3L) sim_data <- MEMWAS:::.simulate_panel_data( n_id = 30L, n_time = 3L, beta = c(x1 = 0.15, x2 = -0.10), cor_matrix = diag(2L), intercept = -0.2, sigma_eps = 0, sigma_b = 0.35, family = "poisson", autocor = "NONE" ) fit_count <- fit_MEMWAS( y ~ x1 + x2, family = "poisson", data = sim_data, id = "id", time = "time", random = ~ (1 | id), 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", "variational_inference"), formula = y ~ x1 + x2, family = "poisson", data = sim_data, id = "id", time = "time", random = ~ (1 | id), 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: ```{r predictor-serial-example, eval=FALSE} set.seed(4L) sim_data <- MEMWAS:::.simulate_panel_data( n_id = 30L, n_time = 4L, beta = c(x1 = 0.6, x2 = -0.3), cor_matrix = diag(2L), intercept = 0.5, sigma_eps = 0.4, sigma_b = 0.3, autocor = "NONE" ) fit_exposure <- fit_MEMWAS( y ~ x1 + x2, family = "gaussian", data = sim_data, id = "id", time = "time", random = ~ (1 | id), autocor = NULL, serial = list( outcome_persistence = serial_component( structure = "OU", name = "outcome_persistence" ), exposure_effect = serial_component( structure = "expOU", predictor = "x1", name = "exposure_effect" ) ) ) ``` The explicit `serial` list contains an outcome-loaded OU process and a separate predictor-loaded `expOU` process. 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 automatic dependence-aware bootstrap route; - operational and strict convergence, gradient/KKT status, Hessian status, rank diagnostics, boundary flags, and dense fallback use; - whether nonlinear screening, 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.