|
Mixed Effects Models With Autocorrelation Structures |
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.
Or install the current development version from GitHub:
The capability registry is the most direct way to inspect the fitting contract implemented by the installed version:
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 formerset_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.
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.
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.
MEMWAS 0.9.3 fits longitudinal mixed-effects models with:
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.
Let \(r=1,\ldots,N\) index analysis
rows and let \(i(r)\) identify the
primary independent cluster supplied through id. Let:
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.
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.
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 | 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.
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.
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; andLet
\[ 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.
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.
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 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.
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.
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.
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:
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.
Before fitting, verify that:
id identifies independent top-level
clusters;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.
A typical MEMWAS workflow is:
MEMWAS_capabilities() for the installed fitting
contract;fit_MEMWAS();summary() and
diagnose_approximation();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.
| 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 |
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
)After a latent non-Gaussian or nonlinear fit, use:
The returned diagnostics distinguish:
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.
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.
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.
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:
id and uses at least 99 bootstrap replicates.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.
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.
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.
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.
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.
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.
A reproducible MEMWAS analysis should report:
id;expOU
transformed;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.
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.