用R进行Echidna程序结果的pin函数运算

Echidna的分析结果可以通过命令符!RML输出到_e.R文件，然后再读入R中：

> res<-esRT(trace=T)
Loading barley8Read 15 items

 Mon Dec 21 14:36:58 2020 
  Iteration     LogL   Sigma2 NEDF
1         1 -873.746 37956.08  149
2         2 -861.084 15412.44  149
3         3 -865.915 52600.29  149
4         4 -846.333 44066.98  149
5         5 -836.973 37188.23  149
6         6 -828.457 35727.69  149
7         7 -827.608 45413.55  149
8         8 -827.603 45240.25  149
9         9 -827.603 45226.23  149
Mon Dec 21 14:36:59 2020LogL Converged

Sussessfully input Echidna runs results into R.

原始pin结果，从Echidna的vpredict语句得到结果，如下：

> pin0(res,org=T)

 Warning: The analysis was performed on the GAMMA scale.
Check all parameters have been converted to the Sigma scale.
   1 Residual                                            45226.0        16819.0
   2 variety                                             15602.0        5039.8
   3 units                                               5137.5        1768.5    
   4 ar1(column)                                         38257.0       0.68194E-01
   5 ar1(row)                                            31346.0       0.10038
 W  components
 VVP matrix for parameters 1 to 5 written to componv
                    with values and labels in componc
 F Vp variety + units
   6 Vp 2                                                20739.0        5231.1
 F Vp1 variety + units + Residual
   7 Vp1 2                                               65966.0        17944.0    
 H h2p variety Vp
     h2p 2        = variety    2/Vp 2       6 =         0.75228       0.09079
 H h2p1 variety Vp1
     h2p1 2       = variety    2/Vp1 2      7 =         0.23652       0.08497
 Notice: The parameter estimates are followed by
          their approximate standard errors.

Var函数可输出方差分量，pin函数可以计算遗传参数：

> AfEchidna::Var(res)
            Term Level    Gamma       Sigma      std.err Z.ratio %ch const
1        variety    25 0.344976 1.56020e+04 5.032903e+03    3.10   0 P    
2          units   150 0.113595 5.13748e+03 1.771545e+03    2.90   0 P    
3    ar1(column)    15 0.845908 8.45908e-01 6.821839e-02   12.40   0 P    
4       ar1(row)    10 0.693088 6.93088e-01 1.004475e-01    6.90   0 P    
5 Residual_units   150 1.000000 4.52262e+04 1.681271e+04    2.69   0      
> pin(res,h2p~V1/(V1+V2))

   Estimate    SE
h2p   0.752 0.091
> pin(res,h2p1~V1/(V1+V2+V5))

    Estimate    SE
h2p1   0.237 0.085

pin函数计算结果与vpredict结果一致。但在R中的pin函数，可以计算更为复杂的公式，主要是涉及分子式的加减，这是vpredict语句无法处理的。例如，计算culius遗传力的遗传力标准误差，示例如下：

> pin(res,h2c~sqrt(1-57.5^2/V1))

    Estimate    SE
h2c    0.888 0.039

运行pin的批量计算：

> pin.batchS(res,mulp=c(h2~V1/(V1+V2),
+   h2c~sqrt(1-57.5^2/V1)),signif=T)
variance components are as following:
            Term       Sigma      std.err vcS
1        variety 1.56084e+04 5.034968e+03  V1
2          units 5.13988e+03 1.766282e+03  V2
3    ar1(column) 8.45929e-01 6.822008e-02  V3
4       ar1(row) 6.93063e-01 1.004439e-01  V4
5 Residual_units 4.52403e+04 1.681796e+04  V5

pin formula: 
h2 ~ V1/(V1 + V2)
h2c ~ sqrt(1 - 57.5^2/V1)

pin results are as following:

    Esmate     SE siglevel
h2  0.7523 0.0907      ***
h2c 0.8878 0.0385      ***
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1