AFEchidna包里的函数echidna()可以达到在R中运行Echidna软件的目的。
其用法(新版本)如下:
echidna(es0.file,trait,fixed,random,residual,
foldN,delf,mulT,met,cycle,
predict,vpredict,
qualifier,jobqualf)
新旧版本的区别在于:使用了参数es0.file,替代了原来的es0.path;增加了参数met,用于多点试验数据。
其用法(旧版本)如下:
echidna(trait,fixed,random,residual,
es0.path,foldN,delf,
mulT,cycle,predict,qualifer)
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)在新版本的AFEchidna里可以直接生成,但仍需要对数据域进行手动设置。
echidna分析流程:
- 零模板(.es0)与目标数据集在同个目录,然后将该目录路径通过参数es0.path传递给函数echidna();
- 指定分析性状(trait)、固定效应(fixed)、随机效应(random)和残差效应(residual);
- 指定predict等参数;
完成上述步骤后,即可运行echidna()进行分析。
简单示例如下:
library(AFEchidna)
//new syntax
res11<-echidna(fixed=h3~1+Rep,random=~Mum,
residual=NULL,
predict=c('Mum'),
es0.file='fm.es0')
// old syntax
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
新版本的echidna()里将增加data.file、message和mulp等参数,用于增加特定功能。