AFEchidna包主函数echidna更新

AFEchidna包主函数echidna()用法(新版本)如下：

echidna(es0.file,trait,fixed,random,residual,
                   delf,foldN,mulT,mulN,
                   met,cycle,
                   trace, maxit,
                   batch,batch.G,batch.R,
                   predict,vpredict,
                   qualifier,jobqualf)

增加了参数batch，将之前的echidna.batch()函数功能纳入；增加了参数batch.G和batch.R，用于一次性运行多个G结构或R结构。目前，这三个参数只能单独运行。

简单理解：batch针对测量值Y；batch.G针对随机效应；batch.R针对误差效应。

1. 参数batch示范–针对测量值Y

// old style
library(AFEchidna)
path<-"C:/Users/yzhlinscau/Desktop/echi/exam"
setwd(path)

trait<-c('h3','h4','h5')

res21<-echidna(es0.file='fm.es0',trait=trait,
                     fixed=y~1+Rep,random=~Fam,
                     predict='Fam',batch=T,
                     trace=T)       

运行过程如下：

> res21<-echidna(es0.file='fm.es0',trait=trait,
+                      fixed=y~1+Rep,random=~Fam,
+                      predict='Fam',batch=T,
+                      trace=T)

Program starts running batch analysis ------ 

run h3 -- -- --:
Running Echidna for analysis:  h3 

 Sun Feb  7 15:36:45 2021 
  Iteration     LogL eSigma NEDF
1         1 -2346.93   1576  554
2         2 -2346.84   1590  554
3         3 -2346.84   1589  554
4         4 -2346.84   1589  554
Sun Feb  7 15:36:46 2021 LogL Converged 

run h4 -- -- --:
Running Echidna for analysis:  h4 

 Sun Feb  7 15:36:46 2021 
  Iteration     LogL eSigma NEDF
1         1 -2557.35   3542  551
2         2 -2556.98   3616  551
3         3 -2556.95   3599  551
4         4 -2556.95   3599  551
5         5 -2556.95   3599  551
Sun Feb  7 15:36:46 2021 LogL Converged 

run h5 -- -- --:
Running Echidna for analysis:  h5 

 Sun Feb  7 15:36:47 2021 
  Iteration     LogL eSigma NEDF
1         1 -2667.81   5289  551
2         2 -2667.71   5340  551
3         3 -2667.71   5333  551
4         4 -2667.71   5333  551
Sun Feb  7 15:36:47 2021 LogL Converged 

结果查看：

> names(res21)
[1] "res.all" "org.par"
> Var(res21)

V1-Residual; V2-Fam
Converge: 1 means True; 0 means FALSE.

       V1     V2  V1.se   V2.se Converge maxit
h3 1589.0 132.59 100.29  56.459        1     4
h4 3599.1 241.39 227.54 115.510        1     5
h5 5333.4 442.00 337.21 187.170        1     4
> pin(res21,mulp=c(h2~V2*4/(V2+V1),
+                  Vp~V2+V1),all=T)

results as following: 

pin formula: 
h2 ~ V2 * 4/(V2 + V1)
Vp ~ V2 + V1

terms:  V1-Residual; V2-Fam; V3-h2; V4-Vp

       V1     V2    V3       V4    SE1     SE2   SE3     SE4
h3 1589.0 132.59 0.308 1721.577 100.29  56.459 0.125 106.461
h4 3599.1 241.39 0.251 3840.461 227.54 115.510 0.116 235.833
h5 5333.4 442.00 0.306 5775.345 337.21 187.170 0.124 357.704

当然，其它结果的提取同之前函数，例如pin, predict，coef等。

2. 参数batch.G示范–针对随机效应

r2g<-echidna(es0.file="fm.es0",
             fixed=h5~1+Rep,
             random=c(G1~Fam,G2~Fam+Plot),
             residual=~units,
             batch.G=T,
             trace=T)
             

运行过程如下：

> r2g<-echidna(es0.file="fm.es0",
+              fixed=h5~1+Rep,
+              random=c(G1~Fam,G2~Fam+Plot),
+              residual=~units,
+              batch.G=T,
+              trace=T)

Program runs for 2 more G-structure at one time. ------ 

run 1 -- random effects: Fam
Running Echidna for analysis:  h5 

 Sun Feb  7 15:45:15 2021 
  Iteration     LogL eSigma NEDF
1         1 -2667.81   5289  551
2         2 -2667.71   5340  551
3         3 -2667.71   5333  551
Sun Feb  7 15:45:15 2021 LogL Converged 


run 2 -- random effects: Fam   Plot
Running Echidna for analysis:  h5 

 Sun Feb  7 15:45:16 2021 
  Iteration     LogL eSigma NEDF
1         1 -2671.03   5274  551
2         2 -2668.12   5298  551
3         3 -2668.09   5329  551
4         4 -2667.73   5332  551
5         5 -2667.73   5332  551
Sun Feb  7 15:45:16 2021 LogL Converged 

结果查看：

> nr2g<-b2s(r2g)
> m1=nr2g$G1
> m2=nr2g$G2
> Var(m1)
      Term   Sigma     SE   Z.ratio
1 Residual 5333.10 337.20 15.815836
2      Fam  441.96 187.28  2.359889
> Var(m2)
      Term      Sigma         SE   Z.ratio
1 Residual 5.3320e+03 3.3694e+02 15.824776
2     Plot 5.0493e-02 3.1907e-03 15.825054
3      Fam 4.4301e+02 1.8661e+02  2.373989
> model.comp(m1,m2,LRT=T)
Model comparison results as following:

   parNO     LogL     AIC     BIC AIC.State BIC.State
m1     2 -2667.71 5339.42 5348.05    better    better
m2     3 -2667.73 5341.46 5354.39                    

=====================================
Likelihood ratio test (LRT) results:
note:left model before "/" is full model,right is reduced.

Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

      df LR-statistic Pr(Chisq)
m2/m1  1        -0.04       0.5

当然也可以进行其它函数操作。

3. 参数batch.R示范

r2r<-echidna(es0.file="MET.es0", 
            fixed=yield~1+Loc,
            random=~Genotype:Loc,
            residual=c(R1~sat(Loc):ar1(Col):ar1(Row),R2~sat(Loc):units), 
            batch.R=T,
            met=T)

nr2r<-b2s(r2r)
m1<-nr2r$R1;m2<-nr2r$R2            

运行过程如下：

> r2r<-echidna(es0.file="MET.es0", 
+             fixed=yield~1+Loc,
+             random=~Genotype:Loc,
+             residual=c(R1~sat(Loc):ar1(Col):ar1(Row),R2~sat(Loc):units), 
+             batch.R=T,
+             met=T)

Program runs for 2 more R-structure at one time. ------ 

run 1 -- residual effects: sat(Loc):ar1(Col):ar1(Row)
Running Echidna for analysis:  yield 

 Sun Feb  7 15:56:47 2021 
  Iteration     LogL NEDF
1         1 -1289.61  630
2         2 -1248.28  630
3         3 -1213.33  630
4         4 -1205.09  630
5         5 -1204.11  630
6         6 -1204.05  630
7         7 -1204.05  630
Sun Feb  7 15:56:47 2021 LogL Converged 


run 2 -- residual effects: sat(Loc):units
Running Echidna for analysis:  yield 

 Sun Feb  7 15:56:48 2021 
  Iteration     LogL NEDF
1         1 -1310.30  630
2         2 -1250.82  630
3         3 -1239.42  630
4         4 -1238.61  630
5         5 -1238.60  630
6         6 -1238.60  630
Sun Feb  7 15:56:48 2021 LogL Converged 

结果查看

> Var(m1)
                                    Term     Sigma       SE     Z.ratio
1                               Residual  1.000000 0.000000         Inf
2                           Genotype.Loc  3.032800 0.754500  4.01961564
3                               ar1(Col)  0.010667 0.112900  0.09448184
4  sat(Loc,1).ar1(Col).ar1(Row);ar1(Row) -0.023540 0.124060 -0.18974690
5  sat(Loc,1).ar1(Col).ar1(Row);ar1(Row) 13.143000 2.091800  6.28310546
6                               ar1(Col)  0.173700 0.109940  1.57995270
7  sat(Loc,2).ar1(Col).ar1(Row);ar1(Row)  0.268190 0.123470  2.17210658
8  sat(Loc,2).ar1(Col).ar1(Row);ar1(Row) 20.948000 3.460700  6.05311064
9                               ar1(Col)  0.558990 0.082029  6.81454120
10 sat(Loc,3).ar1(Col).ar1(Row);ar1(Row) -0.058038 0.137610 -0.42175714
11 sat(Loc,3).ar1(Col).ar1(Row);ar1(Row) 23.010000 4.360300  5.27715983
12                              ar1(Col)  0.501670 0.101850  4.92557683
13 sat(Loc,4).ar1(Col).ar1(Row);ar1(Row)  0.139040 0.117820  1.18010525
14 sat(Loc,4).ar1(Col).ar1(Row);ar1(Row) 28.498000 5.512700  5.16951766
15                              ar1(Col) -0.074562 0.136530 -0.54612173
16 sat(Loc,5).ar1(Col).ar1(Row);ar1(Row)  0.014553 0.127140  0.11446437
17 sat(Loc,5).ar1(Col).ar1(Row);ar1(Row)  9.145600 1.465700  6.23974893
18                              ar1(Col)  0.331200 0.100730  3.28799762
19 sat(Loc,6).ar1(Col).ar1(Row);ar1(Row)  0.122330 0.139670  0.87585022
20 sat(Loc,6).ar1(Col).ar1(Row);ar1(Row)  8.085400 1.390300  5.81557937
> Var(m2)
          Term   Sigma      SE  Z.ratio
1     Residual  1.0000 0.00000      Inf
2 Genotype.Loc  2.8914 0.82762 3.493632
3        units 13.2700 2.11590 6.271563
4        units 20.5640 3.17370 6.479503
5        units 24.0250 3.71870 6.460591
6        units 27.7710 4.19240 6.624129
7        units  9.2094 1.47380 6.248745
8        units  8.0193 1.27040 6.312421

模型比较：

> model.comp(m1,m2,LRT=T)

Model comparison results as forllowing:
   parNO     LogL     AIC     BIC AIC.State BIC.State
m2     7 -1238.60 2491.19 2522.31              better
m1    19 -1204.05 2446.10 2530.57    better          
=====================================
Likelihood ratio test (LRT) results:
note:left model before "/" is full model,right is reduced.

Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

      df LR-statistic Pr(Chisq)    
m1/m2 12         69.1 4.446e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

同样，还可以使用其它函数。

今年将陆续完善和更新AFEchidna包。

AFEchidna包最新版：V0.1.0
更新: 2021-02-07