R/get_B_model_selection.R
select_best_model.RdSelects the best linear regression model for predicting effect size B based on minor allele frequency (MAF). This function fits 8 different models with various transformations of X (MAF) and Y (effect size measure) and selects the best model based on the specified criterion.
select_best_model(
X,
Y,
criterion = "R2",
custom_models = NULL,
return_all = FALSE,
verbose = TRUE
)Numeric vector of minor allele frequencies (MAF values)
Numeric vector of effect size measures (typically BETA^2 or related)
Character string specifying selection criterion:
"R2" (default): R-squared, proportion of variance explained
"adj_R2": Adjusted R-squared, penalizes model complexity
"CV_R2": Cross-validated R-squared (leave-one-out)
See Details section for guidance on choosing a criterion. Note: "CV" is deprecated but still accepted (use "CV_R2" instead). AIC and BIC have been removed - see Details for explanation
Optional list of custom model formulas to include in candidate set. Each element should be a formula object
Logical. If TRUE, returns list of all fitted models with their fit statistics. If FALSE (default), returns only the best model
Logical. If TRUE (default), prints detailed model comparison table and selection summary. Set to FALSE to suppress output in loops or simulations. Default is TRUE for transparency in model selection.
If return_all=FALSE: A fitted linear model object (class "lm") representing the best model based on the specified criterion. If return_all=TRUE: A list with elements:
best_model: The selected fitted lm object
best_name: Character name of the selected model
all_models_info: A list with
models: A named list with one element per candidate model, each a list containing model (the fitted lm), formula, formula_text, R2, adj_R2, CV_R2, sigma, and coefficients
summary_table: A data frame with columns model_name, formula, R2, adj_R2, CV_R2, and is_best (one row per model)
criterion_used: The criterion used for selection
training_data: A list with the original X and Y
The function considers the following transformations:
X transformations: X, f(X) = X(1-X), log(X), log(f(X))
Y transformations: Y, log(Y)
This results in 8 possible models:
Y ~ X
Y ~ f(X)
Y ~ log(X)
Y ~ log(f(X))
log(Y) ~ X
log(Y) ~ f(X)
log(Y) ~ log(X)
log(Y) ~ log(f(X))
The model with the best criterion value is returned. This approach captures the relationship between MAF and effect size, which can vary across different genetic architectures.
Selection Criteria - When to Use Which:
What it measures: Proportion of variance explained in training data
Higher is better
Best for: Large samples (n > 100), exploratory analysis
Advantages: Fast, intuitive, widely used
Disadvantages: Can overfit with small samples, always favors complex models
Default for: Historical compatibility
What it measures: R^2 adjusted for number of predictors, penalizes complexity
Higher is better
Best for: Small-medium samples (n < 100), comparing models with different complexity
Advantages: Fast, balances fit and parsimony, more conservative than R^2
Disadvantages: Still uses in-sample fit (not true out-of-sample validation)
Good choice: When you want simplicity preference without CV computational cost
Comparison with R^2: Will prefer simpler models when fit is similar
What it measures: Out-of-sample prediction accuracy (LOOCV)
Higher is better
Best for: Small-medium samples (10 < n < 200), when generalization matters
Advantages: Most robust against overfitting, best generalization estimate
Disadvantages: Slower (n model fits per candidate), can be noisy with very small n
Recommended for: B estimation from literature data (typical n=20-100)
Quick Guide:
Small sample (n < 50): Use CV_R2 (most robust)
Medium sample (50 < n < 200): Use CV_R2 or adj_R2
Large sample (n > 200): Any criterion works, R^2 is fastest
When in doubt: Use CV_R2 (most robust)
For speed: Use adj_R2 (no iteration)
Balance of speed and robustness: Use adj_R2
AIC and BIC are not available as selection criteria because they are not appropriate for comparing models with different response transformations (e.g., Y vs log(Y) vs sqrt(Y)). Use R^2, adj_R^2, or CV_R^2 instead, which properly account for scale differences between transformations.
R^2, adj_R2: O(k) where k = number of candidate models (typically 8)
CV_R2: O(k * n) where n = sample size (LOOCV requires n fits per model)
For typical use: n=50, k=8 -> ~400 model fits for CV_R2. Fast with lm() (< 1 second)
X and Y have the same length
X values are in (0, 1) for MAF
Y values are positive for log transformation
X values are not exactly 0 or 1 for log transformations
get_B for B estimation using the selected model
compare_B_models for comparing fitted B models
# Simulate MAF and effect sizes
set.seed(123)
MAF <- runif(100, 0.01, 0.5)
BETA <- sqrt(0.5 * MAF * (1 - MAF)) + rnorm(100, 0, 0.1)
Y <- BETA^2
# Select best model using R^2 (default, fastest)
best_model_r2 <- select_best_model(MAF, Y)
#>
#> ========================================
#> Model Selection for B Estimation
#> ========================================
#> Fitting 8 candidate models...
#>
#> Model Comparison:
#> ID Formula R^2 Adj.R^2 CV R^2
#> ---- -------------------- -------- -------- ----------
#> 1 Y ~ X 0.2363 0.2286 ---
#> 2 Y ~ fX 0.2359 0.2281 ---
#> 3 Y ~ logX 0.2370 0.2292 ---
#> 4 Y ~ logfX 0.2191 0.2111 ---
#> 5 logY ~ X 0.2982 0.2910 ---
#> 6 logY ~ fX 0.3486 0.3420 ---
#> 7 logY ~ logX 0.3993 0.3932 ---
#> *8 logY ~ logfX 0.4072 0.4011 --- <- SELECTED
#>
#> Selection Criterion: R^2
#> Best Model: logY ~ logfX (R^2 = 0.4072)
#> ========================================
#>
summary(best_model_r2)
#>
#> Call:
#> lm(formula = fmla, data = environment())
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.7513 -0.4190 0.1353 0.4814 1.9403
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.6690 0.2711 -2.468 0.0153 *
#> logfX 1.1085 0.1351 8.205 9.27e-13 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 0.8357 on 98 degrees of freedom
#> Multiple R-squared: 0.4072, Adjusted R-squared: 0.4011
#> F-statistic: 67.32 on 1 and 98 DF, p-value: 9.268e-13
#>
# Select using CV_R2 (recommended for small samples)
best_model_cv <- select_best_model(MAF, Y, criterion = "CV_R2")
#>
#> ========================================
#> Model Selection for B Estimation
#> ========================================
#> Fitting 8 candidate models...
#>
#> Model Comparison:
#> ID Formula R^2 Adj.R^2 CV R^2
#> ---- -------------------- -------- -------- ----------
#> 1 Y ~ X 0.2363 0.2286 0.2048
#> 2 Y ~ fX 0.2359 0.2281 0.2084
#> 3 Y ~ logX 0.2370 0.2292 0.2104
#> 4 Y ~ logfX 0.2191 0.2111 0.1921
#> 5 logY ~ X 0.2982 0.2910 0.2645
#> 6 logY ~ fX 0.3486 0.3420 0.3138
#> 7 logY ~ logX 0.3993 0.3932 0.3688
#> *8 logY ~ logfX 0.4072 0.4011 0.3776 <- SELECTED
#>
#> Selection Criterion: CV R^2
#> Best Model: logY ~ logfX (CV R^2 = 0.3776)
#> ========================================
#>
# Select using adjusted R^2 (balances fit and simplicity)
best_model_adjr2 <- select_best_model(MAF, Y, criterion = "adj_R2")
#>
#> ========================================
#> Model Selection for B Estimation
#> ========================================
#> Fitting 8 candidate models...
#>
#> Model Comparison:
#> ID Formula R^2 Adj.R^2 CV R^2
#> ---- -------------------- -------- -------- ----------
#> 1 Y ~ X 0.2363 0.2286 ---
#> 2 Y ~ fX 0.2359 0.2281 ---
#> 3 Y ~ logX 0.2370 0.2292 ---
#> 4 Y ~ logfX 0.2191 0.2111 ---
#> 5 logY ~ X 0.2982 0.2910 ---
#> 6 logY ~ fX 0.3486 0.3420 ---
#> 7 logY ~ logX 0.3993 0.3932 ---
#> *8 logY ~ logfX 0.4072 0.4011 --- <- SELECTED
#>
#> Selection Criterion: Adj. R^2
#> Best Model: logY ~ logfX (Adj. R^2 = 0.4011)
#> ========================================
#>
# Suppress verbose output (default is verbose = TRUE)
best_model_quiet <- select_best_model(MAF, Y, criterion = "CV_R2", verbose = FALSE)
# Get all models for comparison (returns a richer structure, see Value)
all_models <- select_best_model(MAF, Y, return_all = TRUE)
#>
#> ========================================
#> Model Selection for B Estimation
#> ========================================
#> Fitting 8 candidate models...
#>
#> Model Comparison:
#> ID Formula R^2 Adj.R^2 CV R^2
#> ---- -------------------- -------- -------- ----------
#> 1 Y ~ X 0.2363 0.2286 0.2048
#> 2 Y ~ fX 0.2359 0.2281 0.2084
#> 3 Y ~ logX 0.2370 0.2292 0.2104
#> 4 Y ~ logfX 0.2191 0.2111 0.1921
#> 5 logY ~ X 0.2982 0.2910 0.2645
#> 6 logY ~ fX 0.3486 0.3420 0.3138
#> 7 logY ~ logX 0.3993 0.3932 0.3688
#> *8 logY ~ logfX 0.4072 0.4011 0.3776 <- SELECTED
#>
#> Selection Criterion: R^2
#> Best Model: logY ~ logfX (R^2 = 0.4072)
#> ========================================
#>
# The ready-made summary table has R2, adj_R2, and CV_R2 for every candidate
print(all_models$all_models_info$summary_table)
#> model_name formula R2 adj_R2 CV_R2 is_best
#> 1 Y ~ X Y ~ X 0.2363465 0.2285542 0.2047647 FALSE
#> 2 Y ~ f(X) Y ~ fX 0.2359031 0.2281062 0.2083588 FALSE
#> 3 Y ~ log(X) Y ~ logX 0.2370212 0.2292357 0.2103874 FALSE
#> 4 Y ~ log(f(X)) Y ~ logfX 0.2190681 0.2110995 0.1920675 FALSE
#> 5 log(Y) ~ X logY ~ X 0.2981762 0.2910148 0.2645285 FALSE
#> 6 log(Y) ~ f(X) logY ~ fX 0.3486218 0.3419750 0.3138277 FALSE
#> 7 log(Y) ~ log(X) logY ~ logX 0.3993027 0.3931731 0.3687912 FALSE
#> 8 log(Y) ~ log(f(X)) logY ~ logfX 0.4071924 0.4011433 0.3776398 TRUE
# Equivalently, build the comparison from the per-model list
model_comparison <- sapply(all_models$all_models_info$models, function(m) {
c(R2 = m$R2, adj_R2 = m$adj_R2, CV_R2 = m$CV_R2)
})
print(model_comparison)
#> Y ~ X Y ~ f(X) Y ~ log(X) Y ~ log(f(X)) log(Y) ~ X log(Y) ~ f(X)
#> R2 0.2363465 0.2359031 0.2370212 0.2190681 0.2981762 0.3486218
#> adj_R2 0.2285542 0.2281062 0.2292357 0.2110995 0.2910148 0.3419750
#> CV_R2 0.2047647 0.2083588 0.2103874 0.1920675 0.2645285 0.3138277
#> log(Y) ~ log(X) log(Y) ~ log(f(X))
#> R2 0.3993027 0.4071924
#> adj_R2 0.3931731 0.4011433
#> CV_R2 0.3687912 0.3776398
# Add custom model
custom_formula <- formula(Y ~ poly(X, 2)) # Polynomial model
best_custom <- select_best_model(MAF, Y,
custom_models = list(custom_formula))
#>
#> ========================================
#> Model Selection for B Estimation
#> ========================================
#> Fitting 9 candidate models...
#>
#> Model Comparison:
#> ID Formula R^2 Adj.R^2 CV R^2
#> ---- -------------------- -------- -------- ----------
#> 1 Y ~ X 0.2363 0.2286 ---
#> 2 Y ~ fX 0.2359 0.2281 ---
#> 3 Y ~ logX 0.2370 0.2292 ---
#> 4 Y ~ logfX 0.2191 0.2111 ---
#> 5 logY ~ X 0.2982 0.2910 ---
#> 6 logY ~ fX 0.3486 0.3420 ---
#> 7 logY ~ logX 0.3993 0.3932 ---
#> *8 logY ~ logfX 0.4072 0.4011 --- <- SELECTED
#> 9 Y ~ poly(X, 2) 0.2398 0.2241 ---
#>
#> Selection Criterion: R^2
#> Best Model: logY ~ logfX (R^2 = 0.4072)
#> ========================================
#>