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.