Main function for computing optimal weights for variant-set tests in the GLOW framework. Computes weights under different hypotheses (H0 vs H1) and approximation methods (normal vs sparse), tailored to the specific test transformation function g.
Optimal_Weights_M(g, Bstar, PI, M, is.posi.wts = TRUE)Function, transformation applied to Z-scores (e.g., identity for Burden, GFisher transformation for Fisher-type tests)
Numeric vector, estimated effect sizes for each variant
Numeric vector, estimated variant-importance scores
Numeric matrix, correlation matrix of marginal Z-scores
Logical, if TRUE, negative weights are set to zero (default TRUE)
A list with optimal weights. The structure depends on the transformation function:
For Burden/Liptak tests (when g is identity):
Weights optimized for Best Effect under H0
Weights optimized for Average Power under Ensemble under H1
For other tests (Fisher, SKAT, etc.):
BE weights using normal approximation
APE weights using normal approximation
BE weights using sparse mixture approximation
APE weights using sparse mixture approximation
This function orchestrates the optimal weight calculation for GLOW tests. The procedure differs based on the transformation function g:
For Burden/Liptak tests (g is identity):
The weights have closed-form solutions: $$wts_{BE} = MU$$ $$wts_{APE} = (I + diag(V)M)^{-1} MU$$ where \(MU = Bstar * PI\) and \(V = Bstar^2 * PI * (1-PI)\).
For other tests (GFisher, SKAT, etc.):
Two approximation methods are used:
Normal Approximation: Assumes Z-scores are approximately normal under the alternative, \(\tilde{Z} \sim N(MMU, M + MVM)\), where:
MMU = M \
MVM = t(M * V) \
The moments of \(g(\tilde{Z}_i)\) are computed on the standardized scale. With
\(s_i = \sqrt{(M + MVM)_{ii}}\) and the correlation matrix
\(R = \mathrm{cov2cor}(M + MVM)\), \(\tilde{Z}_i = MMU_i + s_i U_i\) with
\(U \sim N(0, R)\): get_r_tilde receives the standard deviations
\(s_i\), and CovM_gXgY receives \(R\) together with the scales
\(s_i\), because Mehler's expansion is an expansion in a correlation (bounded by
one, which the truncation at order 64 relies on) with the Hermite coefficients of
\(g(MMU_i + s_i U)\). GLOWr up to 0.3.2, like the legacy code, passed the
covariance \(M + MVM\) into the expansion and its diagonal, a variance, as the
standard deviation; the _N weights of those versions differ from the current ones.
Sparse Mixture Approximation: Explicitly models the sparse signal structure Z_i = Z_0i + Bstar_i * C_i where C_i ~ Bernoulli(PI_i).
Numerical integration and the expansion order:
Every expectation of \(g\) is a one-dimensional integral against the standard normal
density. It is computed by a composite 10-point Gauss-Legendre rule on \([-12, 12]\)
with panels of width 0.1, cut at the point where the argument of \(g\) is zero (the
kink of the two-sided transforms, which depend on \(|x|\)). The Hermite coefficients
\(b_k = E[g(\mu + sU)\,he_k(U)]\) come from one evaluation of \(g\) per distinct
\((\mu, s)\) pair, and Mehler's expansion
\(\mathrm{Cov}\{g(\mu_i + s_iU_i), g(\mu_j + s_jU_j)\} = \sum_k b_k(i)b_k(j)r_{ij}^k\)
runs to order 64. The truncated sum is the covariance of the projected transform, and its
diagonal is the truncated variance; the sparse covariance carries the exact mixture
variance on its diagonal instead, because the mixture formula does not apply to a variant
with itself. Against a direct integration of every pairwise covariance on a 20-variant
window with a near-duplicate pair, the weights of the normal approximation agree to
\(10^{-3}\) of the largest weight and the sparse weights to \(5 \times 10^{-3}\)
(the package tests). GLOWr up to 0.3.3 used integrate() at its default tolerance
and order 8; on that window its _N weights were off by \(9 \times 10^{-2}\).
The rule is validated for the Fisher and polynomial transforms over the scales GLOW meets.
A custom g must be vectorized, square-integrable against the normal density and
smooth except possibly at zero; a transform with another kink, a cusp, a discontinuity or
rapid growth needs its own convergence check.
For each approximation, weights are computed under two hypotheses:
BE (Best under H0): Assumes no signal, optimizes for power under worst-case
APE (Average Power under Ensemble): Accounts for signal, optimizes for average power
Transformation Function Check:
The function checks if g is the identity by comparing function bodies:
identical(body(g), body(function(x) x))
Burden/Liptak: O(n^3) for matrix inversion
Other tests: O(n_pairs * n_nodes * K + 4 * n^2 * K + n^3) where:
n = number of variants
n_pairs = distinct (mean, scale) pairs, at most 2n + 1
n_nodes = quadrature nodes per pair (2,410)
K = number of Hermite polynomial orders (64)
Factor of 4 from the four Mehler sums (the null, the normal-approximation alternative, and two sparse components)
Bstar, PI, and M have compatible dimensions
M is a valid correlation matrix (positive definite, diagonal = 1)
0 <= PI_i <= 1 for all i
Bstar values are positive (standardized effect sizes)
For non-Burden tests, g(mu + s U) is square-integrable under the normal density and smooth except possibly at zero (see the details)
Larger weights indicate variants that contribute more to the test statistic
All weights are normalized by mean(abs(weights)) for comparability
BE weights are conservative (good when signal is uncertain)
APE weights are adaptive (good when signal estimates are reliable)
Sparse weights are more accurate for rare causal variants
Normal weights are computationally simpler approximations
The function uses several optimizations for improved performance:
H0 covariance matrices are cached to speed up repeated calls
Fast path computation for independent variants (identity correlation)
Caching can be controlled: options(GLOWr.use_cache = TRUE/FALSE)
Cache can be cleared: clear_glow_cache()
Cache status: glow_cache_info()
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
if (FALSE) { # \dontrun{
# Example 1: Burden test (identity transformation)
g_identity <- function(x) x
Bstar <- c(0.5, 1.0, 0.3)
PI <- c(0.01, 0.05, 0.02)
M <- diag(3)
weights_burden <- Optimal_Weights_M(g_identity, Bstar, PI, M)
names(weights_burden) # wts_BE, wts_APE
# Example 2: GFisher test
g_GFisher_two <- function(x, df = 2) {
qchisq(log(2) + pnorm(abs(x), lower.tail = FALSE, log.p = TRUE),
df = df, lower.tail = FALSE, log.p = TRUE)
}
weights_fisher <- Optimal_Weights_M(
function(x) g_GFisher_two(x, df = 2),
Bstar, PI, M
)
names(weights_fisher) # wts_BE_N, wts_APE_N, wts_BE_sparse, wts_APE_sparse
} # }