Statistical formula

🎯 Summary-of-equations

Genotype-Properties

1.1 MAF estimate for SNP $j$:

• Meaning: the estimated frequency of the MAF for SNP $j$, assuming that genotypes are coded as 0/1/2.
• R script:

p_hat_j = mean(x[, j]) / 2

1.2 Expected MAF sampling variance:

• Meaning: the estimated variance of MAF under sample size $N$.
• R script:

se2_p_hat_j = p_hat_j*(1 - p_hat_j) / (2 * N)

1.3 Expected variance of a SNP:

• Meaning: the genotype variance of SNP $j$, theoretically, consistent with the variance of the variable encoded as 0/1/2.
• If we set $p$ as MAF and $q$ as another allele frequency, then $var(X_j) = 2\hat p_j * \hat q_j$
• R script:

var_x_j = 2 * p_hat_j * (1 - p_hat_j)

1.4 LD matrix:

• Meaning: the normalized SNP covariance matrix (LD matrix).
• R script:

L = t(x) %*% x / N

SMR

SMR 里 HEIDI heterogeneity：

1）符号与输入（你手里有什么）

假设某个基因 cis 区域有 n 个 SNP（nsnp = n），SNP 索引 $(k = 1,\dots,n)$。

对于每个 SNP $(k)$：

• exposure（eQTL）summary：效应 $( \hat b_{x,k} )$，标准误 $( \hat s_{x,k} )$
• outcome（GWAS）summary：效应 $( \hat b_{y,k} )$，标准误 $( \hat s_{y,k} )$
• LD 相关矩阵 $(R)$：元素 $(R_{kl} = r_{kl})$

代码层面常见对应：

• bx[k] = expo[k,1], sx[k] = expo[k,2]
• by[k] = outco[k,1], sy[k] = outco[k,2]
• R = corr（n x n）
• lead = which.max(expo[,3])（最强 eQTL 作为 lead SNP）

2）SMR ratio（每个 SNP 一个比值）

SMR 的核心比值（Wald ratio）：

• 每个 SNP $k$ 的 ratio：

把所有 SNP 的 ratio 堆成向量：

形状检查：

• bxy：长度 n
• Sigma_b = Cov(bxy)：n x n

3）SMR ratio 的 SE（每个 SNP 一个标准误）

前面已经定义了每个 SNP $k$ 的 ratio：

那么他对应的抽样方差（delta method 一阶近似）则推导如下

把 $u=\hat b_{y,k}$、$v=\hat b_{x,k}$、$f(u,v)=u/v$，则：

因此一般形式（保留 $Cov(\hat b_{y,k},\hat b_{x,k})$）：

SMR 常用两样本设定（GWAS 与 eQTL 来自独立样本）下，近似 $Cov(\hat b_{y,k},\hat b_{x,k})\approx 0$，于是：

最终：

代码层面常见对应：

• bxy[k] = by[k] / bx[k]
• se_bxy[k] = sqrt( sy[k]^2 / bx[k]^2 + (by[k]^2 * sx[k]^2) / bx[k]^4 )

形状检查：

• bxy：长度 n
• se_bxy：长度 n

备注：

• 如果 $\hat b_{x,k}$ 很小，ratio 和 SE 会被放大（弱工具变量问题），实际流程通常会先按 eQTL 显著性/强度筛 SNP。

4）基于 SMR 框架：SMR 输出的 b / se 代表什么？

SMR 把 cis-eQTL top SNP 当作工具变量 $z$，把表达当 exposure $x$，把性状当 outcome $y$。

（a）SMR 的效应估计 b（输出里的 b_SMR）

SMR 默认用 lead SNP 的 Wald ratio 作为 SMR 效应估计：

它表示的是：在工具变量 $z$ 支撑下，表达 $x$ 每增加 1 个单位（以 eQTL 数据中表达的刻度为单位），性状 $y$ 改变多少。

单位怎么理解取决于你输入 summary 的定义：

• 如果 outcome 是定量性状：$\hat b_{y,0}$ 是性状单位上的效应，则 $b_{SMR}$ 是“表达单位 $\rightarrow$ 性状单位”的斜率估计。
• 如果 outcome 是病例-对照且 GWAS 用 $\log(OR)$：则 $b_{SMR}$ 是“表达单位 $\rightarrow$ \log(OR)”的斜率估计（可指数化解释为 OR 的倍数变化）。

关键点：

• b_SMR 不是“某个等位基因对性状的效应”，而是把 SNP 作为工具变量后换算出来的“表达对性状”的效应刻度。
• 在只用单一工具变量（lead SNP）时，这个估计反映的是一种“遗传支撑下的表达-性状关联”（常被称为 pleiotropic association），并不自动等价于“严格因果”。

（b）SMR 的标准误 se（输出里的 se_SMR）

se_SMR 就是同一个 lead SNP ratio 的标准误：

它代表的是：$\hat b_{y,0}$（GWAS）和 $\hat b_{x,0}$（eQTL）的抽样误差，经过 ratio（除法）传播后，$b_{SMR}$ 的抽样不确定性。

（c）SMR 的显著性检验（输出里的 p_SMR）

SMR 用 1 df 的 Wald/chi-square 检验（两种写法等价）：

解释：

• p_SMR 小：说明在 lead SNP 作为工具变量的意义下，表达与性状存在显著的遗传支撑关联（值得进一步用 HEIDI 排除 linkage）。
• p_SMR 大：说明在该工具变量下，没有足够证据支持表达-性状的遗传关联。

（后面的 HEIDI 就是在这个基础上，用多个 SNP 的 ratio 一致性来检验是否更像“同一因果信号”还是“linkage/多信号混合”。）

HEIDI

1）HEIDI 的核心：差值向量 d（用来测“异质性”）

先选一个 lead SNP（索引记作 $0$，代码里就是 lead）。

对每个非 lead SNP $i \neq 0$，定义差值：

把所有 $i \neq 0$ 的差值堆成向量：

其中 $m = n-1$

直觉：

• 如果一个基因区域里所有 SNP 反映的是同一个因果信号，那么所有 $\hat b_{xy,i}$ 应该一致
  => $\hat d_i \approx 0$（只剩抽样误差）
• 如果出现连锁/多信号混合，不同 SNP 的 $\hat b_{xy}$ 会系统性不一致
  => $\hat d$ 会明显偏离 0（异质性）

形状检查：

• d：长度 m = n-1
• V = Cov(d)：m x m

2）MVN 是什么？为什么会出现？

MVN = Multivariate Normal，多元正态分布。

HEIDI 使用近似：

含义：

• $\mathbf d = E(\hat{\mathbf d})$：差值向量的真实均值
• $V = Cov(\hat{\mathbf d})$：差值向量的协方差矩阵

零假设（无异质性）：

为什么 $V$ 不是对角矩阵？

• SNP 之间常有 LD（ $R$ 非对角不为 0）
• 这会导致不同 SNP 的估计误差相关
  => $\hat b_{xy,i}$ 之间相关
  => $\hat d$ 的不同分量也相关
  => 必须用完整协方差 $V$，不能只用方差（对角线）

3）最关键：怎么从 Sigma_b = Cov(bxy) 构造 V = Cov(d)？

先定义：

元素是 $\Sigma_b[k,l] = Cov(\hat b_{xy,k}, \hat b_{xy,l})$

因为 $\hat d_i = \hat b_{xy,i} - \hat b_{xy,0}$，所以对任意非 lead SNP $i,j \neq 0$：

用协方差恒等式 $Cov(A-B, C-B)=Cov(A,C)+Var(B)-Cov(A,B)-Cov(C,B)$ 得到：

这就是构造 V=covdxy 的根本原因：d 是 bxy 的线性变换，协方差可以直接从 Sigma_b 推出来。

4）HEIDI 统计量（概念上怎么用 V）

常见二次型（不显式求逆，用 solve 更稳）：

近似下：

• $Q \sim \chi^2_\nu$（ν 近似为 m，实际实现可能因筛 SNP/约束略有调整）

解释：

• p 小（Q 大）：d 偏离 0 明显 => 异质性显著 => 更像连锁/多信号 => 通常拒绝该 SMR 关联
• p 大：一致性较好 => 支持单信号解释（注意：不等价于“证明因果”）

Eigen Decomposition

1.基本概念

1.1 数学定义 对于实对称矩阵 $R$（如LD矩阵），存在特征值分解：

其中：

• $U$：特征向量矩阵（正交矩阵）
• $\Lambda$：特征值对角矩阵
• $U^T$：$U$ 的转置

1.2 各成分解释

2.几何意义

变换三部曲

1. $U^T \times \text{向量}$：旋转到特征向量坐标系

2. $\Lambda \times (\text{结果})$：沿新坐标轴缩放

3. $U \times (\text{结果})$：旋转回原始坐标系

3. R语言实现

3.1 基础分解

# 生成对称矩阵
set.seed(123)
m <- 5
R <- matrix(runif(m^2, 0.1, 0.9), m, m)
R <- (R + t(R)) / 2  # 强制对称
diag(R) <- 1         # 对角线设为1（LD矩阵）

# 特征值分解
eig <- eigen(R, symmetric = TRUE)

# 提取结果
U <- eig$vectors      # 特征向量矩阵
values <- eig$values  # 特征值向量
Lambda <- diag(values) # 特征值对角矩阵

# 验证分解正确性
R_recon <- U %*% Lambda %*% t(U)
all.equal(R, R_recon, tolerance = 1e-10)  # 应返回TRUE

3.2 降维近似（低秩近似）

# 只保留前k个主成分
k <- 2  # 保留的主成分数量

# 选择前k个特征值和特征向量
U_k <- U[, 1:k]                # m × k
values_k <- values[1:k]        # 长度k
Lambda_k <- diag(values_k)     # k × k

# 低秩近似重建
R_approx <- U_k %*% Lambda_k %*% t(U_k)  # m × m，但秩为k

4. 遗传学中的应用

4.1 LD矩阵分析

# LD矩阵的特征值分解
analyze_ld_structure <- function(R, eigen_ratio = 0.95) {
    eig <- eigen(R, symmetric = TRUE)
    values <- eig$values
    U <- eig$vectors
    
    # 计算累积方差解释比例
    total_var <- sum(values)
    cum_prop <- cumsum(values) / total_var
    
    # 确定保留的主成分数
    k <- which(cum_prop >= eigen_ratio)[1]
    if (is.na(k)) k <- length(values)
    
    # 输出结果
    cat("总方差:", total_var, "\n")
    cat("特征值分布:\n")
    print(head(values, 10))
    cat("\n保留主成分数:", k, "\n")
    cat("解释方差比例:", cum_prop[k], "\n")
    
    return(list(U = U, values = values, k = k))
}

4.2 高效计算技巧

# 技巧1：计算 R × beta（避免显式构造R）
compute_R_times_beta <- function(beta, U, values) {
    # beta: m×1向量
    # U: m×k特征向量矩阵
    # values: k个特征值
    
    # 高效计算：R × beta = U × diag(values) × U^T × beta
    tmp <- crossprod(U, beta)  # = U^T × beta，k×1矩阵
    tmp <- tmp[, 1] * values   # 转换为向量并乘以特征值
    result <- U %*% tmp        # = U × (diag(values) × U^T × beta)
    
    return(result)
}

# 技巧2：计算 R^{-1} × beta（避免求逆）
compute_R_inv_times_beta <- function(beta, U, values, epsilon = 1e-8) {
    # 过滤小特征值，避免数值问题
    keep <- values > epsilon
    U_filtered <- U[, keep]
    values_filtered <- values[keep]
    
    # 计算：R^{-1} × beta = U × diag(1/values) × U^T × beta
    tmp <- crossprod(U_filtered, beta)
    tmp <- tmp[, 1] / values_filtered  # 关键：除以特征值！
    result <- U_filtered %*% tmp
    
    return(result)
}

5. 数值稳定性处理

5.1 正定性校正

# 确保矩阵正定
ensure_positive_definite <- function(R, eig_tol = 1e-8) {
    eig <- eigen(R, symmetric = TRUE)
    values <- eig$values
    U <- eig$vectors
    
    # 检查最小特征值
    min_eig <- min(values)
    if (min_eig < eig_tol) {
        cat(sprintf("警告：最小特征值 = %.2e，进行校正\n", min_eig))
        
        # 方法1：添加岭惩罚
        lambda <- abs(min_eig) + 1e-6
        R_corrected <- R + diag(lambda, nrow(R))
        
        # 重新标准化（如果是相关矩阵）
        s <- 1 / sqrt(diag(R_corrected))
        R_corrected <- diag(s) %*% R_corrected %*% diag(s)
        
        return(R_corrected)
    }
    
    return(R)
}

5.2 nearPD 校正

# 使用Matrix包的nearPD函数
if (!require(Matrix)) install.packages("Matrix")
library(Matrix)

correct_with_nearpd <- function(R) {
    # 寻找最近的正定矩阵
    R_pd <- nearPD(R, 
                   corr = TRUE,      # 保持相关矩阵结构
                   keepDiag = TRUE,  # 保持对角线
                   do2eigen = TRUE,  # 确保正定
                   eig.tol = 1e-6,   # 特征值容忍度
                   conv.tol = 1e-7,  # 收敛容忍度
                   maxit = 100)      # 最大迭代次数
    
    return(as.matrix(R_pd$mat))
}

6. 总结表格