REMLM

Software by Paul VanRaden, 1986

Download remlm.zip (959 KB)

remlm.f90 performs Restricted Maximum Likelihood (REML) estimation of variance components using matrices created by absorbm.f90.

ainv.f90 can be used to calculate the inverse of the numerator relationship matrix (A⁻¹) and sire list to be used by REML program.

The zipped file contains programs, example files, and executable (zipped folder of makefile, ainv.f90, absorbm.f90, remlm.f90), input folder (phenoreml.dat, pedigree.dat, traits.txt) and output folder (ainv.coef, sire.list, pass.dat, remlm.out and others)

References

VanRaden, P. M. 1986. Computational strategies for estimation of variance components. PhD thesis. Department of Animal Science, Iowa State University, Ames. https://doi.org/10.31274/rtd-180813-11154

Neupane, M, Hutchison, JL, Cole, JB, Van Tassell, CP, VanRaden, P.M. 2023. Genomic evaluation of late-term abortion in cows recorded through Dairy Herd Improvement test plans. JDS communications 4 (5), 354-357. [example of current use]

Instructions

Sire model multi-trait REML:

remlm.f90 performs Restricted Maximum Likelihood (REML) estimation

of variance components using matrices created by absorbm.f90.

ainv.f90 can be used to calculate the inverse of the numerator

relationship matrix (A⁻¹) and sire list to be used by REML program.

It estimates single and multitrait variance components,

correlations (genetic, phenotypic), sire PTAs and reliability.

These programs have been used in recent peer-reviewed studies

to estimate heritability for heifer livability and late-term abortion

in dairy cattle (Neupane et al., 2021 and 2023).

DOWNLOAD

Programs, example files, and executable (zipped folder of makefile,

ainv.f90, absorbm.f90, remlm.f90), input folder (phenoreml.dat,

pedigree.dat, traits.txt) and output folder (ainv.coef, sire.list,

pass.dat, remlm.out and others)

PROGRAMS

Make executable files for ainv.f90, absorbm.f90, and remlm.f90 using makefile

# Make file for Paul's 1986 REML programs

# To speed processing, use -O2 instead of -g option in compile

all: absorbm remlm ainv

absorbm: absorbm.f90

ifort absorbm.f90 -o absorbm -s -check bounds

remlm: remlm.f90

ifort remlm.f90 -o remlm -O1 -s

ainv: ainv.f90

ifort ainv.f90 -o ainv -s -check bounds

Step 1: ainv: The Fortran program ainv.f90 calculates the inverse of

the numerator relationship matrix (A⁻¹) using pedigree information.

Input file: pedigree.dat

Output files: ainv.coef (A⁻¹ matrix), sire.list

Step 2: absorbm: The Fortran program absorbm.f90 is designed to absorb

fixed effects (herd-year-season) from phenotypic data and produce output

equations needed by another program, REMLM for variance component estimation.

It preprocesses data for REML analysis by simplifying the design matrix and reducing

computational cost through absorption techniques.

Input files: phenoreml.dat, traits.txt, sire.list (from step 1 output)

Output file: pass.dat (binary)

Step 3: remlm: The Fortran program remlm.f90 performs REML estimation

of variance components and supports both single and multi-trait analyses. It utilizes

preprocessed input from absorbm and ainv programs.

Input files: pass.dat (binary output from absorbm), ainv.coef (output from ainv),

sire.list (same as above)

Output files: remlm.out, sire.rel, sire.pta, envrionment.cor,

genetic.cor, genphen.cor

Summary

Programs Input Output

ainv.f90 pedigree.dat ainv.coef, sire.list

absorbm.f90 phenoreml.dat, traits.txt, sire.list pass.dat (binary)

remlm.f90 pass.dat, ainv.coef, sire.list remlm.out (REML estimates)

# Example Data Set

(also included in zip files as individual files phenoreml.dat, pedigree.dat, traits.txt)

Phenoreml.dat

hys bull y(1) y(2)

1 1 21.452 99.252

1 1 21.193 102.940

1 1 20.102 97.053

1 1 15.610 101.115

1 4 21.753 99.498

1 4 17.635 98.151

1 4 17.301 98.208

1 5 18.115 103.320

1 5 19.691 99.040

1 5 18.096 96.624

1 7 21.695 99.011

1 7 19.747 101.699

1 7 18.190 99.839

1 7 18.640 100.475

1 8 21.863 100.338

1 8 20.278 102.966

1 8 19.036 95.209

1 9 22.753 107.425

2 6 20.535 97.043

2 6 22.340 108.602

2 6 21.858 100.212

2 9 20.744 96.114

2 9 25.811 98.612

2 9 20.141 100.857

3 2 19.942 100.198

3 2 17.885 98.922

3 3 20.602 101.337

3 4 22.583 98.136

3 4 20.959 98.389

3 5 13.799 95.213

3 7 20.633 101.923

3 7 19.461 96.552

3 8 20.785 100.059

3 8 20.173 96.874

3 8 20.248 100.059

3 8 19.835 100.559

4 8 18.479 99.616

4 9 21.499 102.735

4 9 20.597 105.478

4 10 23.001 98.764

4 10 21.717 98.226

5 3 17.992 97.635

5 3 18.390 105.155

5 6 17.793 102.643

5 6 19.624 96.256

5 6 21.908 100.518

5 6 20.833 98.929

5 8 21.593 99.286

5 9 20.695 101.417

5 9 22.094 96.795

5 9 21.574 101.871

5 9 22.054 100.140

6 2 22.435 98.995

6 2 21.791 97.937

6 2 17.319 96.914

6 5 17.613 105.889

6 5 17.380 99.893

6 6 17.250 106.584

6 6 19.635 101.783

6 8 15.466 99.050

6 8 20.455 100.039

6 8 19.435 104.977

6 10 22.648 95.979

6 10 26.222 96.139

7 3 26.296 100.342

7 6 20.806 97.503

7 6 20.763 98.159

7 10 22.614 98.861

7 10 23.448 102.157

8 3 18.982 98.954

8 3 20.153 96.993

8 3 21.037 103.463

8 3 19.591 98.970

8 4 21.092 93.855

8 4 22.272 97.457

8 4 21.337 97.699

8 4 20.634 101.801

8 5 20.079 100.248

8 5 16.552 100.185

8 5 18.584 105.355

8 5 15.844 98.968

8 8 23.231 100.079

8 8 17.003 102.937

8 9 19.926 100.963

9 5 19.556 94.021

9 5 18.436 101.534

9 5 20.937 98.286

9 5 22.219 94.895

9 6 18.650 101.436

9 7 22.398 103.256

9 7 21.380 99.087

9 7 21.741 102.051

9 8 19.918 102.196

9 8 18.795 97.540

9 8 16.233 101.117

9 8 18.818 102.566

9 9 20.229 99.627

9 10 21.042 100.282

10 2 22.067 93.227

10 2 21.902 92.366

10 5 17.589 101.728

10 7 20.218 98.640

10 8 21.940 96.117

10 8 19.244 99.651

10 8 18.784 100.362

10 9 19.100 102.945

10 9 22.482 105.683

--------------------------------------------------------------------

pedigree.dat

bull sire mgs mggs sex birth_year

1 5 6 0 M 1974

2 7 5 0 M 1980

3 5 7 0 M 1975

4 8 0 0 M 1975

5 0 0 0 M 1968

6 0 0 0 M 1962

7 0 0 0 M 1964

8 0 0 0 M 1962

9 0 0 0 M 1965

10 0 0 0 M 1965

LICENSE

This software is public domain and was developed with U.S. taxpayer

funding. Accurate results are not guaranteed. Please report any bugs

to jason.graham@uscdcb.com or mahesh.neupane@usda.gov. You may modify,

improve, use, and redistribute the code to anyone for any purpose.

Or, you can ask Jason Graham or Mahesh Neupane to make changes that could

benefit U.S. evaluations and other users.

REFERENCE

VanRaden, P. M. 1986. Computational strategies for estimation of variance components.

PhD thesis. Department of Animal Science, Iowa State University, Ames.

https://doi.org/10.31274/rtd-180813-11154

Neupane, M. et al. 2021. Genomic evaluation of dairy heifer livability. Journal of

Dairy Science, Volume 104, Issue 8, 8959-8965. https://doi.org/10.3168/jds.2020-19687

Neupane, M. et al. 2023. Genomic evaluation of late-term abortion in cows recorded

through Dairy Herd Improvement test plans. JDS Communications, Volume 4, Issue 5,

354-357. https://doi.org/10.3168/jdsc.2022-0341

Project team members - 2025

Paul M. VanRaden (paul.vanraden@usda.gov)

Jason Graham (Jason.graham@uscdcb.com)

Mahesh Neupane (Mahesh.neupane@usda.gov)

#### Original email with details about program and example files ####

SUBJECT: UPDATE OF REML PROGRAMS

December, 1994

Dear Researcher

This letter will be followed by two additional mailings

which will contain the source codes for the updated programs

ABSORBD and REMLD, which compute multiple-trait restricted

maximum likelihood estimates of genetic variances and covariances.

The programs are more portable than the 7-86 version because they

contain their own subroutines and no longer depend on IMSL. They

are in double precision, should run faster, and are easier to use.

An example data set is included for demonstration and to allow

verification that the programs are working properly on your

system. The example data set, list of sires, and inverse of the

relationship matrix among the sires follow the example results.

If you detect problems please let me know about them. Also, if

you have suggestions for improvements I do from time to time update

the programs. Instructions for using the programs are next and in the

FORTRAN code, near the top of the programs. Results expected from

running the programs as sent on the example data follow.

Sincerely,

Paul VanRaden

USDA, Beltsville, MD 20705

e-mail paul@aipl.arsusda.gov

------------------------------------------------------------------------

C PROGRAM ABSORBD, CREATES AND OUTPUTS MATRICES

C NEEDED BY REMLD, PAUL VANRADEN 6-88

C Fortran program to absorb herd-year-season equations and to

C output sire and group equations and other information needed

C by program REMLD. Program works for any number of traits

C provided they are all measured on every animal. Each trait

C is expected to behave according to a model of the form:

C y = mu + herd-year-season + genetic group + sire(within

C genetic group) + error.

C NOTES ON USE OF ABSORBD:

C 1) If genetic groups are used, the genetic group to which each

C sire belongs should accompany the sire ID's in unit 11.

C Group numbers must be consecutive integers starting with 1.

C Do not set NGM1 equal to 0 or program may not compile.

C 2) Group equations are assumed to be full rank after deleting the

C last row. This requires connectedness among all the groups.

C 3) Format statements 4 and 13 must be specified for reading the

C sire list and the data.

C 4) Other fixed effects known to influence the data may be estimated

C ahead of time and subtracted from the data as correction factors

C as if these were known. This works well for fixed effects

C having few degrees of freedom.

C 5) Be sure that the traits included are not linear functions of

C each other. This would cause a singular error variance-covar.

C matrix which program REMLD is not prepared to deal with.

C 6) Storage space for matrices output to unit 12 is proportional to

C number of sires squared: Mbytes = (number of sires / 300)**2

C---------------------------------------------------------------------

C NOTES ON USE OF REMLD:

C 1) First run program ABSORBD to create the input file.

C 2) Set values of the parameters. Number of sires, traits, and

C groups in REMLD must be the same as declared in ABSORBD.

C 3) The relationship matrix among sires is read twice at

C statements 17 and 163

C 4) To obtain exact REML estimates, set JTYPE = 1.

C 5) Memory required is proportional to the number of sires squared

C Mbytes = (number of sires / 250)**2

C 6) Time required is proportional to number of sires cubed

C Hours = (number of sires / 1500)**3

C----------------------------------------------------------------------

Example Output from ABSORBD:

TOTAL OF 0 OBSERVATIONS DROPPED BECAUSE SIRE WAS NOT FOUND IN LIST

RAW MEANS OF Y VARIABLES = 20.214000 99.860084

NO. OF RECORDS, FILLED SUBCLASSES, HYS, AND SIRES

107 45 10 10

N SIRE ID GROUP #PROG #HYS EFF. NUM.

1 1 1 4 1 3.11

2 2 1 7 3 5.47

3 3 1 8 4 6.29

4 4 1 9 3 7.10

5 5 1 15 6 11.76

6 6 2 12 5 7.84

7 7 2 10 4 8.02

8 8 2 21 8 15.72

9 9 2 14 7 9.61

10 10 2 7 4 5.00

--------------------------------------------------------------------

Example of Output from REMLD, with RELATD = 1

WITHIN SMALLEST SUBCLASS ERROR VARIANCES

1 2.998816

2 7.564156

CHOLESKY DECOMPOSITION OF A-INVERSE DONE

INVERSION OF CHOLESKY MATRIX DONE

LARGE MATRIX MULTIPLICATION 1 DONE

LARGE MATRIX MULTIPLICATION 2 DONE

EIGENVECTORS OF L ZSZ L' OBTAINED

RATIOS FOR DIAGONAL

9.000 9.000

ROUND 1 COMPLETED.

RATIOS FOR DIAGONAL

10.272 3.547

ROUND 2 COMPLETED.

...

RATIOS FOR DIAGONAL

25.119 3.000

ROUND 12 COMPLETED.

RATIOS FOR DIAGONAL

25.118 3.000

ROUND 13 COMPLETED.

FINAL ESTIMATES:

TRAIT ERROR VAR. SIRE VAR. HERITABILITY

1 3.673173 1.030101 0.876072

2 8.632654 0.924886 0.387081

STANDARD ERRORS OF CORRESPONDING PARAMETERS ABOVE

1 0.659721 0.881239 0.598082

2 1.550470 1.208737 0.461219

GENETIC CORRELATIONS ARE ABOVE DIAGONAL

PHENOTYPIC CORRELATIONS ARE BELOW DIAGONAL

TRAIT 1 2

1 1.000-0.748

2 -0.160 1.000

GROUP EFF.NUM. GROUP SOLUTIONS

1 1 20.05 -0.002878 -2.046445

2 2 20.05 0.000000 0.000000

STANDARD ERRORS OF GROUP SOLUTIONS

1 1 20.05 0.470746 1.089458

2 2 20.05 0.000000 0.000000

NOTICE: PTA'S CONTAIN GENETIC GROUP SOLUTIONS

SIRE ID GRP EFF.NUM. PREDICTED TRANSMITTING ABILITIES

1 1 2.95 -0.730242 -1.421113

2 1 4.64 0.289997 -2.493432

3 1 5.59 -0.581313 -1.311632

4 1 6.61 0.333423 -2.540243

5 1 9.09 -1.422976 -0.960622

6 2 7.59 -0.843896 0.631660

7 2 7.05 0.302742 -0.264419

8 2 11.82 -0.378670 -0.116608

9 2 9.15 0.245945 0.401392

10 2 4.88 1.477345 -1.305173

SIRE ID GRP EFF.NUM. RELIABILITIES

1 1 2.95 0.535164 0.349211

2 1 4.64 0.616606 0.411268

3 1 5.59 0.653892 0.449778

4 1 6.61 0.678446 0.472149

5 1 9.09 0.749174 0.551064

6 2 7.59 0.684489 0.455021

7 2 7.05 0.691266 0.476561

8 2 11.82 0.779292 0.584096

9 2 9.15 0.719531 0.494974

10 2 4.88 0.577808 0.343338

--------------------------------------------------------------------

Example of Output from REMLD, with RELATD = 0

RATIOS FOR DIAGONAL

23.393 3.087

ROUND 17 COMPLETED.

FINAL ESTIMATES:

TRAIT ERROR VAR. SIRE VAR. HERITABILITY

1 3.624819 0.998751 0.864052

2 8.644632 0.912359 0.381860

STANDARD ERRORS OF CORRESPONDING PARAMETERS ABOVE

1 0.651037 0.793555 0.551810

2 1.552621 1.090239 0.417385

GENETIC CORRELATIONS ARE ABOVE DIAGONAL

PHENOTYPIC CORRELATIONS ARE BELOW DIAGONAL

TRAIT 1 2

1 1.000-0.725

2 -0.158 1.000

GROUP EFF.NUM. GROUP SOLUTIONS

1 1 20.05 -0.523385 -1.622945

2 2 20.05 0.000000 0.000000

STANDARD ERRORS OF GROUP SOLUTIONS

1 1 20.05 0.470619 1.080333

2 2 20.05 0.000000 0.000000

NOTICE: PTA'S CONTAIN GENETIC GROUP SOLUTIONS

SIRE ID GRP EFF.NUM. PREDICTED TRANSMITTING ABILITIES

1 1 2.95 -0.585373 -1.456923

2 1 4.64 0.115101 -2.412290

3 1 5.59 -0.619525 -1.117504

4 1 6.61 0.191307 -2.260999

5 1 9.09 -1.740281 -0.849488

6 2 7.59 -0.950409 0.690714

7 2 7.05 0.092486 -0.076823

8 2 11.82 -0.575484 0.025303

9 2 9.15 0.150776 0.509325

10 2 4.88 1.363753 -1.213567

SIRE ID GRP EFF.NUM. RELIABILITIES

1 1 2.95 0.417085 0.228152

2 1 4.64 0.513431 0.311145

3 1 5.59 0.548770 0.347791

4 1 6.61 0.590912 0.387358

5 1 9.09 0.657214 0.461200

6 2 7.59 0.614750 0.416278

7 2 7.05 0.581846 0.390974

8 2 11.82 0.696574 0.517120

9 2 9.15 0.651230 0.458918

10 2 4.88 0.532258 0.324567

----------------------------------------------------------------------

Example Data Set:

hys sire y(1) y(2)

1 1 21.452 99.252

1 1 21.193 102.940

1 1 20.102 97.053

1 1 15.610 101.115

1 4 21.753 99.498

1 4 17.635 98.151

1 4 17.301 98.208

1 5 18.115 103.320

1 5 19.691 99.040

1 5 18.096 96.624

1 7 21.695 99.011

1 7 19.747 101.699

1 7 18.190 99.839

1 7 18.640 100.475

1 8 21.863 100.338

1 8 20.278 102.966

1 8 19.036 95.209

1 9 22.753 107.425

2 6 20.535 97.043

2 6 22.340 108.602

2 6 21.858 100.212

2 9 20.744 96.114

2 9 25.811 98.612

2 9 20.141 100.857

3 2 19.942 100.198

3 2 17.885 98.922

3 3 20.602 101.337

3 4 22.583 98.136

3 4 20.959 98.389

3 5 13.799 95.213

3 7 20.633 101.923

3 7 19.461 96.552

3 8 20.785 100.059

3 8 20.173 96.874

3 8 20.248 100.059

3 8 19.835 100.559

4 8 18.479 99.616

4 9 21.499 102.735

4 9 20.597 105.478

4 10 23.001 98.764

4 10 21.717 98.226

5 3 17.992 97.635

5 3 18.390 105.155

5 6 17.793 102.643

5 6 19.624 96.256

5 6 21.908 100.518

5 6 20.833 98.929

5 8 21.593 99.286

5 9 20.695 101.417

5 9 22.094 96.795

5 9 21.574 101.871

5 9 22.054 100.140

6 2 22.435 98.995

6 2 21.791 97.937

6 2 17.319 96.914

6 5 17.613 105.889

6 5 17.380 99.893

6 6 17.250 106.584

6 6 19.635 101.783

6 8 15.466 99.050

6 8 20.455 100.039

6 8 19.435 104.977

6 10 22.648 95.979

6 10 26.222 96.139

7 3 26.296 100.342

7 6 20.806 97.503

7 6 20.763 98.159

7 10 22.614 98.861

7 10 23.448 102.157

8 3 18.982 98.954

8 3 20.153 96.993

8 3 21.037 103.463

8 3 19.591 98.970

8 4 21.092 93.855

8 4 22.272 97.457

8 4 21.337 97.699

8 4 20.634 101.801

8 5 20.079 100.248

8 5 16.552 100.185

8 5 18.584 105.355

8 5 15.844 98.968

8 8 23.231 100.079

8 8 17.003 102.937

8 9 19.926 100.963

9 5 19.556 94.021

9 5 18.436 101.534

9 5 20.937 98.286

9 5 22.219 94.895

9 6 18.650 101.436

9 7 22.398 103.256

9 7 21.380 99.087

9 7 21.741 102.051

9 8 19.918 102.196

9 8 18.795 97.540

9 8 16.233 101.117

9 8 18.818 102.566

9 9 20.229 99.627

9 10 21.042 100.282

10 2 22.067 93.227

10 2 21.902 92.366

10 5 17.589 101.728

10 7 20.218 98.640

10 8 21.940 96.117

10 8 19.244 99.651

10 8 18.784 100.362

10 9 19.100 102.945

10 9 22.482 105.683

--------------------------------------------------------------------

List of bulls (sires of cows in data and tie ancestors if needed)

pedtyp biryr group bull sire mat.gsire

1110 74 1 1 5 6

1110 80 1 2 7 5

1110 75 1 3 5 7

1100 75 1 4 8 0

1000 68 1 5 0 0

1000 62 2 6 0 0

1000 64 2 7 0 0

1000 62 2 8 0 0

1000 65 2 9 0 0

1000 65 2 10 0 0

---------------------------------------------------------------------

Inverse of the relationship matrix among the sires,

row, column, coefficient

1 1 1.454545

1 5 -0.727273

1 6 -0.363636

2 2 1.454545

2 5 -0.363636

2 7 -0.727273

3 3 1.454545

3 5 -0.727273

3 7 -0.363636

4 4 1.333333

4 8 -0.666667

5 1 -0.727273

5 2 -0.363636

5 3 -0.727273

5 5 0.363636

5 5 0.090909

5 5 0.363636

5 5 1.000000

5 6 0.181818

5 7 0.181818

6 1 -0.363636

6 5 0.181818

6 6 0.090909

6 6 1.000000

7 2 -0.727273

7 3 -0.363636

7 5 0.181818

7 7 0.363636

7 7 0.090909

7 7 1.000000

8 4 -0.666667

8 8 0.333333

8 8 1.000000

9 9 1.000000

10 10 1.000000