用R计算多基因病遗传力

可试试网页版的计算工具。

根据Falconer方法, 多基因病的遗传力计算方法如下：

一些教材上使用查表来获取X_g和a_g以及X_r的值，然而现在电脑已经普及，使用电脑可以非常容易的计算这些值，无需再去查表。这里介绍怎么使用R来计算，可以将公式写成函数，简单带入一般人群发病率、患者亲属发病率以及亲缘系数即可一步出结果。

关于正态分布，基本的分布曲线以及累计概率曲线如下图，可以很直观地看到我们所需要求的值使用哪个函数（图中红色曲线为标准正态分布）：

这里不要被这些函数吓到，了解点理论，然后知道怎么用电脑算即可。

从上边的图中可以看出，将易患性分布标准化之后，X实际上是发病率为q时对应的易患性的值，也就是累计概率为1-q时的X值，求X的值，使用的是累计概率曲线函数的反函数，φ^-1, 对应R里的qnorm函数。
比如，q_g为0.1%，那么我们可以这样来求X_g：

q_g <- 0.1/100
X_g <- qnorm(1 - q_g)
X_g

运行示例如下：


同样的办法求X_r即可。

对于a_g的计算，需要先算出X_g，根据X算Z使用的是dnorm函数，示例如下：

q_g <- 0.1/100
X_g <- qnorm(1 - q_g)
z_g <- dnorm(X_g)
a_g <- z_g/q_g

运行示例如下：

现在，我们可以把这些连起来，对于指定的q_g、q_r以及r来计算h²：

q_g <- 0.1/100
q_r <- 22/669    # 669个一级亲属中有22个患病
r   <- 1/2       # 一级亲属的亲缘系数是1/2
X_g <- qnorm(1 - q_g)
z_g <- dnorm(X_g)
a_g <- z_g/q_g
X_r <- qnorm(1 - q_r)
b   <- (X_g - X_r)/a_g
h2  <- b/r   
h2 

这里仅展示一下运行结果：

如果要多次计算，直接把计算过程封装成为函数，那么可以省去很多的步骤，如下

h <- function(q_g, q_r, r=0.5){   #r=0.5是说r的默认值为0.5，这样在算一级亲属的情况的时候，可以不提供这个参数
    X_g = qnorm(1-q_g)
    X_r = qnorm(1-q_r)
    z_g = dnorm(X_g)
    a_g = z_g/q_g
    b   = (X_g - X_r)/a_g
    return(b/r)
}
 
h(0.1/100, 22/669)

看看运行的效果：

这样我们就能够很方便的计算别的情况了，比如：

对于没有一般人群的发病率数据时，使用对照组，比如使用对照组一级亲属数据的情况如下：

试试看，写个函数出来。

---

自2021年9月22日以来，本网页已被访问5567次。