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

Interested in publishing a one-time post on R-bloggers.com? Press here to learn how.

---

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