以MCMC分析双性状为例：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 | library(AAfun) library(MCMCglmm) data(PrSpa) df<-subset(PrSpa,Spacing=='3') ######## bi-trait model ######## df$dj<-1000*df$dj phen.var<-matrix(c(var(df$dj,na.rm=TRUE),0,0, var(df$h5,na.rm=TRUE)),2,2) prior<-list(G=list(G1=list(V=phen.var/2,n=2)), R=list(V=phen.var/2,n=2)) set.seed(1234) m2.glmm <- MCMCglmm(cbind(dj,h5)~ trait-1+trait:Rep, random=~ us(trait):Fam, rcov=~ us(trait):units, data=df, family=c("gaussian", "gaussian"), nitt=130000,thin=100,burnin=30000, prior=prior,verbose=F,pr=T) posterior.mode(m2.glmm$VCV) HPDinterval(m2.glmm$VCV) #### count se for variance component mc.se(m2.glmm$VCV) #### count se for heritability A.h2.glmm<-4*m2.glmm$VCV[,'dj:dj.Fam']/(m2.glmm$VCV[,'dj:dj.Fam']+m2.glmm$VCV[,'dj:dj.units']) mc.se(A.h2.glmm) #### count se for corr gCorr.glmm<-m2.glmm$VCV[,'h5:dj.Fam']/sqrt(m2.glmm$VCV[,'dj:dj.Fam']*m2.glmm$VCV[,'h5:h5.Fam']) mc.se(gCorr.glmm,sigf=1) |
运行的结果输出为：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | ########## MCMC的默认输出结果 ########## > posterior.mode(m2.glmm$VCV) #方差分量 dj:dj.Fam h5:dj.Fam dj:h5.Fam h5:h5.Fam dj:dj.units 87.5 53.9 53.9 797.2 492.7 h5:dj.units dj:h5.units h5:h5.units -175.0 -175.0 5434.9 ########## mc.se()调整后的输出结果 ########## > mc.se(m2.glmm$VCV) var se z.ratio dj:dj.Fam 87.545 24.344 3.596 h5:dj.Fam 53.944 54.253 0.994 dj:h5.Fam 53.944 54.253 0.994 h5:h5.Fam 797.224 218.427 3.650 dj:dj.units 492.689 29.299 16.816 h5:dj.units -174.974 75.992 -2.303 dj:h5.units -174.974 75.992 -2.303 h5:h5.units 5434.884 343.162 15.838 ######### 调整前的遗传力 ######### > posterior.mode(A.h2.glmm) var1 0.596 ######### 调整后的遗传力 ######### > mc.se(A.h2.glmm) var se z.ratio [1,] 0.596 0.149 4 ######### 调整前的遗传相关 ######### > posterior.mode(gCorr.glmm) var1 0.201 ######### 调整后的遗传相关 ######### > mc.se(gCorr.glmm,sigf=1) var se z.ratio sig.level [1,] 0.201 0.169 1.19 Not signif --------------- Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1 |
