AFEchidna包的简单示例2--单变量的批量分析

新版AFEchidna里，已将echidna.batch()纳入echidna()函数里，通过参数batch、batch.G、batch.R来实现。

其用法如下：

// old syntax
echidna.batch(trait,fixed,random,residual,
               es0.path,mulT,mulN,predict)

简单示例如下：

library(AFEchidna)

# trait=c('h3','h4','h5')
# new syntax
res21<-echidna(trait=~h3+h4+h5,
  fixed=~1+Rep,random=~Mum,
  es0.file='fm.es0',batch=TRUE)

# old syntax
es0.path=r'(E:\课件\新建文件夹)'
res21<-echidna.batch(trait=trait,fixed='Rep',random='Mum',es0.path=es0.path)

运行过程：

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

Program starts running batch analysis ------ 

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

 Mon Feb 22 14:43:26 2021 
  Iteration     LogL eSigma NEDF
1         1 -2346.93   1576  554
2         2 -2346.84   1590  554
3         3 -2346.84   1589  554
Mon Feb 22 14:43:26 2021 LogL Converged 


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

 Mon Feb 22 14:43:26 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
Mon Feb 22 14:43:27 2021 LogL Converged 


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

 Mon Feb 22 14:43:27 2021 
  Iteration     LogL eSigma NEDF
1         1 -2667.81   5289  551
2         2 -2667.71   5340  551
3         3 -2667.71   5333  551
Mon Feb 22 14:43:27 2021 LogL Converged 

运行结果：

> 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.466        1     3
h4 3599.1 241.39 227.54 115.510        1     4
h5 5333.1 441.96 337.20 187.280        1     3
> pin(res21,mulp=c(h2~V2*4/(V1+V2)))
results as following: 
pin formula: 
h2 ~ V2 * 4/(V1 + V2)

terms:  V1-h2
      V1    SE
h3 0.308 0.125
h4 0.251 0.116
h5 0.306 0.124
> res21n<-b2s(res21)
> names(res21n)
[1] "h3" "h4" "h5"
> h3.res<-res21n[["h3"]]
> Var(h3.res)
  NO     Term   Sigma      SE   Z.ratio
1  1 Residual 1589.00 100.290 15.844052
2  2      Mum  132.59  56.466  2.348139
> pin(h3.res,org=T)
  NO     Term   Sigma      SE vcS
1  1 Residual 1589.00 100.290  V1
2  2      Mum  132.59  56.466  V2
> pin(h3.res,mulp=c(h2~V2*4/(V1+V2)))
variance components are as following:
      Term   Sigma      SE vcS
1 Residual 1589.00 100.290  V1
2      Mum  132.59  56.466  V2

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

      Term  Estimate       SE
1 Residual 1589.0000 100.2900
2      Mum  132.5900  56.4660
3       h2    0.3081   0.1254

> fix.eff<-coef(h3.res)$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(h3.res)$random
> head(ran.eff)
  Term Level Effect    SE
1  Mum 70048  3.793 7.781
2  Mum 70017 -0.903 8.575
3  Mum 70002  2.552 8.399
4  Mum 70010  3.489 7.786
5  Mum 70041  2.947 8.761
6  Mum 70031  6.638 8.394
> IC(h3.res)
   DF     LogL     AIC     BIC
1 554 -2346.84 4697.68 4706.31
> converge(h3.res)
TRUE