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分析流程:
零模板(.es0)与目标数据集在同个目录,然后将该目录路径通过参数es0.path传递给函数echidna();
指定分析性状(trait)、固定效应(fixed)、随机效应(random)和残差效应(residual);
指定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
转载本文请联系原作者获取授权,同时请注明本文来自林元震科学网博客。 链接地址: https://blog.sciencenet.cn/blog-1114360-1274677.html
上一篇:
AFEchidna包增加函数echidnaT 下一篇:
AFEchidna包的简单示例2--单变量的批量分析