育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

asreml 设定初始值 固定初始值

已有 1016 次阅读 2019-7-23 21:58 |个人分类:数量遗传学|系统分类:科研笔记

1. 背景

一个朋友问我,如何固定asreml的初始值,现在分为单性状和多性状进行说明。

为何要固定初始值:
1,由于群体较小,估算的方差组分不准确,需要手动设定初始值,直接进行求解
2,有些群体数据,估算方差组分不收敛,需要手动固定初始值

为何要设定初始值:
1,从头进行估算,模型运行时间较长,根据先验信息,手动设定初始值,迭代收敛速度更快
2,多性状分析中,模型不容易收敛,手动设定初始值,更容易收敛和迭代

2. 单性状设定初始值和固定初始值

以asreml包中自带的数据harvey为例,进行演示。

> 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

数据前三列为系谱数据,Line为固定因子,ageOfDam为协变量,y1,y2,y3为三个性状。

2.1 运行单性状动物模型

# 计算A逆矩阵
ainv <- asreml.Ainverse(harvey[,1:3])$ginv
head(ainv)

# 1. 单性状模型
mod1 <- asreml(y1 ~ Line,random =~ ped(Calf),ginverse = list(Calf=ainv),data=harvey)
summary(mod1)$varcomp

结果如下:

> summary(mod1)$varcomp
                 gamma component std.error   z.ratio constraint
ped(Calf)!ped 2.144929 108.83588 106.37372 1.0231463   Positive
R!variance    1.000000  50.74101  86.63851 0.5856635   Positive

可以看到Va为108.83,Ve为50.74,模型收敛。

2.2 单性状动物模型设定初始值

设定初始值,是为了更好的收敛,不影响结果。

# 1.1. 单性状设定初始值
mod <- asreml(y1 ~ Line,random =~ ped(Calf),
              ginverse = list(Calf=ainv),
              start.values = T,
              data=harvey)
vc = mod$gammas.table
vc
vc$Value = c(100,50)
vc
mod1.1 <- asreml(y1 ~ Line,random =~ ped(Calf),
              ginverse = list(Calf=ainv),
              G.param = vc,R.param = vc,
              data=harvey)
summary(mod1.1)$varcomp

结果:

> summary(mod1.1)$varcomp
                  gamma component std.error   z.ratio constraint
ped(Calf)!ped 108.83606 108.83606 106.37146 1.0231697   Positive
R!variance      1.00000   1.00000        NA        NA      Fixed
R!units.var    50.74109  50.74109  86.63707 0.5856742   Positive

2.3 单性状动物模型固定初始值

固定初始值,直接求解,asreml的结果方差组分状态为Fixed

# 1.2. 单性状固定方差组分
mod <- asreml(y1 ~ Line,random =~ ped(Calf),
              ginverse = list(Calf=ainv),
              start.values = T,
              data=harvey)
vc = mod$gammas.table
vc
vc$Value = c(100,50)
vc$Constraint = rep("F",2)
vc
mod1.2 <- asreml(y1 ~ Line,random =~ ped(Calf),
                 ginverse = list(Calf=ainv),
                 G.param = vc,R.param = vc,
                 data=harvey)
summary(mod1.2)$varcomp

结果:

> summary(mod1.2)$varcomp
              gamma component std.error z.ratio constraint
ped(Calf)!ped   100       100        NA      NA      Fixed
R!variance       50        50        NA      NA      Fixed

结果可以看出,方差组分变为了100,50,同时状态是Fixed,说明是固定方差组分的结果,这样计算的BLUP值就是我们想要的。

3. 多性状固定方差组分

3.1 运行多性状模型

# 2. 多性状模型
mod2 <- asreml(cbind(y1,y3) ~ trait + trait:Line,
               random =~ us(trait):ped(Calf),
               rcov = ~ (units):us(trait),
               ginverse = list(Calf=ainv),data=harvey)
summary(mod2)$varcomp
> summary(mod2)$varcomp
                                gamma component std.error    z.ratio constraint
trait:ped(Calf)!trait.y1:y1 108.83746 108.83746 106.37437  1.0231549   Positive
trait:ped(Calf)!trait.y3:y1 -51.25056 -51.25056 166.86351 -0.3071406   Positive
trait:ped(Calf)!trait.y3:y3 499.55701 499.55701 500.53419  0.9980477   Positive
R!variance                    1.00000   1.00000        NA         NA      Fixed
R!trait.y1:y1                50.73993  50.73993  86.63929  0.5856457   Positive
R!trait.y3:y1               -21.53905 -21.53905 136.25598 -0.1580778   Positive
R!trait.y3:y3               273.13654 273.13654 410.03528  0.6661294   Positive

3.2 多性状模型固定方差组分

# 2.2 固定初始值
Va = c(108,-51,499)
Ve = c(50,-21,273)

mod2.2 <- asreml(cbind(y1,y3) ~ trait + trait:Line,
                 random =~ us(trait,init=Va):ped(Calf),
                 rcov = ~ units:us(trait,init=Ve),
                 start.values = TRUE,
                 ginverse = list(Calf=ainv),data=harvey)
vc = mod2.2$gammas.table
vc
vc$Value = c(Va,1,Ve)
vc$Constraint = c(rep("F",7))
vc

mod2.3 <- asreml(cbind(y1,y3) ~ trait + trait:Line,
                 random =~ us(trait,init=Va):ped(Calf),
                 rcov = ~ units:us(trait,init=Ve),
                 G.param = vc,R.param = vc,
                 ginverse = list(Calf=ainv),data=harvey)
summary(mod2.3)$varcomp

结果:

> summary(mod2.3)$varcomp
                            gamma component std.error z.ratio constraint
trait:ped(Calf)!trait.y1:y1   108       108        NA      NA      Fixed
trait:ped(Calf)!trait.y3:y1   -51       -51        NA      NA      Fixed
trait:ped(Calf)!trait.y3:y3   499       499        NA      NA      Fixed
R!variance                      1         1        NA      NA      Fixed
R!trait.y1:y1                  50        50        NA      NA      Fixed
R!trait.y3:y1                 -21       -21        NA      NA      Fixed
R!trait.y3:y3                 273       273        NA      NA      Fixed

4. 结论

  • 1,固定方差组分和设置方差组分方法类似, 不同的是constraintFixed
  • 2,设定方差组分时,先要运行start.values=T,这样就可以生产一个表格,进行修改value和contraint即可
  • 3,单性状和多性状设定方法类似

更多知识,欢迎关注我的公众号:R-breeding




http://blog.sciencenet.cn/blog-2577109-1190844.html

上一篇:学习Julia与弯道超车
下一篇:R语言如何绘制小提琴图violin plot

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-3-2 07:56

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部