R语言
asreml-r软件包
导入数据和模型library(asreml) # load the packagedata(“harvey”) head(harvey)str(harvey) ped <- harvey[,1:3] ainv <- asreml.Ainverse(ped)$ginvhead(ainv)
运行单性状模型y2m1 <- 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
运行单性状模型y3m2 <- 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，直接运行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，增大迭代次数（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
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
1，首先运行单性状模型，获得方差组分初始值，然后再运行二性状模型，设置初始值
2，asreml软件可以免费申请试用，详见VSNC官网