-10.4.2.1 单性状

原书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

......

- 10.4.2.2 双性状

原书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.