当前位置：网站首页>What if asreml-r does not converge in operation?
What if asreml-r does not converge in operation?
2022-06-27 02:06:00 【Analysis of breeding data】
stay asreml-r Running , It is often encountered that the model does not converge , Here we introduce 3 There are two ways to solve the problem of non convergence .
Import data and models
library(asreml) # load the package
data(“harvey”)
head(harvey)
str(harvey)
ped <- harvey[,1:3]
ainv <- asreml.Ainverse(ped)$ginv
head(ainv)
Running a single character model y2
m1 <- asreml(y2 ~ Line, random = ~ ped(Calf), ginverse =list(Calf=ainv),data=harvey)
summary(m1)$varcomp # ped is 0, R is 287
- 1.
- 2.
ASReml: Thu Apr 06 17:23:38 2017
LogLik S2 DF wall cpu
-264.4319 1474.0895 62 17:23:38 0.0 (1 restrained)
-264.3059 1592.5437 62 17:23:38 0.0 (1 restrained)
-264.2982 1600.8401 62 17:23:38 0.0 (1 restrained)
-264.2977 1601.3686 62 17:23:38 0.0 (1 restrained)
-264.2977 1601.4020 62 17:23:38 0.0
-264.2977 1601.4020 62 17:23:38 0.0
-264.2977 1601.4020 62 17:23:38 0.0
Finished on: Thu Apr 06 17:23:38 2017
LogLikelihood Converged
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
gamma | component | std.error | z.ratio | constraint | |
ped(Calf)!ped | 1.6e-06 | 2.562243e-03 | 4.601925e-04 | 5.567764 | Boundary |
R!variance | 1.0e+00 | 1.601402e+03 | 2.876203e+02 | 5.567764 | Positive |
You can see ,ped The components of variance are basically 0, residual R by 287
plot(m1)
- 1.
Running a single character model y3
m2 <- asreml(y3 ~ Line, random = ~ ped(Calf), ginverse =list(Calf=ainv),data=harvey)
summary(m2)$varcomp # ped is 500, R is 410
- 1.
- 2.
ASReml: Thu Apr 06 17:22:03 2017
LogLik S2 DF wall cpu
-239.6964 663.7311 62 17:22:03 0.0
-239.3469 591.8505 62 17:22:03 0.0
-239.0326 489.6240 62 17:22:03 0.0
-238.8604 365.1480 62 17:22:03 0.0
-238.8324 297.4237 62 17:22:03 0.0
-238.8306 275.3345 62 17:22:03 0.0
-238.8305 273.1609 62 17:22:03 0.0
-238.8305 273.1669 62 17:22:03 0.0
Finished on: Thu Apr 06 17:22:03 2017
LogLikelihood Converged
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
- 15.
gamma | component | std.error | z.ratio | constraint | |
ped(Calf)!ped | 1.828629 | 499.5207 | 500.5071 | 0.9980292 | Positive |
R!variance | 1.000000 | 273.1669 | 410.0170 | 0.6662330 | Positive |
** You can see ped by 500, residual R by 410
plot(m2)
- 1.
Run the androgynous model
1, Direct operation
m_12 <- asreml(cbind(y2,y3)~trait, random = ~ us(trait):ped(Calf),
ginverse = list(Calf=ainv),
rcov = ~ units:us(trait),data=harvey)
summary(m_12)$varcomp
- 1.
- 2.
- 3.
- 4.
ASReml: Thu Apr 06 17:22:07 2017
LogLik S2 DF wall cpu
-597.9669 1.0000 128 17:22:07 0.0 (3 restrained)
-571.8584 1.0000 128 17:22:07 0.0 (3 restrained)
-543.2659 1.0000 128 17:22:07 0.0
-522.7403 1.0000 128 17:22:07 0.0
-517.3376 1.0000 128 17:22:07 0.0
-516.5146 1.0000 128 17:22:07 0.0
-516.4598 1.0000 128 17:22:07 0.0
-516.4590 1.0000 128 17:22:07 0.0
-516.4590 1.0000 128 17:22:07 0.0
Finished on: Thu Apr 06 17:22:07 2017
LogLikelihood Converged
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
- 15.
- 16.
gamma | component | std.error | z.ratio | constraint | |
trait:ped(Calf)!trait.y2:y2 | 335.91633 | 335.91633 | 640.2366 | 0.5246753 | Positive |
trait:ped(Calf)!trait.y3:y2 | -76.79937 | -76.79937 | 383.5906 | -0.2002118 | Positive |
trait:ped(Calf)!trait.y3:y3 | 546.65342 | 546.65342 | 459.3477 | 1.1900645 | Positive |
R!variance | 1.00000 | 1.00000 | NA | NA | Fixed |
R!trait.y2:y2 | 1387.06671 | 1387.06671 | 635.5144 | 2.1825889 | Positive |
R!trait.y3:y2 | 219.37425 | 219.37425 | 343.8650 | 0.6379663 | Positive |
R!trait.y3:y3 | 238.55889 | 238.55889 | 382.3906 | 0.6238618 | Positive |
You can see , Model convergence
2, Increase the number of iterations （maxit=1000）
m_12one <- asreml(cbind(y2,y3)~trait, random = ~ us(trait):ped(Calf),
ginverse = list(Calf=ainv),
rcov = ~ units:us(trait),data=harvey,
maxit = 1000)
summary(m_12one)$varcomp
3, use update function
m_12two <- update(m_12)
summary(m_12two)$varcomp
- 1.
- 2.
ASReml: Thu Apr 06 17:22:12 2017
LogLik S2 DF wall cpu
-516.4590 1.0000 128 17:22:12 0.0
-516.4590 1.0000 128 17:22:12 0.0
-516.4590 1.0000 128 17:22:12 0.0
-516.4590 1.0000 128 17:22:12 0.0
Finished on: Thu Apr 06 17:22:12 2017
LogLikelihood Converged
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
gamma | component | std.error | z.ratio | constraint | |
trait:ped(Calf)!trait.y2:y2 | 335.91632 | 335.91632 | 640.2367 | 0.5246752 | Positive |
trait:ped(Calf)!trait.y3:y2 | -76.79938 | -76.79938 | 383.5904 | -0.2002119 | Positive |
trait:ped(Calf)!trait.y3:y3 | 546.65342 | 546.65342 | 459.3473 | 1.1900656 | Positive |
R!variance | 1.00000 | 1.00000 | NA | NA | Fixed |
R!trait.y2:y2 | 1387.06671 | 1387.06671 | 635.5145 | 2.1825887 | Positive |
R!trait.y3:y2 | 219.37425 | 219.37425 | 343.8649 | 0.6379664 | Positive |
R!trait.y3:y3 | 238.55889 | 238.55889 | 382.3903 | 0.6238623 | Positive |
4, use init Set initial value
# Note that there ,us(trait,init=c(0,0.1,499)) in ,0 yes y2 in ped The variance components of ,0.1 It's covariance （ Unknown , I'm going to set it to 0.1）,499 yes y3 The variance components of ,
# Empathy rcov It's the variance component of the bisexual residuals
m_12 <- asreml(cbind(y2,y3)~trait, random = ~ us(trait,init=c(0,0.1,499)):ped(Calf),
ginverse = list(Calf=ainv),
rcov =~ units:us(trait,init=c(1601,0.1,273)),data=harvey)
summary(m_12)$varcomp
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
ASReml: Thu Apr 06 17:22:13 2017
LogLik S2 DF wall cpu
-520.4297 1.0000 128 17:22:13 0.0 (3 restrained)
-518.3307 1.0000 128 17:22:13 0.0 (2 restrained)
-517.1351 1.0000 128 17:22:13 0.0
-516.4907 1.0000 128 17:22:13 0.0
-516.4591 1.0000 128 17:22:13 0.0
-516.4590 1.0000 128 17:22:13 0.0
-516.4590 1.0000 128 17:22:13 0.0
Finished on: Thu Apr 06 17:22:13 2017
LogLikelihood Converged
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
gamma | component | std.error | z.ratio | constraint | |
trait:ped(Calf)!trait.y2:y2 | 335.91632 | 335.91632 | 640.2368 | 0.5246751 | Positive |
trait:ped(Calf)!trait.y3:y2 | -76.79937 | -76.79937 | 383.5905 | -0.2002119 | Positive |
trait:ped(Calf)!trait.y3:y3 | 546.65342 | 546.65342 | 459.3474 | 1.1900655 | Positive |
R!variance | 1.00000 | 1.00000 | NA | NA | Fixed |
R!trait.y2:y2 | 1387.06672 | 1387.06672 | 635.5146 | 2.1825884 | Positive |
R!trait.y3:y2 | 219.37425 | 219.37425 | 343.8649 | 0.6379663 | Positive |
R!trait.y3:y3 | 238.55889 | 238.55889 | 382.3903 | 0.6238622 | Positive |
边栏推荐
- Reading a book in idea is too much!
- 1.44 inch TFT-LCD display screen mold taking tutorial
- “所有专业都在劝退”，对大学生最友好的竟然是它？
- LeetCode 785：判断二分图
- [the path of system analyst] Chapter 6: duplicate demand engineering (case paper)
- Don't be brainwashed. This is the truth about the wages of 90% of Chinese people
- NOKOV动作捕捉系统使多场协同无人机自主建造成为可能
- UVM in UVM_ config_ Use of DB in sequence
- Due to the invalidation of the prospectus of bori technology, CICC has stopped providing guidance to it and abandoned the listing on the Hong Kong stock exchange?
- Oracle/PLSQL: Rpad Function
猜你喜欢
随机推荐
按键控制LED状态翻转
dat.gui.js星星圆圈轨迹动画js特效
three.js多米诺骨牌js特效
Nokov motion capture system makes it possible for multi field cooperative UAV to build independently
ConstraintLayout（约束布局）开发指南
使用命令行安装达梦数据库
idea 插件开发一些异常处理
CVPR2022 | PointDistiller:面向高效紧凑3D检测的结构化知识蒸馏
福元医药上市在即：募资净额将达到16亿元，胡柏藩为实际控制人
memcached基础10
Consumers pursue the iPhone because its cost performance exceeds that of domestic mobile phones
Look! In June, 2022, the programming language ranking list was released! The first place is awesome
两个页面之间传参方法
snakemake 使用的注意事项
Recursion will make strtok more attractive
Oracle/PLSQL: Lpad Function
Parameter estimation -- Chapter 7 study report of probability theory and mathematical statistics (point estimation)
Oracle/PLSQL: Cast Function
memcached基础14
p5.js死亡星球