SIRE REMLM
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
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---------------------------------------------------------------------
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
   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   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

