N2 | 一年多点数据如何计算BLUP值

邓飞

发布于 2021-01-12 15:02:06

2.4K0

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

一年多点数据如何计算BLUP值

上一节，介绍了什么是BLUP值（N1 | 什么是BLUP值？），鸽了这么多天，今天水一篇。

话说，「工欲善其事，必先利其器」，我搞定了Typora写markdown设置免费图库之后（良心教程 | 如何在Typora中设置免费的图床），这写作体验，杠杠的。

为何要用BLUP值

首先是确定模型，BLUP值是随机因子的效应值，所以计算某个因素的BLUP值，将其作为随机因子放到模型中即可。

一般我们所说的BLUP育种值，主要是指个体ID的BLUP值。

对于不考虑系谱关系的个体，将其作为随机因子，计算BLUP值，将其作为排序的依据，当数据出现缺失或者不平衡试验时，BLUP更靠谱。

对于考虑亲缘关系的个体，将其作为随机因子，可以矫正亲缘关系的影响，计算出的BLUP值，更靠谱。

其实，不仅是动物育种里面的动物模型（animal model）使用BLUP值，林木，水产，作物都用BLUP值，使用BLUP值作为品种的排名，比平均值更好。百利而无一害，值得替换。

一年多点数据探索性分析

数据来源于我编写的R包：learnasreml中的MET数据，回头我写篇博客介绍一下这个R包，learnasreml包的安装方法：

if (!requireNamespace("devtools")) install.packages("devtools")

library(devtools)

install_github("dengfei2013/learnasreml")

原数据是多年多点的数据，这里选择一年的数据演示如何操作：

library(learnasreml)
library(tidyverse)
data(MET)
MET %>% filter(Year == 2009) %>% select(-1) -> dat
summary(dat)
str(dat)

数据概况：

> summary(dat)
 Location Rep                    Cul         Yield        
 CI:40    1:50   CalhounGray       :20   Min.   :  7.213  
 FL:40    2:50   CrimsonSweet      :20   1st Qu.: 45.049  
 KN:40    3:50   EarlyCanada       :20   Median : 69.216  
 SC:40    4:50   FiestaF1          :20   Mean   : 68.644  
 TX:40           GeorgiaRattlesnake:20   3rd Qu.: 89.303  
                 Legacy            :20   Max.   :172.065  
                 (Other)           :80   NA's   :2        
> str(dat)
'data.frame': 200 obs. of  4 variables:
 $ Location: Factor w/ 5 levels "CI","FL","KN",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ Rep     : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ Cul     : Factor w/ 10 levels "CalhounGray",..: 3 1 9 2 5 4 7 10 6 8 ...
 $ Yield   : num  56.2 74.2 32.6 74.2 64.8 ...

可以看到，数据有5个地点，每个地点有4个重复，应该是完全随机区组设计，共有10个品种，观测数据为产量。

看一下每个地点的概况：

dat %>% ggplot(.,aes(x = Location,y = Yield,colour = Location)) + geom_boxplot()

可以看到地点FL整体表现最好，TX地点，整体表现最差。

然后在看一下每个品种的表现：

图中，可以看到品种StarbriteF1最好，SugarBaby表现最差。

计算常规的模型

所谓常规的模型，是指通用的普通的模型，上面的数据结构，常用的模型是

**固定因子：**Location + Location:Rep

**随机因子：**Cul + Cul:Location

注意，这里的Rep不是独立的，而是镶嵌在Location中的。如果直接写为：Location + Rep是错误的！

免费的R包：

library(lme4)
m1 =  lmer(Yield ~ Location + Location:Rep +(1|Cul) + (1|Cul:Location),data=dat)
print(m1)
anova(m1)  ####固定因子显著性
ranef(m1)  ####求随机效应的BLUP值
fixef(m1)  ####求固定效应的BLUE值

结果有诡异的地方：boundary (singular) fit: see ?isSingular，这种报错真的很无语，完全不知所云……

> library(lme4)
> m1 =  lmer(Yield ~ Location + Location:Rep +(1|Cul) + (1|Cul:Location),data=dat)
boundary (singular) fit: see ?isSingular
> summary(m1)
Linear mixed model fit by REML ['lmerMod']
Formula: Yield ~ Location + Location:Rep + (1 | Cul) + (1 | Cul:Location)
   Data: dat

REML criterion at convergence: 1623.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.29009 -0.56680 -0.09692  0.56107  2.97208 

Random effects:
 Groups       Name        Variance Std.Dev.
 Cul:Location (Intercept)   0.0     0.00   
 Cul          (Intercept) 150.7    12.27   
 Residual                 370.1    19.24   
Number of obs: 198, groups:  Cul:Location, 50; Cul, 10

Fixed effects:
                Estimate Std. Error t value
(Intercept)       61.740      7.509   8.222
LocationFL        31.342      8.851   3.541
LocationKN        -3.344      8.851  -0.378
LocationSC        21.063      8.851   2.380
LocationTX        -4.564      8.851  -0.516
LocationCI:Rep2   -6.237      8.851  -0.705
LocationFL:Rep2   19.035      8.604   2.212
LocationKN:Rep2    3.953      8.604   0.459
LocationSC:Rep2   -1.138      8.604  -0.132
LocationTX:Rep2  -14.340      8.604  -1.667
LocationCI:Rep3  -16.180      8.851  -1.828
LocationFL:Rep3    1.086      8.604   0.126
LocationKN:Rep3    6.805      8.604   0.791
LocationSC:Rep3   -9.526      8.604  -1.107
LocationTX:Rep3  -26.521      8.604  -3.083
LocationCI:Rep4   -8.315      8.851  -0.940
LocationFL:Rep4    3.542      8.604   0.412
LocationKN:Rep4   19.642      8.604   2.283
LocationSC:Rep4   10.983      8.604   1.277
LocationTX:Rep4  -28.289      8.851  -3.196

Correlation matrix not shown by default, as p = 20 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

结果可以看到，应该是品种和地点的交互方差组分为0Cul:Location (Intercept) 0.0 0.00，才会出现上面的报警。

用商业包：asreml比较一下：

> library(asreml)
> m2 = asreml(Yield ~ Location/Rep, random = ~ Cul + Cul:Location, data=dat)
Model fitted using the gamma parameterization.
ASReml 4.1.0 Sat Jan  9 20:18:23 2021
          LogLik        Sigma2     DF     wall    cpu
 1      -653.007       384.055    178 20:18:23    0.0
 2      -650.251       389.806    178 20:18:23    0.0 (1 restrained)
 3      -648.890       383.009    178 20:18:23    0.0 (1 restrained)
 4      -648.233       373.708    178 20:18:23    0.0 (1 restrained)
 5      -648.163       370.592    178 20:18:23    0.0 (1 restrained)
 6      -648.162       370.130    178 20:18:23    0.0
> summary(m2)$varcomp
                component std.error  z.ratio bound %ch
Cul          1.506628e+02  79.77070 1.888699     P 0.1
Cul:Location 9.715953e-05        NA       NA     B 0.0
units!R      3.701301e+02  40.26423 9.192528     P 0.0

就是品种与地点交互的方差组分为0。

比较一下品种的BLUP值：完全一致。

高级一点的模型

上面的普通模型，假定地点的方差一致，实际情况不知道，所以我们可以用更通用的模型去拟合，这样无论地点的方差是否一致，该模型都是合适的。lme4包不能定义残差的结构，下面用asreml进行演示。

这里定义地点间残差异质的函数是：dsum(units|Location)，

地点同质的矩阵结构：

用矩阵表示是：

\begin{bmatrix} \sigma_a^2& 0& 0\\ 0 &\sigma_a^2 &0\\ 0 &0& \sigma_a^2\\ \end{bmatrix}

这里

\sigma_a^2

为每个地点的残差方差，地点的残差方差都是一样的，可以进行联合方差分析。

地点异质的矩阵结构：

\begin{bmatrix} \sigma_a^2& 0& 0\\ 0 &\sigma_b^2 &0\\ 0 &0& \sigma_c^2\\ \end{bmatrix}

在地点异质的模型中，每个地点都会估计出一个方差组分，代码如下：

library(asreml)
m3 = asreml(Yield ~ Location/Rep, random = ~ Cul + Cul:Location, residual = ~ dsum(~units|Location),data=dat)
summary(m3)$varcomp

结果：

> m3 = update(m3)
Multi-section model using the sigma parameterization.
ASReml 4.1.0 Sat Jan  9 20:32:20 2021
          LogLik        Sigma2     DF     wall    cpu
 1      -647.287           1.0    178 20:32:20    0.0
 2      -647.287           1.0    178 20:32:20    0.0
> summary(m3)$varcomp
                 component std.error  z.ratio bound %ch
Cul           1.420246e+02  75.65677 1.877222     P   0
Cul:Location  2.370980e-05        NA       NA     B   0
Location_KN!R 3.049822e+02  75.35989 4.047010     P   0
Location_CI!R 3.238231e+02  80.74030 4.010675     P   0
Location_SC!R 4.604491e+02 112.08537 4.108021     P   0
Location_FL!R 3.821732e+02  94.01243 4.065135     P   0
Location_TX!R 3.830825e+02  94.86128 4.038344     P   0

可以看到，不同地点的方差组分为：

KN：304.982
CI：323.8
……

地点间方差是否同质呢？`LRT检验`

asreml软件中，提供了LRT检测的程序：lrt.asreml

> lrt.asreml(m2,m3)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

      df LR-statistic Pr(Chisq)
m3/m2  1   8.9145e-06    0.4988

可以看到，P值为0.4988，不显著，说明地点齐次，将所有地点假定齐次的模型m2没有问题。否者，m3模型会更好。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-01-09，如有侵权请联系 cloudcommunity@tencent.com 删除

网站

本文分享自育种数据分析之放飞自我微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

网站

登录后参与评论

0 条评论

热度