在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 

可以看到，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 

可以看到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 

可以看到，模型收敛

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 

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