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

数据分析 农业技术
邓飞    2021-01-13    1093

一年多点数据如何计算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()

图片.png

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

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

图片.png

图中,可以看到品种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值:完全一致。

图片.png

高级一点的模型

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

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

地点同质的矩阵结构:

用矩阵表示是:

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

地点异质的矩阵结构:

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

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模型会更好。

本文来自微信公众号【育种数据分析之放飞自我】公众号ID:R-breeding;未经许可谢绝二次转载至其他网站。

本文为专栏作者授权科易网发表,版权归原作者所有。文章系作者个人观点,不代表科易网立场,转载请联系原作者。如有任何疑问,请联系ky@1633.com。
热门观点