Estimate Effect Sizes (B) for Variants Based on MAF

Source: R/get_B.R

Trains the best regression model from training data and predicts allelic effect sizes (B) for target MAF values in one call. This is a convenience wrapper around train_B_model (training) and predict_B (prediction).

Two estimation methods are supported: (1) direct beta method using effect sizes, and (2) p-value/Z-score method using trait-independent h-squared. When both are applicable, the function compares them and selects the better-fitting model as the primary result.

For workflows that train once and predict many times (e.g., per-chromosome analysis), use train_B_model and predict_B directly instead.

get_B(
  training_trait,
  training_MAF,
  training_BETA = NULL,
  training_P = NULL,
  training_P_mlog10 = NULL,
  training_N = NULL,
  target_trait = NULL,
  target_MAF,
  target_SE = NULL,
  target_case_prop = NULL,
  method = "auto",
  selection_criterion = "R2",
  custom_models = NULL,
  outlier_method = "none",
  outlier_action = "flag",
  cook_threshold = 4,
  biological_rules = NULL,
  show_model_selection = TRUE,
  return_full = FALSE,
  verbose = 1
)

Arguments

training_trait

Character string or NULL, trait type of training data. Can be "binary", "continuous", "mixed" (when training data contains both trait types), or NULL (unspecified). When "mixed" or NULL, only the p-value/Z-score method is available since the direct beta method requires a single trait type

training_MAF

Numeric vector of minor allele frequencies from training data, used to fit the B prediction model. Values should be in (0, 0.5]. Values > 0.5 are automatically folded to 1 - MAF with a warning

training_BETA

Numeric vector of effect sizes (regression coefficients) from training data. Required for direct beta method. Should have the same length as training_MAF

training_P

Numeric vector of p-values from training data. Used for p-value/Z-score method. Values should be in (0, 1). Either training_P OR training_P_mlog10 should be provided (not both) for p-value method

training_P_mlog10

Numeric vector of -log10 transformed p-values from training data. Use this instead of training_P when dealing with extremely small p-values (e.g., from GWAS) to avoid numerical underflow. For example, if p = 1e-300, use training_P_mlog10 = 300. Either training_P OR training_P_mlog10 should be provided (not both) for p-value method

training_N

Numeric vector of sample sizes from training data. Required for p-value/Z-score method

target_trait

Character string or NULL, trait type of target data. Must be either "binary", "continuous", or NULL (default). For the beta method, target_trait can be omitted if you are confident that training and target data share the same trait type. For the p-value method, target_trait is required. When NULL, appropriate warnings will be issued to ensure scientific validity

target_MAF

Numeric vector of minor allele frequencies from target data for which to predict effect sizes. Values should be in (0, 0.5]. Values > 0.5 are automatically folded to 1 - MAF with a warning

target_SE

Numeric vector of standard errors from target data. Required when target_trait == "continuous" for p-value/Z-score method

target_case_prop

Numeric value or vector specifying the proportion of cases in the target case-control study. This is calculated as: number_of_cases / (number_of_cases + number_of_controls). For example, in a balanced study with 500 cases and 500 controls, target_case_prop = 500/1000 = 0.5. IMPORTANT: This is NOT the population disease prevalence. For a rare disease with 0.1\ have target_case_prop = 0.5, not 0.001. Required when target_trait == "binary" for p-value/Z-score method. Values must be in (0, 1)

method

Character string specifying which method to use: "auto" (default), "beta", "pvalue", or "both". "auto" uses both methods when same trait and beta available, otherwise uses pvalue method. "both" runs both methods for comparison when both are applicable

selection_criterion

Character string specifying model selection criterion: "R2" (default), "adj_R2", or "CV_R2". See select_best_model for detailed guidance. Note: "CV" is deprecated, use "CV_R2" instead

custom_models

Optional list of custom model formulas to include in candidate model set. Each element should be a formula object (e.g., formula(Y ~ poly(X, 2))). Custom models are fit alongside the 8 standard models

outlier_method

Character string specifying outlier detection method: "none" (default, no detection), "statistical" (Cook's distance), "biological" (implausible MAF-effect combinations), or "both"

outlier_action

Character string specifying action for detected outliers: "flag" (default, report but don't remove) or "remove" (exclude from analysis)

cook_threshold

Numeric threshold multiplier for Cook's distance outlier detection (default: 4). Outliers are defined as points with Cook's D > cook_threshold/n

biological_rules

Optional list of custom rules for biological outlier detection. Default uses common variants (MAF > 0.05) with large effects (|effect| > 10)

show_model_selection

Logical. If TRUE (default), prints the model selection comparison table when estimating B. Set to FALSE to suppress output in loops or simulations

return_full

Logical. If FALSE (default), returns only the primary B estimates as a numeric vector. If TRUE, returns a comprehensive result object of class glow_B_estimate

verbose

Integer controlling verbosity: 0=silent, 1=warnings (default), 2=info messages, 3=debug output

Value

If return_full = FALSE (default): Numeric vector of predicted effect sizes (B) for each variant in target_MAF.

If return_full = TRUE: An object of class glow_B_estimate containing:

B

Primary B predictions (numeric vector)

B_beta_method

B predictions from beta method (NULL if not used)

B_pvalue_method

B predictions from p-value method (NULL if not used)

model

A glow_B_model object (from train_B_model) containing all training results: fitted models, outlier info, training data, and training-time comparison

target_summary

List with n_predictions, trait_type, MAF_range

When both methods were used, the model$comparison field is augmented with prediction-based comparison statistics (correlation, RMSE, percent difference).

Details

This function is a convenience wrapper equivalent to:


  b_model <- train_B_model(...)
  B <- predict_B(b_model, target_MAF, ...)

For "train once, predict many times" workflows, call train_B_model and predict_B separately.

Two Estimation Methods:

Method 1: Direct Beta Method

  • Uses training_BETA values directly

  • Fits regression model: \(BETA^2 \sim f(MAF)\)

  • Best when training and target have the same trait type

  • Provides most direct relationship between MAF and effect size

Method 2: P-value/Z-score Method

  • Converts p-values to Z-squared using inverse chi-squared CDF

  • Calculates trait-independent function: \(h^2(q) = Z^2 / (2 * N * q * (1-q))\)

  • Fits regression model: \(h^2 \sim f(MAF)\)

  • Converts back to trait-specific B:

    • Continuous: \(B = h(q) * SE\)

    • Binary: \(B = h(q) / \sqrt{p_0 * (1-p_0)}\) where \(p_0\) is the case proportion in the case-control study

  • Can transfer information between different trait types

Method Selection Logic:

  • method="auto" (default): Uses "both" if same trait + BETA available, else "pvalue"

  • method="beta": Only direct beta method (requires training_BETA)

  • method="pvalue": Only p-value/Z-score method (requires training_P/training_N)

  • method="both": Runs both methods for comparison

Primary Method Selection (when method="both"):

When both methods are run, the function selects the better method as primary based on model quality. The method with the higher selection_criterion value (R2/adj_R2/CV_R2) is chosen. The comparison is done on the training data models (BETA^2 vs h^2).

Outlier Detection:

When outlier_method is not "none", outliers are detected using statistical (Cook's distance), biological (implausible combinations), or both methods. The outlier_action parameter controls whether outliers are flagged or removed.

Computational Complexity

O(8 * n_train + n_test) where:

  • n_train = length(training_MAF) (fitting 8 models)

  • n_test = length(target_MAF) (prediction step)

References

Zhang, H., Liu, M., Landers, J. E., and Wu, Z. (2026). Integrated weighted association test with application to genetic association studies. The Annals of Applied Statistics, 20(3), 2074-2097. doi:10.1214/26-AOAS2223

See also

train_B_model for training models separately predict_B for applying trained models to new data select_best_model for the model selection algorithm Optimal_Weights_M for calculating optimal weights using B

Examples

if (FALSE) { # \dontrun{
# Example 1: Direct beta method (same trait type)
set.seed(123)
training_MAF <- runif(50, 0.001, 0.3)
training_BETA <- sqrt(0.5 * training_MAF * (1 - training_MAF)) * 0.1
target_MAF <- runif(10, 0.001, 0.5)

B_estimates <- get_B(
  training_trait = "binary",
  training_MAF = training_MAF,
  training_BETA = training_BETA,
  target_trait = "binary",
  target_MAF = target_MAF,
  method = "beta"
)

# Example 2: P-value/Z-score method
training_P <- runif(50, 0.0001, 0.1)
training_N <- rep(1000, 50)

B_pval <- get_B(
  training_trait = "binary",
  training_MAF = training_MAF,
  training_P = training_P,
  training_N = training_N,
  target_trait = "binary",
  target_MAF = target_MAF,
  target_case_prop = 0.1,
  method = "pvalue"
)

# Example 3: Full result with both methods
result_full <- get_B(
  training_trait = "binary",
  training_MAF = training_MAF,
  training_BETA = training_BETA,
  training_P = training_P,
  training_N = training_N,
  target_trait = "binary",
  target_MAF = target_MAF,
  target_case_prop = 0.1,
  method = "both",
  return_full = TRUE,
  verbose = 2
)
print(result_full)
print(result_full$model)  # nested glow_B_model

# Example 4: Transfer across trait types (continuous to binary)
training_P_cont <- runif(50, 0.0001, 0.1)
training_N_cont <- rep(2000, 50)

B_cross_trait <- get_B(
  training_trait = "continuous",
  training_MAF = training_MAF,
  training_P = training_P_cont,
  training_N = training_N_cont,
  target_trait = "binary",
  target_MAF = target_MAF,
  target_case_prop = 0.5,
  method = "pvalue"  # Only pvalue method works for different traits
)

# Example 5: Suppress model selection output (useful in loops)
for (i in 1:5) {
  B <- get_B(
    training_trait = "binary",
    training_MAF = training_MAF,
    training_BETA = training_BETA,
    target_trait = "binary",
    target_MAF = target_MAF,
    show_model_selection = FALSE  # Suppress table in loop
  )
}

# Example 6: Extremely small p-values using -log10 transformation
# This is useful for GWAS data with very significant hits
training_P_mlog10 <- runif(50, 10, 350)  # -log10(p) from 10 to 350
# Equivalent to p-values from 1e-10 to 1e-350
training_N_gwas <- rep(50000, 50)  # Large GWAS sample size

B_extreme_pval <- get_B(
  training_trait = "binary",
  training_MAF = training_MAF,
  training_P_mlog10 = training_P_mlog10,  # Use -log10(p) for numerical stability
  training_N = training_N_gwas,
  target_trait = "binary",
  target_MAF = target_MAF,
  target_case_prop = 0.5,
  method = "pvalue"
)

# Example 7: Use in optimal weight calculation
PI_estimates <- runif(10, 0.001, 0.1)
M <- diag(10)  # Independent variants
g_identity <- function(x) x

weights <- Optimal_Weights_M(g_identity, B_estimates, PI_estimates, M)
} # }