社区首页 >专栏 >表型相关系数与标准误

表型相关系数与标准误

邓飞

发布于 2020-03-31 08:56:32

1.2K00

代码可运行

文章被收录于专栏：育种数据分析之放飞自我育种数据分析之放飞自我

运行总次数：0

代码可运行

一个同学发信给我：

“请问我想计算几个表型数据的相关系数和标准误，如何用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