DMU-单性状重复力模型-学习笔记3

单性状重复力模型

本次主要是演示如何使用DMU分析单性状重复力模型.

重复力模型和动物模型的区别:
不是所有的性状都可以分析重复力模型, 首先重复力模型是动物模型的拓展, 它适合一个个体多个观测值的情况.

  • 比如猪的产仔数, 一个母猪有多个胎次
  • 比如鸡的产蛋, 不同时间段, 鸡都有产蛋量
  • 牛的产奶量, 不同的测定日, 产奶量不同
  • 猪的饲料消耗, 也是重复测量的数据

只有这样的数据才可以将永久环境效应剖分出来.

重复力是遗传力的上限:
教科书上这样说, 这句话怎么理解呢?
首先, 我们认为
P = G + E P = G + E P=G+E
遗传力为
h 2 = V g / ( V g + V e ) h^2 = Vg/(Vg+Ve) h2=Vg/(Vg+Ve)
这里的Vg如果有重复测量的数据, 可以剖分为可以遗传的部分, 和不可以遗传的部分(永久环境效应), 那么遗传力的计算公式为:
h 2 = V g V g + V p e + V e h^2 = \frac{Vg}{Vg+Vpe+Ve} h2=Vg+Vpe+VeVg​
重复力的计算公式为:
h p e 2 = V g + V p e V g + V p e + V e h_{pe}^2 = \frac{Vg+Vpe}{Vg+Vpe+Ve} hpe2​=Vg+Vpe+VeVg+Vpe​
当Vpe为0时, 重复力=遗传力, 当Vpe>0时, 重复力>遗传力, 所以说重复力是遗传力的上限!

数据使用learnasreml包中的数据

learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的repeatmodel.dat和repeatmodel.ped的数据.

如果没有软件包, 首先安装:

setwd("d:/dmu-test/")
library(devtools)
# install_github("dengfei2013/learnasreml")
library(learnasreml)
data("repeatmodel.dat")
data("repeatmodel.ped")

dat = repeatmodel.dat
ped = repeatmodel.ped

summary(dat)
summary(ped)

看一下数据:

> summary(dat)
     ANIMAL         BYEAR      AGE          YEAR         LAYDATE     
 1      :   5   1000   : 109   2:308   1004   :  79   Min.   : 0.00  
 3      :   5   1001   :  98   3:322   1005   :  78   1st Qu.:20.00  
 9      :   5   999    :  86   4:339   1003   :  69   Median :24.00  
 17     :   5   1002   :  85   5:315   1006   :  64   Mean   :23.54  
 42     :   5   987    :  70   6:323   1002   :  60   3rd Qu.:27.00  
 50     :   5   989    :  66           988    :  54   Max.   :41.00  
 (Other):1577   (Other):1093           (Other):1203                  
> summary(ped)
       ID           FATHER           MOTHER      
 Min.   :   1   Min.   :   0.0   Min.   :   0.0  
 1st Qu.: 328   1st Qu.:   0.0   1st Qu.: 135.0  
 Median : 655   Median :   0.0   Median : 538.0  
 Mean   : 655   Mean   : 261.5   Mean   : 547.4  
 3rd Qu.: 982   3rd Qu.: 458.0   3rd Qu.: 932.0  
 Max.   :1309   Max.   :1304.0   Max.   :1306.0  

数据介绍:

  • 有因子4个: 分别是 ANIMAL BYEAR AGE YEAR
  • 有变量1个: LAYDATE
  • 缺失值用0表示
  • 系谱中有三列数据, 无出生时间一列, 缺失值为0

需要做的处理

  • 系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化
  • 在保存数据时, 去掉行头
  • 编辑DIR文件

使用R语言清洗数据, 并保存数据到D盘dmu-test

dat = repeatmodel.dat
ped = repeatmodel.ped

summary(dat)
summary(ped)

dmuped = ped
dmuped$Birth = 2018

head(dat)
library(data.table)
# write.table(dat,"animal-model.txt",row.names = F,col.names = F)
fwrite(dat,"repeat-model.txt",sep = " ",col.names = F)
fwrite(dmuped,"repeat-ped.txt",sep = " ",col.names = F)

编写DIR文件

想要分析的模型:
观测值: LAYDATE (第四列)
固定因子: 第二列, 第三列, 第四列
随机因子: ID, ide(ID)

所以这里编写DIR
第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

第二部分是分析方法, 默认是AI

$ANALYSE 1 1 0 0

第三部分是定义因子数和变量书, 以及文件位置:

$DATA  ASCII (4,1,0) d:/dmu-test/repeat-model.txt

第四部分, 定义变量名, 也是为了方便结果输出, 相当于数据的行头名

$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT

第五部分, 有6行, 定义模型
整体来说是:
第一行: 单性状 # 1
第二行: 无吸收 # 0
第三行: 主要定义y变量, 固定因子, 随机因子

  • 分析的是第一个变量 # 1
  • 无权重考虑 # 0
  • 共五个因子(固定+随机, 固定写前面, 随机写后面) # 5
  • 第一个固定因子是第二列因子 #2 #BYEAR
  • 第二个固定一致是第三列因子 #3 #AGE
  • 第三个固定因子是第四列 #4 #YEAR
  • 第四个随机因子是第一列 #1 #ANIMAL
  • 第五个随机因子是第一列 #1 #ANIMAL
  • 所以, 5个因子, 三个固定因子:2,3,4, 两个随机因子:1,1 #1 0 5 2 3 4 1 1

第四行: 有两个随机因子, 他们的关系是独立的, 所以是2 1 2

1
0
1 0 5 2 3 4 1 1
2 1 2
0 
0

第六部分: 指定系谱

$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt

注意, 如果想要输出BLUP值, 定义:$DMUAI

$DMUAI
10
1D-7
1D-6
1

完整DIR文件, 放入model.txt中, 然后重命名为: Uni_repeatmodel.DIR

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

$ANALYSE 1 1 0 0
$DATA  ASCII (4,1,0) d:/dmu-test/repeat-model.txt
$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT
$MODEL
1
0
1 0 5 2 3 4 1 1
2 1 2
0 
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt


$DMUAI
10
1D-7
1D-6
1

执行DIR文件

这里运行的run_dmuai.bat, 将DMU安装路径下的文件run_dmuai.bat拷贝到d:/dmu-test文件夹, 在终端cmd界面键入:

run_dmuai.bat Uni_repeatmodel

执行结果:

D:\dmu-test>run_dmuai.bat Uni_repeatmodel

D:\dmu-test>Echo OFF
Starting DMU using Uni_repeatmodel.DIR as directive file

查看结果

在文件*lst中有估算的方差组分, 结果如下:

                             SUMMARY OF MINIMIZATION PROCESS
                                                            
 Eval Criterion    !!Delta!!   !!Gradient!!                   Parameters
 ---- ---------     ---------  ------------ |------------------------------------------
   1   12629.2     0.8574        4.330      |    1.8100        1.8898        1.8705    
  
   2   8234.59      1.370        6.822      |    2.9917        3.3812        3.2879    
  
   3   6444.28      1.776        8.529      |    4.2397        5.4642        5.1761    
  
   4   5857.47      1.566        6.869      |    4.9013        7.4662        6.8832    
  
   5   5736.36     0.6798        2.497      |    4.9407        8.3737        7.6324    
  
   6   5727.01     0.7387E-01   0.2311      |    4.9325        8.4634        7.7233    
  
   7   5726.91     0.1399E-02   0.1596E-02  |    4.9341        8.4621        7.7245    
  
   8   5726.91     0.7706E-04   0.5119E-04  |    4.9340        8.4622        7.7245    
  
   9   5726.91     0.4564E-05   0.2610E-05  |    4.9340        8.4622        7.7245    
  
  10   5726.91     0.2695E-06   0.1558E-06  |    4.9340        8.4622        7.72

可以看到模型收敛

方差组分为:

                   Estimated (co)-variance components
                          ----------------------------------
                                                            
             Parameter vector for L at convergence          
             Asymptotic SE based on AI-information matrix   
                                                            
               No          Parameter             Asymp. S.E.
                                                            
                1           4.93404               1.76364    
                2           8.46217               1.63818    
                3           7.72445              0.329943    

遗传力需要手动计算, 这里还没有找到解决方案.

对比asreml的结果:

代码:

library(asreml)

head(dat)
dat[dat$LAYDATE==0,]$LAYDATE=NA
ainv = asreml.Ainverse(ped)$ginv
mod = asreml(LAYDATE ~ BYEAR + AGE + YEAR, random = ~ ped(ANIMAL)+ide(ANIMAL), ginverse = list(ANIMAL=ainv),data=dat)
summary(mod)$varcomp
pin(mod,h2 ~ V1/(V1+V2+V3))

方差组分:

> summary(mod)$varcomp
                    gamma component std.error   z.ratio constraint
ped(ANIMAL)!ped 0.6387559  4.934041 1.7636385  2.797649   Positive
ide(ANIMAL)!id  1.0955038  8.462169 1.6381812  5.165588   Positive
R!variance      1.0000000  7.724454 0.3299432 23.411466   Positive

遗传力:

> pin(mod,h2 ~ V1/(V1+V2+V3))
   Estimate         SE
h2 0.233612 0.07907261

DMU和asreml比较

两者方差组分一致.

相关阅读:

DMU-参数介绍-学习笔记1
DMU-单性状动物模型-学习笔记2
DMU-单性状重复力模型-学习笔记3
DMU-多性状动物模型-学习笔记4
DMU-单性状动物模型-母体效应–学习笔记5
DMU软件 语法高亮 vim设置–学习笔记6

关注我的公众号:R-breeding

上一篇:Javscript 基础


下一篇:JavaSE-11.2.2【抽象类的成员特点】