||
原书ASReml代码:
// asreml style
fm <- asreml( h5 ~ 1 + Rep,
random = ~ Fam ,
data = df,
subset = Spacing==3,
maxit = 30 )
AFEchidna代码:
// AFEchidna style
fm <- echidna( h5 ~ 1 + Rep,
random = ~ Fam ,
es0.file = 'dfm.es0',
qualifier = '!filter Spacing !select 1',
maxit = 30 )
AFEchidna代码结果:
> fm <- echidna( h5 ~ 1 + Rep,
+ random = ~ Fam ,
+ es0.file = 'dfm.es0',
+ qualifier = '!filter Spacing !select 1',
+ maxit = 30 )
Running Echidna for analysis: h5
Thu Aug 12 09:58:26 2021
Iteration LogL eSigma NEDF
1 1 -2667.81 5289 551
2 2 -2667.71 5340 551
3 3 -2667.71 5333 551
Thu Aug 12 09:58:26 2021 LogL Converged
> waldT(fm,term=c('mu','Rep'))
Wald tests for fixed effects:
DF Wald-F Pr(Chisq)
mu 1 19320.6 < 2.2e-16 ***
Rep 4 70.2 2.065e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Var(fm)
Term Sigma SE Z.ratio
1 Residual 5333.10 337.20 15.815836
2 Fam 441.96 187.28 2.359889
> coef(fm)$fixed
Term Level Effect SE
1 Rep 1 0.000 0.000
2 Rep 2 2.023 10.420
3 Rep 3 109.833 10.138
4 Rep 4 93.284 10.114
5 Rep 5 -11.506 10.095
6 mu 1 547.770 8.038
> coef(fm)$random[1:5,]
Term Level Effect SE
1 Fam 70048 14.508 14.236
2 Fam 70017 8.979 15.683
3 Fam 70002 7.075 15.363
4 Fam 70010 -4.220 14.244
5 Fam 70041 -5.099 16.021
> pin(fm, mulp=c(h2~4* V2/(V1+ V2)))
variance components are as following:
Term Sigma SE vcS
1 Residual 5333.10 337.20 V1
2 Fam 441.96 187.28 V2
pin formula:
h2 ~ 4 * V2/(V1 + V2)
Term Estimate SE
1 h2 0.306 0.124
......
原书ASReml代码:
// asreml style
fm2<-asreml(cbind(dj,h5)~ trait+trait:Rep,
random=~ us(trait):Fam,
residual=~ units:us(trait),
data = df,
subset = Spacing==3 )
AFEchidna代码:
// AFEchidna style
fm2<-echidna(cbind(dj,h5)~ Trait+Trait:Rep,
random=~ us(Trait):Fam,
residual=~ units:us(Trait),
qualifier = '!filter Spacing !select 1',
mulT=T,
es0.file='dfm.es0')
AFEchidna代码结果:
> fm2<-echidna(cbind(dj,h5)~ Trait+Trait:Rep,
+ random=~ us(Trait):Fam,
+ residual=~ units:us(Trait),
+ qualifier = '!filter Spacing !select 1',
+ mulT=T,
+ es0.file='dfm.es0')
Running Echidna for analysis: dj h5
Thu Aug 12 10:04:43 2021
Iteration LogL NEDF
1 1 -5806.49 1105
2 2 -5724.48 1105
3 3 -5676.36 1105
4 4 -5682.81 1105
5 5 -5683.72 1105
6 6 -5683.74 1105
7 7 -5683.74 1105
Thu Aug 12 10:04:44 2021 LogL Converged
> Var(fm2)
Term Sigma SE Z.ratio
1 Residual 1.000 0.000 Inf
2 us(Trait):Fam;us(Trait) 63.930 22.095 2.893415
3 us(Trait):Fam;us(Trait) 75.472 46.458 1.624521
4 us(Trait):Fam;us(Trait) 456.040 189.740 2.403500
5 units:us(Trait);us(Trait) 497.880 31.424 15.843941
6 units:us(Trait);us(Trait) -212.900 73.616 -2.892034
7 units:us(Trait);us(Trait) 5325.500 336.490 15.826622
> waldT(fm2,term=c('Trait','Trait.Rep'))
Wald tests for fixed effects:
DF Wald-F Pr(Chisq)
Trait 2 35334 < 2.2e-16 ***
Trait.Rep 8 47 1.828e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> head(coef(fm2)$random)
Term Level Effect SE
1 us(Trait):Fam 1.70048 5.164 4.602
2 us(Trait):Fam 1.70017 -0.462 5.162
3 us(Trait):Fam 1.70002 -1.895 5.034
4 us(Trait):Fam 1.7001 -8.938 4.605
5 us(Trait):Fam 1.70041 3.865 5.301
6 us(Trait):Fam 1.70031 2.450 5.030
> pin(fm2,mulp=c(h2_dj~4*V2/(V2+V5),
+ h2_h5~4*V4/(V4+V7)))
variance components are as following:
Term Sigma SE vcS
1 Residual 1.000 0.000 V1
2 us(Trait):Fam;us(Trait) 63.930 22.095 V2
3 us(Trait):Fam;us(Trait) 75.472 46.458 V3
4 us(Trait):Fam;us(Trait) 456.040 189.740 V4
5 units:us(Trait);us(Trait) 497.880 31.424 V5
6 units:us(Trait);us(Trait) -212.900 73.616 V6
7 units:us(Trait);us(Trait) 5325.500 336.490 V7
pin formula:
h2_dj ~ 4 * V2/(V2 + V5)
h2_h5 ~ 4 * V4/(V4 + V7)
Term Estimate SE
1 h2_dj 0.455 0.145
11 h2_h5 0.316 0.125
> pin(fm2,mulp=c(gcorr~V3/sqrt(V2*V4),
+ pcorr~(V3+V6)/sqrt((V2+V5)*(V4+V7))),signif=T)
variance components are as following:
Term Sigma SE vcS
1 Residual 1.000 0.000 V1
2 us(Trait):Fam;us(Trait) 63.930 22.095 V2
3 us(Trait):Fam;us(Trait) 75.472 46.458 V3
4 us(Trait):Fam;us(Trait) 456.040 189.740 V4
5 units:us(Trait);us(Trait) 497.880 31.424 V5
6 units:us(Trait);us(Trait) -212.900 73.616 V6
7 units:us(Trait);us(Trait) 5325.500 336.490 V7
pin formula:
gcorr ~ V3/sqrt(V2 * V4)
pcorr ~ (V3 + V6)/sqrt((V2 + V5) * (V4 + V7))
Term Estimate SE Siglevel
1 gcorr 0.4420 0.2570 *
11 pcorr -0.0763 0.0449 *
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
AFEchidna包免费对外开放,仅用于学术研究,不可用于商业研究,否则造成的后果自负。
感兴趣者,可发邮件到yzhlinscau@163.com免费索取程序包。
参考文献:
Zhang WH, Wei RY, Liu Y, Lin YZ. AFEchidna is a R package for genetic evaluation of plant and animal breeding datasets. BioRxiv. DOI: 10.1101/2021.06.24.449740.
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-20 11:24
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社