Fit Linear or Mixed Effects Models Across Multiple Variable Combinations

This function fits linear (`lm`) or mixed effects (`lmer`) models across combinations of outcomes, exposures, effect modifiers, and sensitivity covariates. It supports interaction terms, random effects, and returns a tidy summary of each model. It is designed to simplify model fitting across many variable combinations with consistent covariates and random effects.

For internal testing purposes, a hidden `.test_force` argument can be passed via `...` to simulate interactive, top-level calls. This enables reliable testing of preview and reminder output without relying on `interactive()` or `sys.nframe()`.

Arguments

data: A data frame containing the variables used in the model.
outcome: A character vector of outcome variable names.
exposure: A character vector of exposure variable names.
covariates: A character vector of covariate names to adjust for (e.g., confounders). These are included in all models.
effect_modifier: A character vector of effect modifiers to interact with the exposure. All combinations will be tested.
sensitivity_cov: A character vector of sensitivity covariates to include in the model. All combinations will be tested.
random_effects: A string specifying the random effects structure (e.g., "(1 | group)"). If `NULL`, a linear model is used.
p_values: Logical. If `TRUE` (default), uses `lmerTest` to compute p-values for mixed models. If `FALSE`, uses `lme4` without p-values.
verbose: Logical. If `TRUE` (default), prints progress and model status messages.
return_grid: Logical. If `TRUE`, attaches the model grid as an attribute to the returned object.

Value

A named list of model results. Each element contains:

model: The fitted model object.
tidy: A tidy summary of the model (term, estimate, confidence intervals, p-values, etc.).
residuals: Model residuals.
formula: The model formula used.

If multiple combinations are tested, the list is named using the format `outcome&exposure&modifier&sensitivity`.

A list containing model results. If multiple outcomes or exposures are provided, returns a named list of results.

Note

For internal testing purposes, a `.test_force` flag can be passed via the `.internal` argument (as a named list) to simulate interactive, top-level calls. This is used only in tests and should not be used in production.

Examples

# Examples

data <- mtcars

if (FALSE) { # \dontrun{
# 1. Run a single model with one outcome and one exposure
run_linear_models(data = data,
                  outcome = "mpg",
                  exposure = "cyl")
} # }

# 2. Store the output for a single model
model_output <- run_linear_models(data = data,
                                  outcome = "mpg",
                                  exposure = "cyl")
#> Running model with outcome = mpg, exposure = cyl
model_output$tidy  # access tidy results
#> NULL

# 3. Store output for multiple outcomes and exposures
multi_output <- run_linear_models(data = data,
                                  outcome = c("mpg", "hp"),
                                  exposure = c("cyl", "wt"))
#> Running model with outcome = mpg, exposure = cyl
#> Running model with outcome = mpg, exposure = wt
#> Running model with outcome = hp, exposure = cyl
#> Running model with outcome = hp, exposure = wt
multi_output[["mpg&cyl"]]$tidy  # access a specific model's tidy results
#> # A tibble: 2 × 11
#>   term  estimate conf.low conf.high std.error  p.value error n_obs   BIC formula
#>   <chr>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl> <lgl> <int> <dbl> <chr>  
#> 1 (Int…    37.9     33.6      42.1      2.07  8.37e-18 NA       32  174. mpg ~ …
#> 2 cyl      -2.88    -3.53     -2.22     0.322 6.11e-10 NA       32  174. mpg ~ …
#> # ℹ 1 more variable: note <chr>

# 4. Use named vectors for outcomes and exposures
outcomes <- c("mpg", "hp")
exposures <- c("cyl", "wt")
named_output <- run_linear_models(data = data,
                                  outcome = outcomes,
                                  exposure = exposures)
#> Running model with outcome = mpg, exposure = cyl
#> Running model with outcome = mpg, exposure = wt
#> Running model with outcome = hp, exposure = cyl
#> Running model with outcome = hp, exposure = wt
named_output[["hp&wt"]]$tidy
#> # A tibble: 2 × 11
#>   term   estimate conf.low conf.high std.error p.value error n_obs   BIC formula
#>   <chr>     <dbl>    <dbl>     <dbl>     <dbl>   <dbl> <lgl> <int> <dbl> <chr>  
#> 1 (Inte…    -1.82    -67.8      64.2     32.3  9.55e-1 NA       32  353. hp ~ wt
#> 2 wt        46.2      26.5      65.8      9.63 4.15e-5 NA       32  353. hp ~ wt
#> # ℹ 1 more variable: note <chr>

# 5. Include covariates directly
covar_output <- run_linear_models(data = data,
                                  outcome = "mpg",
                                  exposure = "cyl",
                                  covariates = c("wt", "hp"))
#> Running model with outcome = mpg, exposure = cyl
covar_output$tidy
#> NULL

# 6. Include covariates using a named vector
covars <- c("wt", "hp")
names(covars) <- covars
named_covar_output <- run_linear_models(data = data,
                                        outcome = "mpg",
                                        exposure = "cyl",
                                        covariates = covars)
#> Running model with outcome = mpg, exposure = cyl
named_covar_output$tidy
#> NULL

# 7. Include an effect modifier (e.g., interaction term)
interaction_output <- run_linear_models(data = data,
                                        outcome = "mpg",
                                        exposure = "cyl",
                                        effect_modifier = "am")
#> Running model with outcome = mpg, exposure = cyl, effect_modifier = am
interaction_output$tidy
#> NULL

# 8. Model with a non-existent variable (graceful failure)
error_output <- run_linear_models(data = data,
                                  outcome = "mpg",
                                  exposure = "not_a_column")
#> Running model with outcome = mpg, exposure = not_a_column
#> Model failed for formula: mpg ~ not_a_column
#> Error: object 'not_a_column' not found

error_output$tidy
#> NULL