# MEMWAS

 <table style="border: none; border-collapse: collapse; margin-left: 0; margin-right: auto;">
   <tr style="border: none;">
     <td style="border: none; vertical-align: middle;">
      <img src = "https://github.com/EnochKang/MEMWAS/blob/main/vignettes/MEMWAS_Logo.png?raw=true" align = "left" width = "120" style="margin-left: 20px;" />
     </td>
     <td style="border: none; vertical-align: middle; padding-left: 20px; text-align: left; line-height: 1.5;">
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;"> M</span>ixed<br>
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;">E</span>ffects<br>
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;">M</span>odels<br>
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;">W</span>ith<br>
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;">A</span>utocorrelation<br>
       <span style="color: #00008B; font-size: 1.2em; font-weight: bold;">S</span>tructures
     </td>
   </tr>
      <tr style="border: none;">
     <td style="border: none; vertical-align: middle; padding-left: 20px; text-align: center;">
       Version 0.9.3 <br>
       Aug. 04, 2026
     </td>
     <td style="border: none; vertical-align: middle; padding-left: 20px; text-align: left; line-height: 1.5;">
     </td>
   </tr>
 </table>

&nbsp;

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

```r
# 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

```r
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:

```r
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}(\theta_0) +
      \sum_j \operatorname{diag}(x_{ij}) K_{ji}(\theta_j)
      \operatorname{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.

## Ranking temporal autocorrelation structures

```r
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:

```r
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

```r
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

```r
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.

```r
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)} + 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} + 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

```r
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

\[
\lambda_1 = \alpha\lambda,\qquad
\lambda_2 = (1-\alpha)\lambda,
\]

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

```r
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.


