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

单性状重复力模型

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

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

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

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

重复力是遗传力的上限:
教科书上这样说, 这句话怎么理解呢?
首先, 我们认为
P=G+EP = G + EP=G+E
遗传力为
h2=Vg/(Vg+Ve)h^2 = Vg/(Vg+Ve)h2=Vg/(Vg+Ve)
这里的Vg如果有重复测量的数据, 可以剖分为可以遗传的部分, 和不可以遗传的部分(永久环境效应), 那么遗传力的计算公式为:
h2=VgVg+Vpe+Veh^2 = \frac{Vg}{Vg+Vpe+Ve}h2=Vg+Vpe+VeVg
重复力的计算公式为:
hpe2=Vg+VpeVg+Vpe+Veh_{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.pedsummary(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.pedsummary(dat)
summary(ped)dmuped = ped
dmuped$Birth = 2018head(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_repeatmodelD:\dmu-test>Echo OFF
Starting DMU using Uni_repeatmodel.DIR as directive file

查看结果

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

                             SUMMARY OF MINIMIZATION PROCESSEval 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)$varcompgamma 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比较

两者方差组分一致.