asreml怎样设定初始值,很多新手对此不是很清楚,为了帮助大家解决这个难题,下面小编将为大家详细讲解,有这方面需求的人可以来学习下,希望你能有所收获。
一个朋友问我,如何固定asreml的初始值,现在分为单性状和多性状进行说明。
为何要固定初始值:
1,由于群体较小,估算的方差组分不准确,需要手动设定初始值,直接进行求解
2,有些群体数据,估算方差组分不收敛,需要手动固定初始值
为何要设定初始值:
1,从头进行估算,模型运行时间较长,根据先验信息,手动设定初始值,迭代收敛速度更快
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为三个性状。
# 计算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,模型收敛。
设定初始值,是为了更好的收敛,不影响结果。
# 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
固定初始值,直接求解,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值就是我们想要的。
# 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
# 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
1,固定方差组分和设置方差组分方法类似, 不同的是constraint
为Fixed
2,设定方差组分时,先要运行start.values=T
,这样就可以生产一个表格,进行修改value和contraint即可
3,单性状和多性状设定方法类似
看完上述内容是否对您有帮助呢?如果还想对相关知识有进一步的了解或阅读更多相关文章,请关注亿速云行业资讯频道,感谢您对亿速云的支持。
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。