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:

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

Family-specific conditional models

Likelihood and approximation methods

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

Key parameters

The main direct interface is:

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
)

A fully explicit multi-component serial model can be declared as follows:

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:

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:

A sandwich covariance estimator is deliberately outside this version’s registered inference contract.

Within-subject correlation structures

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:

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.

Minimal educational example

The following simulation has a student-specific random intercept, an AR(1) latent serial process, and independent Gaussian observation error.

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.

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.

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:

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.