MEMWAS

	Mixed Effects Models With Autocorrelation Structures
Version 0.9.3 Aug. 04, 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, and multiple independent named serial processes in one likelihood. 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.3.tar.gz", repos = NULL, type = "source")

# Or install from GitHub:
remotes::install_github("EnochKang/MEMWAS")

Main fitting workflow

library(MEMWAS)

fit <- fit_MEMWAS(
  y ~ x1 + x2,
  family = "gaussian",
  data = dat,
  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.

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 = dat,
  id = "id",
  time = "time",
  random = ~1,
  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}(0) + j (x{ij}) K{ji}(j) (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.

Ranking temporal autocorrelation structures

ranked <- rank_autocorrelation_structures(
  y ~ x1 + x2,
  data = dat,
  id = "id",
  time = "time",
  random = ~1,
  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 primary ID wholly to one fold, uses one common fold assignment for all candidates, and requires every fold to 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.

Nonlinearity screening and assumption diagnostics are disabled by default:

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 nonlinearity screening

baseline <- fit_MEMWAS(
  y ~ x1 + x2,
  data = dat,
  id = "id",
  time = "time",
  random = ~1,
  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 remains available through screen_nonlinear = TRUE in fit_MEMWAS() or tune_MEMWAS().

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 fixed-effect penalties.
• 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}b_i,u_i)} = X_{it}^{(0)} + Z_{it}b_i + S_{it}u_i + o_{it} ]

for the Original Model equation, and

[ g{E(Y_{it}b_i,u_i)} = X_{it} + Z_{it}b_i + S_{it}u_i + o_{it} ]

for the Final equation. Here, beta^(0) is the unpenalized reference vector and tilde(beta) is the adjusted vector from the penalized fit. The comparison table also stores adjusted - original for every coefficient.

Penalty tuning

tuned <- tune_MEMWAS(
  y ~ x1 + x2,
  data = dat,
  id = "id",
  time = "time",
  random = ~1,
  autocor = "AR(1)",
  K = 5L,
  metric = "RMSE",
  n_initial = 10L,
  max_iter = 100L,
  refit_final = TRUE,
  verbose = FALSE
)

Primary IDs are assigned as complete units to folds. A candidate remains eligible only when every fold satisfies the selected convergence policy and produces a finite validation metric. Search candidates use

[ _1 = , _2 = (1-), ]

where lambda = lambda_base^log_lambda. Selection can use grouped cross-validation, an approximation-specific full-data criterion, or an ML Laplace likelihood refit of the leading validation candidates. When refit_final = TRUE, original and adjusted coefficients and both equations are propagated to the tuning object.

Approximation and prediction tools

MEMWAS_capabilities()
diagnose_approximation(fit)

comparison <- compare_approximations(
  c("laplace", "variational_inference"),
  formula = y ~ x1 + x2,
  data = dat,
  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.