asreml 设定初始值 固定初始值
2019-07-23 本文已影响1人
育种数据分析之放飞自我
. 背景
一个朋友问我，如何固定asreml的初始值，现在分为单性状和多性状进行说明。
为何要固定初始值：
1，由于群体较小，估算的方差组分不准确，需要手动设定初始值，直接进行求解
2，有些群体数据，估算方差组分不收敛，需要手动固定初始值
为何要设定初始值：
1，从头进行估算，模型运行时间较长，根据先验信息，手动设定初始值，迭代收敛速度更快
2，多性状分析中，模型不容易收敛，手动设定初始值，更容易收敛和迭代
2. 单性状设定初始值和固定初始值
以asreml包中自带的数据harvey为例，进行演示。
> library(asreml)> data(harvey)> head(harvey) Calf Sire Dam Line ageOfDam y1 y2 y31 101 Sire_1 0 1 3 192 390 2242 102 Sire_1 0 1 3 154 403 2653 103 Sire_1 0 1 4 185 432 2414 104 Sire_1 0 1 4 183 457 2255 105 Sire_1 0 1 5 186 483 2586 106 Sire_1 0 1 5 177 469 267数据前三列为系谱数据，Line为固定因子，ageOfDam为协变量，y1,y2,y3为三个性状。
2.1 运行单性状动物模型
# 计算A逆矩阵ainv <- asreml.Ainverse(harvey[,1:3])$ginvhead(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 constraintped(Calf)!ped 2.144929 108.83588 106.37372 1.0231463 PositiveR!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.tablevcvc$Value = c(100,50)vcmod1.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 constraintped(Calf)!ped 108.83606 108.83606 106.37146 1.0231697 PositiveR!variance 1.00000 1.00000 NA NA FixedR!units.var 50.74109 50.74109 86.63707 0.5856742 Positive2.3 单性状动物模型固定初始值
固定初始值，直接求解，asreml的结果方差组分状态为Fixed
# 1.2. 单性状固定方差组分mod <- asreml(y1 ~ Line,random =~ ped(Calf), ginverse = list(Calf=ainv), start.values = T, data=harvey)vc = mod$gammas.tablevcvc$Value = c(100,50)vc$Constraint = rep("F",2)vcmod1.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 constraintped(Calf)!ped 100 100 NA NA FixedR!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 constrainttrait:ped(Calf)!trait.y1:y1 108.83746 108.83746 106.37437 1.0231549 Positivetrait:ped(Calf)!trait.y3:y1 -51.25056 -51.25056 166.86351 -0.3071406 Positivetrait:ped(Calf)!trait.y3:y3 499.55701 499.55701 500.53419 0.9980477 PositiveR!variance 1.00000 1.00000 NA NA FixedR!trait.y1:y1 50.73993 50.73993 86.63929 0.5856457 PositiveR!trait.y3:y1 -21.53905 -21.53905 136.25598 -0.1580778 PositiveR!trait.y3:y3 273.13654 273.13654 410.03528 0.6661294 Positive3.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.tablevcvc$Value = c(Va,1,Ve)vc$Constraint = c(rep("F",7))vcmod2.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 constrainttrait:ped(Calf)!trait.y1:y1 108 108 NA NA Fixedtrait:ped(Calf)!trait.y3:y1 -51 -51 NA NA Fixedtrait:ped(Calf)!trait.y3:y3 499 499 NA NA FixedR!variance 1 1 NA NA FixedR!trait.y1:y1 50 50 NA NA FixedR!trait.y3:y1 -21 -21 NA NA FixedR!trait.y3:y3 273 273 NA NA Fixed4. 结论
- 1，固定方差组分和设置方差组分方法类似， 不同的是
constraint为Fixed - 2，设定方差组分时，先要运行
start.values=T，这样就可以生产一个表格，进行修改value和contraint即可 - 3，单性状和多性状设定方法类似
更多知识，欢迎关注我的公众号：R-breeding
