AFEchidna: a new package solving mixed linear model for plant and animal datasets

[This article was first published on R-posts.com, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Mixed linear models(MLM) are linear models with a combination of fixed and random effects to explain the degree of variation in interest traits, such as milk yield in cows or volume growth in forest trees.  MLM are widely used in the analysis in the progeny test data of plants and animals. . Nowadays, most software uses the Restricted Maximum Likelihood (REML) method to estimate the variance components of random effects, and then estimate the fixed effects and predict the random effects. Such genetic analysis software includes ASReml, Echidna, SAS, BLUPF90, and R packages sommer, breedR, etc. Echidna is a free software developed in 2018 by Professor Gilmour, the main developer of ASReml. It also uses REML method to estimate parameter values, and its syntax and function is very close to that of ASReml. It is the most powerful free software for animal and plant genetic assessment, but its usage is a little complicated, which may be difficult for ordinary users. 

Here, I released  AFEchidna, an R package based on Echidna software, and demonstrated how to use a mixed linear model to generate solutions for variance components, genetic parameters, and random effects BLUPs.

The main functions in AFEchidna:
  • get.es0.file: generate es0 file
  • echidna(): specify mixed linear model
  • Var(): output variance components
  • plot(): model diagnose plots
  • pin(): calculate genetic parameters
  • predict(): model predictions
  • coef(): model equation solutions
  • model.comp():  compare different models
  • update(): run new mode
library(AFEchidna) 
setwd("D:\\Rdata") 
get.es0.file(dat.file=" Provenance.csv")  # generate .es file
get.es0.file(es.file=" Provenance.es")    # generate .es0 file

Specified a mixed model:

m1.esr <- echidna( fixed=height~1+Prov,
                   random=~Block*Female,
                   residual=~units,
                   es0.file='Provenance.es0')

Output related results:
> Var(m1.esr)
          Term   Sigma       SE   Z.ratio
1     Residual 2.52700 0.131470 19.221115
2        Block 0.10749 0.089924  1.195343
3       Female 0.18950 0.083980  2.256490
4 Block:Female 0.19762 0.086236  2.291618
> pin(m1.esr, mulp=c(Va~4*V3,
+ Vp~V1+V3+V4,
+ h2~4*V3/(V1+V3+V4)), digit=5)
 Term Estimate      SE
1  Va  0.75801 0.33592
11 Vp  2.91415 0.14973
12 h2  0.26011 0.10990
> coef(m1.esr)$fixed
  Term Level    Effect         SE
1 Prov    11 0.0000000  0.0000000
2 Prov    12 -1.6656325 0.3741344
3 Prov    13 -0.6237406 0.3701346
4 Prov     0 -1.2201920 0.3769892
5   mu     1 11.5120637 0.3619919
> coef(m1.esr)$random %>% head
    Term Level    Effect        SE
1  Block     2 0.3110440 0.1854718
2  Block     3 0.1268193 0.1858290
3  Block     4 0.2055624 0.1858158
4  Block     5 -0.2918516 0.1866871
5  Block     0 -0.3515741 0.1878287
6 Female   191 -0.1633745 0.3433451

A easy way to run batch analysis:
mt.esr <- update(m1.esr,trait=~height+diameter+volume,batch=TRUE)
> Var(mt.esr)

V1-Residual; V2-Block; V3-Female; V4-Block.Female
Converge: 1 means True; 0 means FALSE.

              V1       V2      V3      V4    V1.se V2.se V3.se V4.se Converge maxit
height    2.5270 0.107490 0.18950 0.19762 0.13147 0.089924 0.08398 0.08624 1 4
diameter 16.8810 0.167130 0.80698 0.69100 0.87659 0.198100 0.41470 0.50067 1 5
volume    0.0037 0.000084 0.00022 0.00016 0.00019 0.000077 0.00010 0.00011 1 6
AFEchidna: a new package solving mixed linear model for plant and animal datasets was first posted on May 16, 2022 at 7:27 am.
To leave a comment for the author, please follow the link and comment on their blog: R-posts.com.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)