表型相关系数与标准误
一个同学发信给我：
“请问我想计算几个表型数据的相关系数和标准误，如何用R语言操作？”
我回答：“R中默认的函数有cor计算相关系数，标准误的话估计要用重抽样去操作？，但是很少有人会计算标准误这个数值。”
于是，该同学发过来一个截图：
上三角为表型相关系数，下三角为表型相关系数的标准误。
于是，我陷入了沉思，这里的表型相关系数，应该是和遗传相关系数对应的，属于数理遗传学的范畴，而我说的相关系数，应该是统计学的范畴。
然后我告诉该同学，你等着，我写篇公众号，解释一下这个误区。
统计学范畴的相关系数计算方法：
首先我们生成两个数据：
> set.seed(123)
> x=rnorm(1000)
> y= 0.3*x + rnorm(1000)
> dat=cbind(x,y)
> head(dat)
x y
[1,] -0.56047565 -1.1639414
[2,] -0.23017749 -1.1090083
[3,] 1.55870831 0.4496323
[4,] 0.07050839 -0.1110226
[5,] 0.12928774 -2.5105565
[6,] 1.71506499 1.5550930然后使用cor函数计算x和y的相关系数：
> cor(dat)
x y
x 1.0000000 0.3573148
y 0.3573148 1.0000000检验x和y的相关系数的显著性：
> cor.test(dat)
Error in cor.test.default(dat) : 缺少参数"y",也没有缺省值
> cor.test(dat$x,dat$y)
Pearson's product-moment correlation
data: dat$x and dat$y
t = 12.086, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3020116 0.4102211
sample estimates:
cor
0.3573148 可以看出，两者显著性达到极显著。
数量遗传学范畴的相关系数计算方法：
一般是混合线性模型，有随机因子，比如父本，比如家系，比如个体动物模型，然后计算性状间的方差组分和协方差组分，最后计算遗传相关系数和表型相关系数。
> library(asreml)
> data("harvey")
> head(harvey)
Calf Sire Dam Line ageOfDam y1 y2 y3
1 101 Sire_1 0 1 3 192 390 224
2 102 Sire_1 0 1 3 154 403 265
3 103 Sire_1 0 1 4 185 432 241
4 104 Sire_1 0 1 4 183 457 225
5 105 Sire_1 0 1 5 186 483 258
6 106 Sire_1 0 1 5 177 469 267这里，将Sire作为随机因子，分析y1，y3两个性状的遗传相关和表型相关，以及他们的标准误。
> mod = asreml(cbind(y1,y3) ~ trait + trait:Line ,
+ random = ~ us(trait):Sire, residual = ~ units:us(trait), data=harvey)
Model fitted using the sigma parameterization.
ASReml 4.1.0 Sun Mar 29 17:46:47 2020
LogLik Sigma2 DF wall cpu
1 -437.786 1.0 124 17:46:47 0.0
2 -433.479 1.0 124 17:46:47 0.0
3 -429.450 1.0 124 17:46:47 0.0
4 -427.635 1.0 124 17:46:47 0.0
5 -427.223 1.0 124 17:46:47 0.0
6 -427.206 1.0 124 17:46:47 0.0
7 -427.206 1.0 124 17:46:47 0.0
> summary(mod)$varcomp
component std.error z.ratio bound %ch
trait:Sire!trait_y1:y1 27.20936 26.59364 1.0231532 P 0
trait:Sire!trait_y3:y1 -12.81265 41.71667 -0.3071350 P 0
trait:Sire!trait_y3:y3 124.88915 125.13874 0.9980055 P 0
units:trait!R 1.00000 NA NA F 0
units:trait!trait_y1:y1 132.36803 25.00112 5.2944842 P 0
units:trait!trait_y3:y1 -59.97696 39.91635 -1.5025662 P 0
units:trait!trait_y3:y3 647.80434 122.32450 5.2957858 P 0注意，这里的方差组分是component一列，前三列为两个性状的遗传方差组分，后三列为环境方差组分。
遗传相关 = cov_sire_12/sqrt(V_sire_11*V_sire_22)
表型相关 = (cov_sire_12 + cov_units_12)/sqrt((V_sire_11 + V_units_11)*(V_sire_22 + V_units_22))
因此相关的计算方法：
> # 遗传相关
> vpredict(mod, rg ~ V2/sqrt(V1*V3))
Estimate SE
rg -0.2197948 0.6668593
> # 表型相关
> vpredict(mod, rp ~ (V2 + V6)/sqrt((V1+V5)*(V3+V7)))
Estimate SE
rp -0.2072908 0.1434625如果只是单纯的使用统计的相关系数，计算如下：
> # 相关系数：统计的方法
> cor(harvey$y1,harvey$y3)
[1] -0.3208209
> cor.test(harvey$y1,harvey$y3)
Pearson's product-moment correlation
data: harvey$y1 and harvey$y3
t = -2.6886, df = 63, p-value = 0.009171
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.52373856 -0.08345172
sample estimates:
cor
-0.3208209可以看到结果偏高，所以注意，如果你是分析的育种数据，还是要按照上面数量遗传的方法，通过计算方差组分和协方差组分，计算出的遗传相关和表型相关才是正确的。
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2020-03-29，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读