AFEchidna包里的函数echidna()可以达到在R中运行Echidna软件的目的。

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

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

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

新旧版本的区别在于：使用了参数es0.file，替代了原来的es0.path；增加了参数met，用于多点试验数据。

echidna()基于零模板(.es0)，通过添加指定性状trait、固定效应fixed、随机效应random和残差效应residual来达到分析目的。

零模板(.es0)，与.es程序类似，但其少了线性模型及其他特定语句或命令符。零模板(.es0)的简单示例如下：

!WORK 2 # !REN  !ARG 4 !out
TITLE: fm # !DOPART $1
 # "TreeID","Mum","Dad","Spacing","Rep","Plot","dj","dm","wd","h1","h2","h3","h4"," ...
 # "80001","70048",NA,"3","1","1",0.334,0.405,0.358,29,130,239,420,630,1,1,1 ...
TreeID  !I    # 80001
Mum !I       # 70048
Dad *        # *
Spacing !I   # 3
Rep  *       # 1
Plot *       # 1
dj  !* 1000      # 0.334
dm        # 0.405
wd        # 0.358
h1        # 29
h2        # 130
h3        # 239
h4        # 420
h5        # 630
dis        # 1
lt         # 1
str        # 1
# Verify data fields are correctly classified as factors (!A !I !P *) or variates
fm.csv !SKIP 1 !maxit 30 !SLN !filter Spacing !select 1

零模板(.es0)在新版本的echidna()里可以直接生成，但仍需要对数据域进行手动设置。

echidna分析流程：

1. 零模板(.es0)与目标数据集在同个目录，然后将该目录路径通过参数es0.path传递给函数echidna()；
   

2. 指定分析性状(trait)、固定效应(fixed)、随机效应(random)和残差效应(residual）；

3. 指定predict等参数；

完成上述步骤后，即可运行echidna()进行分析。

简单示例如下：

// simple example
library(AFEchidna)
es0.path=r"(E:\课件\新建文件夹)" 
res11<-echidna(trait='h3',fixed='Rep',random='Mum',
             residual=NULL,
             predict=c('Mum'),
            es0.path)

运行过程：

// running procedure
> res11<-echidna(trait='h3',fixed='Rep',random='Mum',
+              residual=NULL,
+              predict=c('Mum'),
+             es0.path)
Read 3 items

 Mon Jan 11 12:44:30 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
Mon Jan 11 12:44:30 2021 LogL Converged

运行结果：

// results
> # show result class
> class(res11)
[1] "esR"
> names(res11) # show element details
 [1] "Version"           "job"               "StartTime"         "Iterations"       
 [5] "LogLikelihood"     "Residual.DF"       "AkaikeIC"          "BayesianIC"       
 [9] "ICparameter.Count" "Traits"            "keyres"            "waldT"            
[13] "components"        "FinishAt"          "Converge"          "esr"              
[17] "coef"              "yht"               "pred"              "evp"              
[21] "vpc0"              "vpc"               "vpv.mat"     
> # variance componets
> Var(res11)
  NO     Term   Sigma      SE  Z.ratio
1  1 Residual 1589.00 100.290 15.84405
2  2      Mum  132.59  56.459  2.34843
> # wald table for fixed effects
> wald(res11)

Wald tests for fixed effects.
			Wald F statistics
Source of Variation           NumDF     DenDF     F-inc              P-inc 
 mu                               1            11253.79
 Rep                              4               36.34
> # summary of key results
> summary(res11)
 Akaike Information Criterion   4697.68 (assuming 2 parameters).
 Bayesian Information Criterion 4706.31

          Analysis of h3 

                         Wald F statistics
Source of Variation           NumDF     DenDF     F-inc              P-inc 
 mu                               1            11253.79
 Rep                              4               36.34

 Model_Term                     Order     Gamma         Sigma     Z_ratio  %C
 Mum                               55  0.834411E-01   132.587        2.35   0 P    
 Residual_units                   559   1.00000       1588.99       15.84
  Mum                                     55 effects fitted.

> ## run pin functions
> pin(res11,org=T) # show variance components sign
  NO     Term   Sigma      SE vcS
1  1 Residual 1589.00 100.290  V1
2  2      Mum  132.59  56.459  V2
> pin11(res11,h2~V2*4/(V2+V1),signif=T,all=T)
Loading required package: msm
variance components are as following:
      Term   Sigma      SE vcS
1 Residual 1589.00 100.290  V1
2      Mum  132.59  56.459  V2

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

pin results are as following:

      Term Estimate      SE Siglevel
1 Residual 1589.000 100.290       **
2      Mum  132.590  56.459       **
3       h2    0.308   0.125       **
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

> # for more than 2 parameters
> pin(res11,mulp=c(h2~V2*4/(V2+V1),
+                      Vp~V2+V1),signif=T)
variance components are as following:
      Term   Sigma      SE vcS
1 Residual 1589.00 100.290  V1
2      Mum  132.59  56.459  V2

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


      Term  Estimate       SE Siglevel
1 Residual 1589.0000 100.2900      ***
2      Mum  132.5900  56.4590       **
3       h2    0.3081   0.1254       **
4       Vp 1721.5770 106.4608      ***

> ## get solutions for fixed and random effects
> fix.eff<-coef(res11)$fixed
> head(fix.eff)
  Term Level  Effect    SE
1  Rep     1   0.000 0.000
2  Rep     2 -11.732 5.641
3  Rep     3  46.779 5.534
4  Rep     4  23.807 5.521
5  Rep     5   4.895 5.510
6   mu     1 231.494 4.389
> ran.eff<-coef(res11)$random
> head(ran.eff)
  Term Level Effect    SE
1  Mum 70048  3.792 7.781
2  Mum 70017 -0.903 8.574
3  Mum 70002  2.552 8.399
4  Mum 70010  3.489 7.786
5  Mum 70041  2.947 8.760
6  Mum 70031  6.637 8.393

> ## get predictions
> pred<-predict(res11)
> pred$heads
[1] "   Predicted values of  h3"
> head(pred$pred$pred1)
  Prediction Stnd_Error Ecode   Mum
1    248.037      7.742     E 70048
2    243.342      8.601     E 70017
3    246.796      8.412     E 70002
4    247.733      7.744     E 70010
5    247.191      8.803     E 70041
6    250.881      8.406     E 70031
> pred$ased
$ased1
[1] 12.0795

> # get other results
> IC(res11)  # for AIC and BIC
   DF     LogL     AIC     BIC
1 554 -2346.84 4697.68 4706.31
> converge(res11)
TRUE