multifit: Multivariate Regression Analysis

In summata: Publication-Ready Summary Tables and Forest Plots

Multivariate Regression Analysis

Description

Performs regression analyses of a single predictor (exposure) across multiple outcomes. This function is designed for studies where a single exposure variable is tested against multiple endpoints, such as complication screening, biomarker associations, or phenome-wide association studies. Returns publication-ready formatted results with optional covariate adjustment. Supports interactions, mixed-effects models, stratification, and clustered standard errors.

Usage

multifit(
  data,
  outcomes,
  predictor,
  covariates = NULL,
  interactions = NULL,
  random = NULL,
  strata = NULL,
  cluster = NULL,
  model_type = "glm",
  family = "binomial",
  columns = "adjusted",
  p_threshold = 1,
  conf_level = 0.95,
  show_n = TRUE,
  show_events = TRUE,
  digits = 2,
  p_digits = 3,
  labels = NULL,
  predictor_label = NULL,
  include_predictor = TRUE,
  keep_models = FALSE,
  exponentiate = NULL,
  conf_method = NULL,
  parallel = TRUE,
  n_cores = NULL,
  number_format = NULL,
  verbose = NULL,
  ...
)

Arguments

data	Data frame or data.table containing the analysis dataset. The function automatically converts data frames to data.tables for efficient processing.
outcomes	Character vector of outcome variable names to analyze. Each outcome is tested in its own model with the predictor. For time-to-event analysis, use Surv() syntax for the outcome variable (e.g., c("Surv(time1, status1)", "Surv(time2, status2)")).
predictor	Character string specifying the predictor (exposure) variable name. This variable is tested against each outcome. Can be continuous or categorical (factor).
covariates	Optional character vector of covariate variable names to include in adjusted models. When specified, models are fit as outcome ~ predictor + covariate1 + covariate2 + ..., and only the predictor effect is reported. Default is NULL (unadjusted models).
interactions	Optional character vector of interaction terms to include in adjusted models, using colon notation (e.g., c("predictor:sex")). Interactions involving the predictor will have their effects extracted and reported. Default is NULL.
random	Optional character string specifying random effects formula for mixed effects models (e.g., "(1|hospital)" or "(1|site/patient)"). Required when model_type is "glmer", "lmer", or "coxme" unless random effects are included in the covariates vector. Alternatively, random effects can be included directly in the covariates vector using the same syntax (e.g., covariates = c("age", "sex", "(1|site)")). Default is NULL.
strata	Optional character string naming the stratification variable for Cox or conditional logistic models. Creates separate baseline hazards for each stratum. Default is NULL.
cluster	Optional character string naming the clustering variable for Cox models. Computes robust clustered standard errors. Default is NULL.
model_type	Character string specifying the type of regression model to fit. Options include: "glm" - Generalized linear model (default). Supports multiple distributions via the family parameter including logistic, Poisson, Gamma, Gaussian, and quasi-likelihood models. "lm" - Linear regression for continuous outcomes with normally distributed errors. "coxph" - Cox proportional hazards model for time-to-event survival analysis. Requires Surv() outcome syntax. "clogit" - Conditional logistic regression for matched case-control studies. "negbin" - Negative binomial regression for overdispersed count data (requires MASS package). Estimates an additional dispersion parameter compared to Poisson regression. "glmer" - Generalized linear mixed-effects model for hierarchical or clustered data with non-normal outcomes (requires lme4 package and random parameter). "lmer" - Linear mixed-effects model for hierarchical or clustered data with continuous outcomes (requires lme4 package and random parameter). "coxme" - Cox mixed-effects model for clustered survival data (requires coxme package and random parameter).
family	For GLM and GLMER models, specifies the error distribution and link function. Can be a character string, a family function, or a family object. Ignored for non-GLM/GLMER models. Binary/Binomial outcomes: "binomial" or binomial() - Logistic regression for binary outcomes (0/1, TRUE/FALSE). Returns odds ratios (OR). Default. "quasibinomial" or quasibinomial() - Logistic regression with overdispersion. Use when residual deviance >> residual df. binomial(link = "probit") - Probit regression (normal CDF link). binomial(link = "cloglog") - Complementary log-log link for asymmetric binary outcomes. Count outcomes: "poisson" or poisson() - Poisson regression for count data. Returns rate ratios (RR). Assumes mean = variance. "quasipoisson" or quasipoisson() - Poisson regression with overdispersion. Use when variance > mean. Continuous outcomes: "gaussian" or gaussian() - Normal/Gaussian distribution for continuous outcomes. Equivalent to linear regression. gaussian(link = "log") - Log-linear model for positive continuous outcomes. Returns multiplicative effects. Positive continuous outcomes: "Gamma" or Gamma() - Gamma distribution for positive, right-skewed continuous data (e.g., costs, lengths of stay). When passed as a string, resolves to log link for interpretable multiplicative effects. Gamma(link = "inverse") - Gamma with inverse (canonical) link. Gamma(link = "identity") - Gamma with identity link for additive effects on positive outcomes. "inverse.gaussian" or inverse.gaussian() - Inverse Gaussian for positive, highly right-skewed data. For negative binomial regression (overdispersed counts), use model_type = "negbin" instead of the family parameter. See family for additional details and options.
columns	Character string specifying which result columns to display when both unadjusted and adjusted models are fit (i.e., when covariates is specified): "adjusted" - Show only adjusted (covariate-controlled) results [default] "unadjusted" - Show only unadjusted (crude) results "both" - Show both unadjusted and adjusted results side-by-side Ignored when covariates = NULL.
p_threshold	Numeric value between 0 and 1 specifying a p-value threshold for filtering results. Only outcomes with p-value less than or equal to the threshold are included in the output. Default is 1 (no filtering; all outcomes returned).
conf_level	Numeric confidence level for confidence intervals. Must be between 0 and 1. Default is 0.95 (95% confidence intervals).
show_n	Logical. If TRUE, includes the sample size column in the output table. Default is TRUE.
show_events	Logical. If TRUE, includes the events column in the output table (relevant for survival and logistic regression). Default is TRUE.
digits	Integer specifying the number of decimal places for effect estimates (OR, HR, RR, coefficients). Default is 2.
p_digits	Integer specifying the number of decimal places for p-values. Values smaller than 10^(-p_digits) are displayed as "< 0.001" (for p_digits = 3), "< 0.0001" (for p_digits = 4), etc. Default is 3.
labels	Named character vector or list providing custom display labels for variables. Can include labels for outcomes, predictors, and covariates. Names should match variable names, values are the display labels. Labels are applied to: (1) outcome names in the Outcome column, (2) predictor variable name when displayed, and (3) variable names in formatted interaction terms. Variables not in labels use their original names. Default is NULL.
predictor_label	Optional character string providing a custom display label for the predictor variable. Takes precedence over labels for the predictor. Default is NULL (uses label from labels or original name).
include_predictor	Logical. If TRUE (default), includes the predictor column in the output table showing which level of a factor predictor is being compared. If FALSE, omits the predictor column, which may be useful when the predictor information will be explained in a table caption or figure legend.
keep_models	Logical. If TRUE, stores all fitted model objects in the output as an attribute. Models are accessible via attr(result, "models"). Default is FALSE.
exponentiate	Logical. Whether to exponentiate coefficients (display OR/HR/RR instead of log odds/log hazards). Default is NULL, which automatically displays raw coefficients for linear models and exponentiates for logistic, Poisson, and Cox models.
conf_method	Character string controlling the confidence interval method. If NULL (default), uses getOption("summata.conf_method", "profile"). "profile" - Profile likelihood intervals for GLM and negative binomial models (via MASS::confint.glm()), exact t-distribution intervals for linear models. Falls back to Wald on profiling failure. Quasi-likelihood families always use Wald (no true likelihood). "wald" - Wald intervals (coefficient \pm z \times SE) for all model types. Faster but less accurate near boundary conditions or with small subgroups. Cox and mixed-effects models use Wald intervals regardless of this setting. Set globally with options(summata.conf_method = "wald") to use Wald throughout a session.
parallel	Logical. If TRUE (default), fits models in parallel for improved performance with many outcomes.
n_cores	Integer specifying the number of CPU cores to use for parallel processing. Default is NULL (auto-detect: uses detectCores() - 1). Ignored when parallel = FALSE.
number_format	Character string or two-element character vector controlling thousand and decimal separators in formatted output. Named presets: "us" - Comma thousands, period decimal: 1,234.56 [default] "eu" - Period thousands, comma decimal: 1.234,56 "space" - Thin-space thousands, period decimal: 1 234.56 (SI/ISO 31-0) "none" - No thousands separator: 1234.56 Or provide a custom two-element vector c(big.mark, decimal.mark), e.g., c("'", ".") for Swiss-style: ⁠1'234.56⁠. When NULL (default), uses getOption("summata.number_format", "us"). Set the global option once per session to avoid passing this argument repeatedly: options(summata.number_format = "eu")
verbose	Logical. If TRUE, displays model fitting warnings (e.g., singular fit, convergence issues). If FALSE (default), routine fitting messages are suppressed while unexpected warnings are preserved. When NULL, uses getOption("summata.verbose", FALSE).
...	Additional arguments passed to the underlying model fitting functions.

Details

Analysis Approach:

The function implements a multivariate (multi-outcome) screening workflow that inverts the typical regression paradigm:

1. For each outcome in outcomes, fits a separate model with the predictor as the main exposure

2. If covariates specified, fits adjusted model: outcome ~ predictor + covariates + interactions

3. Extracts only the predictor effect(s) from each model, ignoring covariate coefficients

4. Combines results into a single table for comparison across outcomes

5. Optionally filters by p-value threshold

This is conceptually opposite to uniscreen(), which tests multiple predictors against a single outcome. Use multifit() when you have one exposure of interest and want to screen across multiple endpoints.

When to Use Multivariate Regression Analysis:

• Complication screening: Test one exposure (e.g., operative time, BMI, biomarker level) against multiple postoperative complications

• Treatment effects: Test one treatment against multiple efficacy and safety endpoints simultaneously

• Biomarker studies: Test one biomarker against multiple clinical outcomes to understand its prognostic value

• Phenome-wide association studies (PheWAS): Test genetic variants or exposures against many phenotypes

• Risk factor profiling: Understand how one risk factor relates to a spectrum of outcomes

Handling Categorical Predictors:

When the predictor is a factor variable with multiple levels:

• Each non-reference level gets its own row for each outcome

• Reference category is determined by factor level ordering

• The Predictor column shows "Variable (Level)" format (e.g., "Treatment (Drug A)", "Treatment (Drug B)")

• For binary variables with affirmative non-reference levels (Yes, 1, True, Present, Positive, +), shows just "Variable" (e.g., "Diabetes" instead of "Diabetes (Yes)")

• Effect estimates compare each level to the reference

Adjusted vs. Unadjusted Results:

When covariates is specified, the function fits both models but only extracts predictor effects:

• columns = "adjusted": Reports only covariate-adjusted effects. Column labeled "aOR/aHR," etc.

• columns = "unadjusted": Reports only crude effects. Column labeled "OR/HR," etc.

• columns = "both": Reports both side-by-side. Useful for identifying confounding (large change in effect) or independent effects (similar estimates)

Interaction Terms:

When interactions includes terms involving the predictor:

• Main effect of predictor is always reported

• Interaction effects are extracted and displayed with formatted names

• Format: Variable (Level) × Variable (Level) using multiplication sign notation

• Useful for testing effect modification (e.g., does treatment effect differ by sex?)

Mixed-Effects Models:

For clustered or hierarchical data (e.g., patients within hospitals):

• Use model_type = "glmer" with random = "(1|cluster)" for random intercept models

• Nested random effects: random = "(1|site/patient)"

• Crossed random effects: random = "(1|site) + (1|doctor)"

• For survival outcomes, use model_type = "coxme"

Stratification and Clustering (Cox models):

For Cox proportional hazards models:

• strata: Creates separate baseline hazards for each stratum level. Use when hazards are non-proportional across strata but stratum effects do not need to be estimated

• cluster: Computes robust (sandwich) standard errors accounting for within-cluster correlation. Alternative to mixed effects when only robust SEs are needed

Filtering based on p-value:

The p_threshold parameter filters results after fitting all models:

• Only outcomes with p less than or equal to the threshold are retained in output

• For factor predictors, outcome is kept if any level is significant

• Useful for focusing on significant associations in exploratory analyses

• Default is 1 (no filtering) - recommended for confirmatory analyses

Outcome Homogeneity:

All outcomes in a single multifit() call should be of the same type (all binary, all continuous, or all survival). Mixing outcome types produces tables with incompatible effect measures (e.g., odds ratios alongside regression coefficients), which can mislead readers. The function validates outcome compatibility and issues a warning when mixed types are detected.

For analyses involving multiple outcome types, run separate multifit() calls for each type:

# Binary outcomes
binary_results <- multifit(data, outcomes = c("death", "readmission"),
                           predictor = "treatment", model_type = "glm")

# Continuous outcomes
continuous_results <- multifit(data, outcomes = c("los_days", "cost"),
                               predictor = "treatment", model_type = "lm")

Effect Measures by Model Type:

• Logistic (model_type = "glm", family = "binomial"): Odds ratios (OR/aOR)

• Cox (model_type = "coxph"): Hazard ratios (HR/aHR)

• Poisson (model_type = "glm", family = "poisson"): Rate ratios (RR/aRR)

• Linear (model_type = "lm"): Coefficient estimates

• Mixed effects: Same as fixed-effects counterparts

Memory and Performance:

• parallel = TRUE (default) uses multiple cores for faster fitting

• keep_models = FALSE (default) discards model objects to save memory

• For many outcomes, parallel processing provides substantial speedup

• Set keep_models = TRUE only when you need model diagnostics

Value

A data.table with S3 class "multifit_result" containing formatted multivariate regression results. The table structure includes:

Outcome

Character. Outcome variable name or custom label

Predictor

Character. For factor predictors: formatted as "Variable (Level)" showing the level being compared to reference. For binary variables where the non-reference level is an affirmative value (Yes, 1, True, Present, Positive, +), shows just "Variable". For continuous predictors: the variable name. For interactions: the formatted interaction term (e.g., "Treatment (Drug A) × Sex (Male)")

n

Integer. Sample size used in the model (if show_n = TRUE)

Events

Integer. Number of events (if show_events = TRUE)

OR/HR/RR/Coefficient (95% CI)

Character. Unadjusted effect estimate with CI (if columns = "unadjusted" or "both")

aOR/aHR/aRR/Adj. Coefficient (95% CI)

Character. Adjusted effect estimate with CI (if columns = "adjusted" or "both")

Uni p / Multi p / p-value

Character. Formatted p-value(s). Column names depend on columns setting

The returned object includes the following attributes accessible via attr():

raw_data

data.table. Unformatted numeric results with separate columns for effect estimates, standard errors, confidence intervals, and p-values. Suitable for custom analysis or visualization

models

list (if keep_models = TRUE). Named list of fitted model objects, with outcome names as list names. Each element contains $unadjusted and/or $adjusted models depending on settings

predictor

Character. The predictor variable name

outcomes

Character vector. The outcome variable names

covariates

Character vector or NULL. The covariate variable names

interactions

Character vector or NULL. The interaction terms

random

Character or NULL. The random effects formula

strata

Character or NULL. The stratification variable

cluster

Character or NULL. The clustering variable

model_type

Character. The regression model type used

columns

Character. Which columns were displayed

analysis_type

Character. "multi_outcome" to identify analysis type

significant

Character vector. Names of outcomes with p < 0.05 for the predictor (uses adjusted p-values when available)

See Also

uniscreen for screening multiple predictors against one outcome, multiforest for creating forest plots from multifit results, fit for single-outcome regression with full coefficient output, fullfit for complete univariable-to-multivariable workflow

Other regression functions: compfit(), fit(), fullfit(), print.compfit_result(), print.fit_result(), print.fullfit_result(), print.multifit_result(), print.uniscreen_result(), uniscreen()

Examples

# Load example data
data(clintrial)
data(clintrial_labels)

# Example 1: Basic multivariate analysis (unadjusted)
# Test treatment effect on multiple binary outcomes
result1 <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    labels = clintrial_labels,
    parallel = FALSE
)
print(result1)
# Shows odds ratios comparing Drug A and Drug B to Control



# Example 2: Adjusted analysis with covariates
# Adjust for age, sex, and disease stage
result2 <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    labels = clintrial_labels,
    parallel = FALSE
)
print(result2)
# Shows adjusted odds ratios (aOR)

# Example 3: Compare unadjusted and adjusted results
result3 <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    columns = "both",
    labels = clintrial_labels,
    parallel = FALSE
)
print(result3)
# Useful for identifying confounding effects

# Example 4: Continuous predictor across outcomes
# Test age effect on multiple outcomes
result4 <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "age",
    covariates = c("sex", "treatment", "stage"),
    labels = clintrial_labels,
    parallel = FALSE
)
print(result4)
# One row per outcome for continuous predictor

# Example 5: Cox regression for survival outcomes
library(survival)
cox_result <- multifit(
    data = clintrial,
    outcomes = c("Surv(pfs_months, pfs_status)", 
                 "Surv(os_months, os_status)"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    model_type = "coxph",
    labels = clintrial_labels,
    parallel = FALSE
)
print(cox_result)
# Returns hazard ratios (HR/aHR)

# Example 6: Cox with stratification by site
cox_strat <- multifit(
    data = clintrial,
    outcomes = c("Surv(os_months, os_status)"),
    predictor = "treatment",
    covariates = c("age", "sex"),
    strata = "site",
    model_type = "coxph",
    labels = clintrial_labels,
    parallel = FALSE
)
print(cox_strat)

# Example 7: Cox with clustered standard errors
cox_cluster <- multifit(
    data = clintrial,
    outcomes = c("Surv(os_months, os_status)"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    cluster = "site",
    model_type = "coxph",
    labels = clintrial_labels,
    parallel = FALSE
)
print(cox_cluster)

# Example 8: Interaction between predictor and covariate
# Test if treatment effect differs by sex
result_int <- multifit(
    data = clintrial,
    outcomes = c("surgery", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    interactions = c("treatment:sex"),
    labels = clintrial_labels,
    parallel = FALSE
)
print(result_int)
# Shows main effects and interaction terms with × notation

# Example 9: Linear model for continuous outcomes
linear_result <- multifit(
    data = clintrial,
    outcomes = c("los_days", "biomarker_x"),
    predictor = "treatment",
    covariates = c("age", "sex"),
    model_type = "lm",
    labels = clintrial_labels,
    parallel = FALSE
)
print(linear_result)
# Returns coefficient estimates, not ratios

# Example 10: Poisson regression for equidispersed count outcomes
# fu_count has variance ~= mean, appropriate for standard Poisson
poisson_result <- multifit(
    data = clintrial,
    outcomes = c("fu_count"),
    predictor = "treatment",
    covariates = c("age", "stage"),
    model_type = "glm",
    family = "poisson",
    labels = clintrial_labels,
    parallel = FALSE
)
print(poisson_result)
# Returns rate ratios (RR)
# For overdispersed counts (ae_count), use model_type = "negbin" instead

# Example 11: Filter to significant results only
sig_results <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "stage",
    p_threshold = 0.05,
    labels = clintrial_labels,
    parallel = FALSE
)
print(sig_results)
# Only outcomes with significant associations shown

# Example 12: Custom outcome labels
result_labeled <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    labels = c(
        surgery = "Surgical Resection",
        pfs_status = "Disease Progression",
        os_status = "Death",
        treatment = "Treatment Group"
    ),
    parallel = FALSE
)
print(result_labeled)

# Example 13: Keep models for diagnostics
result_models <- multifit(
    data = clintrial,
    outcomes = c("surgery", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex"),
    keep_models = TRUE,
    parallel = FALSE
)

# Access stored models
models <- attr(result_models, "models")
names(models)

# Get adjusted model for surgery outcome
surgery_model <- models$surgery$adjusted
summary(surgery_model)

# Example 14: Access raw numeric data
result <- multifit(
    data = clintrial,
    outcomes = c("surgery", "os_status"),
    predictor = "age",
    parallel = FALSE
)

# Get unformatted results for custom analysis
raw_data <- attr(result, "raw_data")
print(raw_data)
# Contains exp_coef, ci_lower, ci_upper, p_value, \emph{etc.}

# Example 15: Hide sample size and event columns
result_minimal <- multifit(
    data = clintrial,
    outcomes = c("surgery", "os_status"),
    predictor = "treatment",
    show_n = FALSE,
    show_events = FALSE,
    parallel = FALSE
)
print(result_minimal)

# Example 16: Customize decimal places
result_digits <- multifit(
    data = clintrial,
    outcomes = c("surgery", "os_status"),
    predictor = "age",
    digits = 3,
    p_digits = 4,
    parallel = FALSE
)
print(result_digits)

# Example 17: Force coefficient display (no exponentiation)
result_coef <- multifit(
    data = clintrial,
    outcomes = c("surgery"),
    predictor = "age",
    exponentiate = FALSE,
    parallel = FALSE
)
print(result_coef)

# Example 18: Complete publication workflow
final_table <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage", "grade"),
    columns = "both",
    labels = clintrial_labels,
    digits = 2,
    p_digits = 3,
    parallel = FALSE
)
print(final_table)

# Example 19: Gamma regression for positive continuous outcomes
gamma_result <- multifit(
    data = clintrial,
    outcomes = c("los_days", "recovery_days"),
    predictor = "treatment",
    covariates = c("age", "surgery"),
    model_type = "glm",
    family = Gamma(link = "log"),
    labels = clintrial_labels,
    parallel = FALSE
)
print(gamma_result)
# Returns multiplicative effects on positive continuous data

# Example 20: Quasipoisson for overdispersed counts
quasi_result <- multifit(
    data = clintrial,
    outcomes = c("ae_count"),
    predictor = "treatment",
    covariates = c("age", "diabetes"),
    model_type = "glm",
    family = "quasipoisson",
    labels = clintrial_labels,
    parallel = FALSE
)
print(quasi_result)
# Adjusts standard errors for overdispersion

# Example 21: Generalized linear mixed effects (GLMER)
# Test treatment across outcomes with site clustering
if (requireNamespace("lme4", quietly = TRUE)) {
    glmer_result <- suppressWarnings(multifit(
        data = clintrial,
        outcomes = c("surgery", "pfs_status", "os_status"),
        predictor = "treatment",
        covariates = c("age", "sex"),
        random = "(1|site)",
        model_type = "glmer",
        family = "binomial",
        labels = clintrial_labels,
        parallel = FALSE
    ))
    print(glmer_result)
}

# Example 22: Cox mixed effects with random site effects
if (requireNamespace("coxme", quietly = TRUE)) {
    coxme_result <- multifit(
        data = clintrial,
        outcomes = c("Surv(pfs_months, pfs_status)", 
                     "Surv(os_months, os_status)"),
        predictor = "treatment",
        covariates = c("age", "sex", "stage"),
        random = "(1|site)",
        model_type = "coxme",
        labels = clintrial_labels,
        parallel = FALSE
    )
    print(coxme_result)
}

# Example 23: Multiple interactions across outcomes
multi_int <- multifit(
    data = clintrial,
    outcomes = c("surgery", "pfs_status", "os_status"),
    predictor = "treatment",
    covariates = c("age", "sex", "stage"),
    interactions = c("treatment:stage", "treatment:sex"),
    labels = clintrial_labels,
    parallel = FALSE
)
print(multi_int)
# Shows how treatment effects vary by stage and sex across outcomes

summata documentation built on May 7, 2026, 5:07 p.m.