gsem 0.1.3

Genomic Structural Equation Modeling from GWAS summary statistics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

use faer::{Mat, Side};
use gsem_matrix::error::MatrixError;
use gsem_matrix::vech;
use statrs::distribution::{ChiSquared, ContinuousCDF};

/// Compute Q_SNP heterogeneity statistic for a factor.
///
/// Tests whether SNP effects are homogeneous across indicators of a factor.
///
/// Q_SNP = eta' * P * Eig^{-1} * P' * eta
/// where eta = vech(S_subset - Sigma_subset)
pub fn compute_q_snp(
    s_subset: &Mat<f64>,
    sigma_subset: &Mat<f64>,
    v_subset: &Mat<f64>,
) -> Result<(f64, usize, f64), MatrixError> {
    let p = s_subset.nrows();
    let s_vec = vech::vech(s_subset)?;
    let sig_vec = vech::vech(sigma_subset)?;

    let eta: Vec<f64> = s_vec
        .iter()
        .zip(sig_vec.iter())
        .map(|(s, sig)| s - sig)
        .collect();

    let n = eta.len();

    let Ok(eigen) = v_subset.self_adjoint_eigen(Side::Lower) else {
        return Ok((f64::NAN, 0, 1.0));
    };

    let u = eigen.U();
    let s_diag = eigen.S().column_vector();

    let q: f64 = (0..n)
        .filter_map(|i| {
            let eigenval = s_diag[i];
            (eigenval > 1e-10).then(|| {
                let dot: f64 = eta.iter().enumerate().map(|(j, &e)| u[(j, i)] * e).sum();
                dot * dot / eigenval
            })
        })
        .sum();

    let df = p.saturating_sub(1);
    let p_val = if df > 0 && q.is_finite() {
        ChiSquared::new(df as f64)
            .map(|chi2| 1.0 - chi2.cdf(q))
            .unwrap_or(1.0)
    } else {
        1.0
    };

    Ok((q, df, p_val))
}