/********************** Introduction ************************/

A general reference is Searle,S.R., G. Casella and C.E.
  McCulloch. 1992. Variance Components.  Wiley.

The linear mixed model is:
                    Y = Xa + Zb + e
where b is composed of r subvectors corresponding to r
random factors.  X and Z are known design matrices,
a and b are called the random and fixed effect vectors respectively.

b is stacked as b1 | b2 | .... | br and Z is partitioned conformably
as Z1~Z2~...~Zr.  bi is assumed to be composed of iid N(0,si) random
variables.  The si's are the variance components.  Finally, it is assumed
that ej's are iid N(0,s0).  s0 is called the residual variance.

The EM procs will estimate a, b, s0, and s =  s1 | ... | sr.

/********************** INSTRUCTIONS ************************/

Program structure:

1. A user-defined data proc initializes or reads a data file and sets up
      scalars and matrices.
      The data proc returns:
      r = number of variance components
      Z = the design matrix for the random effects
      q = an r-vector containing the number of random effects associated with
          the r random factors
          Important:
          q lists, in order, the number of columns in each
          of the submatrices Zi where Z = Z1~Z2~...~Zr.
      X = design matrix for the fixed effects.
      y = response vector
      name = character vector, names of the fixed effect vectors
      dname = string used as a header to identify the data

      Three small examples of data procs for published data sets are
      provided in the program file.

2. Parameter estimates are computed via an EM algorithm.
      EM estimates are computed by calling ML to compute the ML estimates or
      by calling REML which computes the restricted ML estimates of s and s0
      and ML estimates of a are also computed.

      Both procs are passed the following variables
      s = an initial estimate of the variance components s1 ... sr (a guess
          should be sufficient)
      r,Z,q,X,y (explained above),
      test = criterion (e.g., 1E-5) used to determine convergence
      lflag = if 1, compute the likelihood at each iteration.  It is not
          necessary to do so, and is time-consuming for large problems.
      Both procs return the parmeter estimates and a prediction of b = E[b | y].

3. Standard errors are computed by one of two algorithms:
      IOML - use if ML was used to compute parameter estimates
      IOREML - use if REML was used to compute parameter estimates

      Both procs return
      ion = observed information matrix
      se = vector of se's associated with a|s0|s;
         The proc prints out the square root of the estimates and the
         large sample (Delta method) estimates of the standard errors.

  Notes: the program is not designed for large problems.  More than
           150 random effects will take a considerable length of time
           until convergence.  If variance component is nearly zero,
           many iterations may be necessary until convergence.

Please contact me if you have any problems.