asreml-r 在运行中不收敛怎么办？

2023-12-13 12:58:21

在asreml-r运行中，经常遇到模型不收敛的情况，这里我们介绍3种方法来解决不收敛的情况。

导入数据和模型

library(asreml) # load the package
data(“harvey”)
head(harvey)
str(harvey)

ped <- harvey[,1:3]
ainv <- asreml.Ainverse(ped)$ginv
head(ainv)

运行单性状模型y2

m1 <- asreml(y2 ~ Line, random = ~ ped(Calf), ginverse =list(Calf=ainv),data=harvey)
summary(m1)$varcomp   # ped is 0, R is 287

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

	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

可以看到，ped的方差组分基本为0，残差R为287

plot(m1)

运行单性状模型y3

m2 <- asreml(y3 ~ Line, random = ~ ped(Calf), ginverse =list(Calf=ainv),data=harvey)
summary(m2)$varcomp # ped is 500, R is 410

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

	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

**可以看到ped为500，残差R为410

plot(m2)

运行双性状模型

1，直接运行

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

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

	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

可以看到，模型收敛

2，增大迭代次数（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，用update函数

m_12two <- update(m_12)
summary(m_12two)$varcomp

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

	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，用init设置初始值

# 注意这里，us(trait,init=c(0,0.1,499))中，0是y2中ped的方差组分，0.1是协方差（未知，这里设置为0.1），499是y3的方差组分，
# 同理rcov是两性状残差的方差组分
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

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

	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

码农公寓