MEMWAS
Mixed-Effects Models with Autocorrelation Structures
Fits longitudinal mixed-effects models through a registered 'C++'
numerical backend. Supported serial covariance structures include first-order
autoregressive (AR(1)), exponential or Ornstein-Uhlenbeck, higher-order
autoregressive (AR(p)), first-order autoregressive moving-average (ARMA(1,1)),
compound symmetry, Toeplitz, and unstructured covariance. Serial processes
can be unified or attached independently to numeric predictor loadings.
Candidate temporal structures can be ranked on a common sample by
dependence-component grouped cross-validation, the Akaike information
criterion, the Bayesian information criterion, or log-likelihood. Clustered,
crossed, and nested random intercepts and slopes are assembled jointly with
diagonal or term-specific unstructured covariance. Available approximation
methods include Laplace, saddlepoint likelihood with latent Laplace
integration, adaptive Gaussian quadrature, full-covariance Gaussian variational
inference, and penalized quasi-likelihood. Penalized smooth mean terms
include ordinary and cyclic P-splines, factor-by and varying-coefficient terms,
tensor products, shrinkage smooths, and whole-term selection. Term-specific
penalties, grouped fold-local smoothing selection, null-space
constraints, and smooth effective degrees of freedom remain separate from
elastic-net coefficient shrinkage while the smooth mean and serial covariance
are fitted jointly. Bootstrap resampling preserves the declared dependence
components. The mixed-effects framework is inspired by Laird and Ware (1982)
MEMWAS
|
Mixed Effects Models With Autocorrelation Structures |
|
Version 0.9.5 Aug. 22, 2026 |
MEMWAS is a base-R package for fitting longitudinal mixed-effects models for Gaussian, Bernoulli/grouped-binomial, Poisson, negative-binomial, Gamma, and exponential outcomes. Models may combine fixed effects, Gaussian random effects, offsets, fixed-effect penalties, penalized smooth mean terms, and multiple independent named serial processes in one likelihood. Native smooths include ordinary and cyclic P-splines, factor-by/treatment-specific and varying-coefficient terms, tensor products, shrinkage smooths, and grouped-CV whole-term selection. Each smooth has its own penalty and smoothing parameter, separate from elastic-net coefficient shrinkage, and is estimated jointly with the random and serial covariance components. Each serial process has a one-column design and its own AR(1), exponential/Ornstein-Uhlenbeck, AR(p), ARMA(1,1), compound-symmetry, Toeplitz, or unstructured covariance, permitting residual and predictor-modulated autocorrelation structures to be estimated simultaneously. Grouped-binomial totals, family links, conditional predictions, and serial-component contributions are retained consistently through fitting and prediction. Besides, the MEMWAS package offers an integrated transparency advantage by embedding configurable assumption screening and structured diagnostic reporting within the model- analysis workflow. This design supports research integrity by making diagnostic choices, findings, limitations, and unavailable tests more visible and auditable.
Environment
- R version R >= 4.1.0
- RTools version 4.5.6768
- C++ 11
Installation
# Install the package from a local source archive:
install.packages("MEMWAS_0.9.5.tar.gz", repos = NULL, type = "source")
Workflow
library(MEMWAS)
set.seed(1L)
sim_data <- MEMWAS:::.simulate_panel_data(
n_id = 50L, n_time = 3L,
beta = c(x1 = 0.7, x2 = -0.4, x3 = 0.2),
cor_matrix = diag(3L), intercept = 0.8,
sigma_eps = 0.5, sigma_b = 0.4, autocor = "AR(1)",
autocor_param = list(rho = 0.35)
)
Ranking temporal autocorrelation structures
ranked <- rank_autocorrelation_structures(
y ~ x1 + x2,
data = sim_data,
id = "id",
time = "time",
random = ~ (1 | id),
candidates = list(
independent = "NONE",
residual_ar1 = list(residual_autocor = "AR(1)"),
predictor_specific = list(
residual_autocor = "AR(1)",
predictor_autocor = c(x1 = "OU", x2 = "AR(1)")
)
),
criterion = "grouped_cv",
K = 5L,
metric = "RMSE",
verbose = FALSE
)
ranked$ranking
ranked$selected
criterion accepts "grouped_cv", "AIC", "BIC", and "logLik".
Grouped CV assigns each connected dependence component wholly to one fold.
Components are induced jointly by the primary ID and every nested or crossed
random-effect grouping factor, so shared latent levels cannot leak between
training and validation rows. One common fold assignment is used for all
candidates, and every fold must fit, predict, and return a finite metric. AIC,
BIC, and log-likelihood comparisons require an available likelihood criterion.
All viable candidates are refitted on the same complete-case rows, so the
ranking never compares different analysis samples.
The selected candidate can be refitted automatically with
refit_best = TRUE.
For large random-effects designs, control = list(dense_fallback = FALSE)
forbids compatibility materialization of the global dense design. The native
pivoted rank diagnostic continues to use compressed row pointers, column
indices, and nonzero values and preserves its pivot and aliased-column report.
Do not combine that prevention request with block_factorization = FALSE,
which explicitly selects the dense reference path.
Standalone nonlinearity screening
baseline <- fit_MEMWAS(
y ~ x1 + x2,
data = sim_data,
id = "id",
time = "time",
random = ~ (1 | id),
autocor = "NONE",
se_method = "none",
verbose = FALSE
)
screen <- screen_MEMWAS_nonlinearity(
baseline,
nonlinear_predictors = c("x1", "x2"),
nonlinear_screening_method = "nuisance_adjusted_score",
nonlinear_bootstrap_reps = 499L,
verbose = FALSE
)
screen$summary
screen$selected_formula
The default score procedure uses one unpenalized ML-Laplace null model,
nuisance-adjusted spline score blocks, and a shared primary-cluster maxT
multiplier bootstrap. It requires random-effect grouping factors to be nested
within the primary ID. The "likelihood_ratio" option rebuilds and refits each
candidate model and can be used with crossed grouping structures.
Integrated screening is also available through screen_nonlinear = TRUE in
fit_MEMWAS().
Select the native P-spline screening engine explicitly when desired:
pscreen <- screen_MEMWAS_nonlinearity(
baseline,
nonlinear_predictors = c("x1", "x2"),
nonlinear_spline = "pspline",
spline_knots = 10L,
smooth_control = list(cv_folds = 5L, metric = "RMSE"),
nonlinear_bootstrap_reps = 499L,
verbose = FALSE
)
The P-spline screen reports penalty-adjusted effective degrees of freedom, not
the number of basis columns. Selected smooths retain automatic, term-specific
smoothing parameters for the final joint mean/covariance fit. The engine is not
run when screen_nonlinear = FALSE, and run_checks_and_screening = FALSE
suppresses integrated screening regardless of engine choice.
Approximation and prediction tools
MEMWAS_capabilities()
diagnose_approximation(fit)
comparison <- compare_approximations(
c("laplace", "variational_inference"),
formula = y ~ x1 + x2,
data = sim_data,
id = "id",
time = "time"
)
Available approximation requests include Laplace, saddlepoint likelihood with latent Laplace integration, adaptive Gaussian quadrature, full-covariance Gaussian variational inference, and penalized quasi-likelihood. The capability table reports valid family/link, random-effect, serial, estimation, inference, and prediction combinations.
Prediction modes include fitted-cluster conditional values, zero-random-effect
values, population marginal means, and new-cluster predictive distributions.
Use predict() with interval, level, and mode to request uncertainty when
the fitted approximation supports it.
Fitting MEMWAS
fit <- fit_MEMWAS(
y ~ x1 + x2,
family = "gaussian",
data = sim_data,
id = "id",
time = "time",
random = ~(1 | id),
autocor = "AR(1)",
L2_penalty = 0.1,
verbose = FALSE
)
fit_MEMWAS() supports Gaussian, binomial, Poisson, negative-binomial, gamma,
and exponential responses with valid family-specific links. Random effects may
be clustered, crossed, or nested, with diagonal or term-specific unstructured
covariance. Serial processes may be outcome-loaded or attached independently to
numeric predictors. Supported structures include no serial process, AR(1), OU,
expOU, AR(p), ARMA(1,1), compound symmetry, Toeplitz, and unstructured
covariance.
For ordinary serial structures, duplicate rows with the same ID and time share one latent state; compound symmetry remains observation-indexed. Raw zero means independence for the default non-negative AR(1), OU, expOU, and compound- symmetry dependence coordinates.
Separate residual autocorrelation components
The explicit autocorrelation interface separates the outcome-loaded residual process from predictor-loaded residual processes:
fit_components <- fit_MEMWAS(
y ~ x1 + x2,
data = sim_data,
id = "id",
time = "time",
random = ~ (1 | id),
residual_autocor = "AR(1)",
predictor_autocor = c(x1 = "OU", x2 = "CS"),
verbose = FALSE
)
residual_autocor defines an optional outcome-loaded component.
predictor_autocor is a uniquely named character vector or list; each entry
loads an independently parameterized serial covariance component by that
numeric predictor. An entry equal to "NONE" is omitted. Do not combine this
explicit interface with non-NULL autocor or serial arguments.
For primary cluster i, the native covariance contribution can be written as
$$ V_i = K_{0i}(\theta_0) + \sum_j \mathrm{diag}(x_{ij}) K_{ji}(\theta_j) \mathrm{diag}(x_{ij}). $$
Covariance construction, factorization, approximation, likelihood evaluation, and optimization remain in the registered C++ backend. Predictor-loaded components describe conditional residual covariance heterogeneity; they are not joint stochastic time-series models for the predictors themselves.
Penalized smooth mean terms
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"
),
verbose = FALSE
)
summary(fit_smooth)$smooth_summary
For term (j), MEMWAS minimizes the fitted objective plus
$$ \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}. $$
This penalty is distinct from L1_penalty and L2_penalty. Smooth
coordinates stay in the profiled fixed-coefficient block, are excluded from L1
shrinkage, and are fitted jointly with the random effects and serial covariance.
The penalty is diagonalized term by term; its range and null space are separated
and given term-appropriate identifiability constraints. Ordinary, cyclic, and
tensor terms are centered/projected against the parametric design. For factor-by
and varying-coefficient terms, residualization diagnoses aliases only: the
fitted forms remain B(x) * I(level) and z * B(x) and do not acquire
unrequested interactions with unrelated parametric covariates. select = TRUE (or
type = "shrinkage") also penalizes the null space, while
selection = "whole_term" uses grouped CV to compare removal of the complete
term. These two selection modes are mutually exclusive.
Available type values cover "pspline", "cyclic_pspline", "factor_by",
"varying_coefficient", and "tensor_product"; a factor by produces
treatment-specific curves and a numeric by produces a varying coefficient.
Tensor terms accept per-margin k, degree, difference order, boundary, and
smoothing parameters. Term-specific lambda_initial and rank_tolerance
override their smooth_control defaults. Cyclic terms predict periodically;
noncyclic terms default to constant boundary extrapolation and may instead use
their stored "error" or "linear" rule.
Subject-specific random smooths are intentionally deferred until sparse latent
integration is redesigned.
Automatic smoothing is dependence-component grouped and fold-local. Training folds rebuild boundaries, bases, constraints, and any requested nonlinear screen from positive-weight rows, while validation data are transformed through the training blueprint. Only positive-weight rows enter the dependence graph, so a zero-weight row cannot bridge otherwise independent components or create a fold. Each smoothing candidate re-estimates the smooth mean and serial covariance jointly. Summaries report term-specific smoothing parameters, penalty/null-space ranks, effective and reference degrees of freedom, statistics, p-values, and inferential status rather than using raw basis counts.
With any fixed-effect, smooth, or serial-covariance penalty or whole-term smooth-selection step active, ordinary Wald and Hessian/delta uncertainty is suppressed. The dependence-preserving bootstrap repeats retained automatic smoothing and whole-term selection while remaining conditional on nonlinear- screen candidate selection. Coordinate-wise percentile intervals are limited to invariant parametric coefficients. Primary-ID case resampling rebuilds a global smooth projection from each resample, so its nominal parametric coefficients are not comparable. The same suppression applies when whole-term replay can add or remove a global projection. MEMWAS then retains the refit and selection diagnostics but explicitly suppresses that coefficient table in favor of common-scale function or prediction contrasts.
Nonlinearity screening and assumption diagnostics can be disabled by:
formals(fit_MEMWAS)$screen_nonlinear
# FALSE
formals(fit_MEMWAS)$check_assumptions
# FALSE
This keeps ordinary fitting focused on estimation. Each optional stage can be run independently after obtaining a suitable fitted model.
Standalone assumption diagnostics
checks <- check_MEMWAS_assumptions(
fit,
autocorrelation_check = "All",
distribution_link_check = "All",
conditional_independence_check = "All",
random_effects_normality_check = "All",
random_effects_predictor_independence_check = "All",
homogeneity_variance_check = "All"
)
checks$results
checks$flagged
checks$unavailable
The six diagnostic categories are selected independently. Holm adjustment is applied jointly to all requested diagnostics with available finite p-values. Structurally unavailable diagnostics remain in the result with an explanation; they are not treated as evidence that an assumption holds.
Integrated checking is available through check_assumptions = TRUE in
fit_MEMWAS(). run_checks_and_screening = FALSE disables both optional stages,
whereas TRUE activates either stage whose individual flag was omitted.
Original and adjusted fixed-effect coefficients
For a penalized fit, MEMWAS distinguishes two coefficient vectors:
- Original coefficients: estimates from the same selected model refitted without elastic-net L1/L2 penalties; any smooth penalties remain.
- Adjusted coefficients: regularized estimates from the final selected penalized model. These are the coefficients used for fitted values and predictions.
coef(fit, type = "original")
coef(fit, type = "adjusted")
coef(fit, type = "both")
fit$original_equation
fit$final_equation
With g denoting the link, the stored equations correspond to
$$ g{E(Y_{it}\mid b_i,u_i)} = X_{it}\beta^{(0)} + \sum_j f_j(x_{it}) + Z_{it}b_i + S_{it}u_i + o_{it} $$
for the Original Model equation, and
$$ g{E(Y_{it}\mid b_i,u_i)} = X_{it}\widetilde{\beta} + \sum_j f_j(x_{it}) + Z_{it}b_i + S_{it}u_i + o_{it} $$
for the Final equation. Here, beta^(0) is the L1/L2-unpenalized reference vector
and tilde(beta) is the adjusted vector from the penalized fit. Smooths are
shown by term label instead of internal basis-coordinate names. The comparison
table also stores adjusted - original for every coefficient.