asreml-r 模型不收敛怎么办

1、导入数据和模型

library(asreml) # load the packagedata(“harvey”)head(harvey)str(harvey)ped <- harvey[,1:3]ainv <- asreml.Ainverse(ped)$ginvhead(ainv)

2、运行单性状模型y2

m1 <- asreml(y2 ~ Line, random = ~ ped(Calf), ginverse =list(Calf=ainv),data=harvey) summary(m1)$varcomp   > summary(m1)$varcomp                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

3、运行单性状模型y3

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

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

1、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

可以看到，模型收敛

2、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、3，用update函数

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

4、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