Fit a Linear Mixed Effects Model
DESCRIPTION:
A method for the generic function lme() for objects inheriting from class formula.
USAGE: 
lme.formula(fixed, random, cluster, data=sys.parent(),
            start, est.method = c("RML","ML"), re.block, re.structure =
            c("unstructured", "diagonal", "identity", "compsymm", "ar1"),
            re.paramtr = c("matrixlog","logcholesky", "cholesky", "spherical",
            "givens"), serial.structure = c("identity", "ar1",
            "ar1.continuous", "compsymm", "ar2", "ma1", "ma2",
            "arma11"), serial.covariate = NULL,
            serial.covariate.transformation = c("rank.within.cluster", "none",
            "round", "global.rank"), var.function = c("identity",
            "power", "expon", "cte.power", "cte.expon"), var.covariate =
             NULL, var.estimate = T, na.action, control)
REQUIRED ARGUMENTS:
fixed:
a formula object, specifying the fixed effects part of the model, with the response on the left of a ~ operator, and the terms, separated by "+" operators, on the right. If data is given, all names used in the formula should be defined as variables in the data frame. A model with the intercept as the only fixed effect can be specified as ~1.
cluster:
an expression or formula object, specifying the experimental units over which the random effects vary. If cluster is given as a formula, it should have no left side to the ~expression.
OPTIONAL ARGUMENTS:
random:
a formula object, specifying the random effects part of the model, with the terms, separated by "+" operators, on the right of a ~ operator. This argument has the same general characteristics as fixed, but there will be no left side to the ~expression. A model with the intercept as the only random effect can be specified as ~1.

NOTE: random effects are always assumed to have mean zero. A nonzero mean can be specified by including an identical term in the fixed effects part of the model.

data:
a data frame in which to interpret the variables named in fixed, random, cluster, serial.covariate, and var.covariate. Default is the calling frame.
start:
an list containing initial values for the scaled variance-covariance matrix of the random effects (D) or the random effects variance-covariance parameter vector (theta) that defines the factorization of D, the serial structure parameters (alpha), and the variance function parameters (delta). Theta has to be consistent with the chosen structure and parametrization (see re.structure and re.paramtr below). See serial.structure and var.function below for definitions of alpha and delta respectively.
est.method:
estimation method; if equal to "ML", Maximum Likelihood is used, otherwise if "RML" is specified, Restricted Maximum Likelihood is used. Partial matching of arguments is used, so only the first character needs to be provided. Default is "RML".
re.block:
a list indicating how the random effects should be blocked. Random effects pertaining to different blocks are assumed to be independent. The random effects can be referenced either by their names, or the order in which they appear in the random formula

(NOTE: unless a -1 is used in random, an (Intercept) term will be included as the first random effect). Elements in the re.block list are vectors containing the names or numbers of the random effects. Within a vector all elements have to be of the same type (i.e. all names or all numbers). By default all random effects are included in the same block.

re.structure:
a character vector describing the variance-covariance structure in each block of the random effects. Predefined structure names are "unstructured" for a general variance-covariance matrix with q(q+1)/2 parameters (q = the number of random effects for the block); "diagonal" for independent random effects with possibly different variances (q parameters); "identity" for independent random effects with the same variance (1 parameter); "compsymm" for a compound symmetry structure, random effects with common variance and common covariance (2 parameters); "ar1" for a common variance and AR(1) autocorrelation structure (2 parameters). Users can define their own variance-covariance structure as a class with methods for the generic functions lme.re.factor() and lme.re.param(). The generic function lme.re.factor() is used to obtain the factorization of D from the parameter vector theta, and lme.re.param() is used to obtain theta from D (see lme.re.factor.ar1() and lme.re.param.ar1() for examples of the use of the generic functions). If re.structure is of length 1 the name will be used for all blocks. Partial matching against the predefined names is used, so only the first character needs to be provided. Default is "unstructured".
re.paramtr:
the parameterization to be used for the unstructured scaled variance-covariance matrices specified in re.structure; possible values are "cholesky" for the Cholesky decomposition, "logcholesky" for Cholesky using logs of the diagonal elements, "spherical" for spherical coordinates of columns of the Cholesky decomposition, "matrixlog" for the matrix logarithm, and "givens" for a Givens rotation approach to defining eigenvalues and eigenvectors. These variance-covariance parametrizations are described in Pinheiro and Bates (1996) - see references section. Partial matching of arguments is used, so only the first character needs to be provided. Default is "matrixlog".
serial.structure:
the variance-covariance structure to be used for the within-cluster errors. Predefined structure names are "identity" for independent errors with constant variance; "ar1" for an AR(1) time series structure with respect to an integer valued time covariate (see serial.covariate below); "ar1.continuous" for a continuous time version of an AR(1) structure with the correlation between two observations d units of time apart being given by (a ^ d) where a is a parameter assumed to be one-unit lag correlation; "compsymm" for a compound symmetry structure, within-cluster errors with common variance and common covariance; "ar2" for an AR(2) structure with respect to an integer valued time covariate; "ma1" for an MA(1) structure with respect to an integer valued time covariate; "ma2" for an MA(2) structure with respect to an integer valued time covariate; and "arma11" for an ARMA(1,1) structure with respect to an integer valued time covariate. Users can define their own serial structure as a class with methods for the generic functions lme.serial.factor(), lme.serial.transform(), and lme.serial.untransform(). The generic function lme.serial.factor() is used to obtain the correlation matrix of the within-cluster errors for a given cluster, lme.serial.transform() is used to transform the original parameter vector into an unrestricted form (used in the optimization algorithms), and lme.serial.untransform() does the opposite operation (see lme.serial.factor.ar1(), lme.serial.transform.ar1(), and lme.serial.untransform.ar1() for examples of such functions). Default is "identity".
serial.covariate:
a formula object specifying the time covariate to be used with serial.structure. This argument has the same characteristics as random. If serial.structure is not equal to "identity" and no serial.covariate is provided, the time covariate will be set to the order of the observations within each cluster.
serial.covariate.transformation:
a transformation to be applied to the time covariate defined in serial.covariate. Possible values are "rank.within.cluster" for the ranks of the time covariate within each cluster; "none" for no transformation; "round" for rounding the time covariate to the nearest integer; and "global.rank" for the overall ranks of the time covariate. Default is "rank.within.cluster".
var.function:
the variance function structure to be used for the within-cluster errors. Predefined structure names are "identity" for constant within-cluster error variance, "power" for a variance function that is proportional to some power of the absolute value of the predictions (var ~ |yhat|^(2a)), "expon" for a variance function that varies exponentially with the absolute value of the predictions (var ~ exp(a*|yhat|)); and "cte.power" for a variance function that varies with some power of the absolute value of the predictions, but approaches a constant value when the predictions go to zero (var ~ a + |yhat|^(2b)). Users can define their own variance functions as a class with methods for the generic functions lme.varfunc.transform(), lme.varfunc.untransform(), and lme.varfunc.weights(). The generic functions lme.varfunc.transform() and lme.varfunc.untransform() are used to transform the original parameters into an unrestricted form (used in the optimization algorithms) and vice-versa. lme.varfunc.weights() is used to obtain the variance function weights, defined as the inverse of the scaled error standard deviations. This function has the general form proposed by Davidian and Giltinan (1995) - see references below: 1/sqrt(g(yhat, parameters, covariates)), where yhat denotes the predictions, parameters is the parameter vector, and covariates are optional covariates defined in var.covariate (see below). See lme.varfunc.transform.power(), lme.varfunc.untransform.power(), and lme.varfunc.weights.power() for examples of such functions). Default is "identity".
var.covariate:
a formula object specifying the covariates to be used in the calculation of the variance function weights. These covariates will be passed to the generic function lme.varfunc.weights() as the model.frame corresponding to var.covariate, using the data frame data. This argument has the same characteristics as random.
var.estimate:
a logical variable indicating whether the parameters in the variance function structure should be estimated. If equal to T the variance function parameters will be estimated together with the other parameters in the model; if F the variance function parameters will be kept fixed at their initial values (which may be given in the start list). Default if T.
subset:
expression saying which subset of the rows of the data should be used in the fit. This can be a logical vector, or a numeric vector indicating which observation numbers are to be included, or a character vector of the row names to be included. All observations are included by default.
na.action:
a function to filter missing data. This is applied to the model.frame after any subset argument has been used. The default (with na.fail) is to create an error if any missing values are found. A possible alternative is na.omit, which deletes observations that contain one or more missing values.
control:
a list of control values for the estimation algorithm (see lme.control()).
VALUE:
Generic functions such as print(), plot() and summary() have methods to show the results of the fit. See lme.object for the components of the fit. The functions residuals(), coefficients() and fitted.values() can be used to extract some of its components.
DETAILS:
The computational methods are described in Lindstrom, M.J. and Bates, D.M. (1988). The general method is described in Laird, N.M. and Ware, J.H. (1982). The variance-covariance parametrizations are described in Pinheiro, J.C. and Bates., D.M. (1996). The different time series structures available for the serial.structure argument are described in Box, G.E.P., Jenkins, G.M., and Reinsel G.C. (1994). The use of variance functions for linear and nonlinear mixed effects models is presented in detail in Davidian, M. and Giltinan, D.M. (1995).
REFERENCES:
Lindstrom, M.J. and Bates, D.M. (1988). Newton-Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data. Journal of the American Statistical Association 83, 1014-1022.

Laird, N.M. and Ware, J.H. (1982). Random-Effects Models for Longitudinal Data. Biometrics 38, 963-974.

Pinheiro, J.C. and Bates., D.M. (1996). Unconstrained Parametrizations for Variance-Covariance Matrices Statistics and Computing to appear.

Box, G.E.P., Jenkins, G.M., and Reinsel G.C. (1994). Time Series Analysis: Forecasting and Control, 3rd Edition, Holden-Day.

Davidian, M. and Giltinan, D.M. (1995). Nonlinear Mixed Effects Models for Repeated Measurement Data. Chapman and Hall.

SEE ALSO:
lme.control , lme.lmList , lme.object , lmList .
EXAMPLES: 
# Example from Lindstrom and Bates (1988) J.A.S.A.
Ovary.fit <- lme(fixed = follicles ~ sin(2*pi*Time) + cos(2*pi*Time),
                 random = ~sin(2*pi*Time) + cos(2*pi*Time),
                 cluster = ~ Mare,
                 data = Ovary, re.block = list(1, 2:3),
                 re.structure = c("i","d"), serial.structure = "ar1",
                 var.function = "power")
Ovary.fit


# Returns the following:
Call:
  Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)
 Random:  ~ sin(2 * pi * Time) + cos(2 * pi * Time)
Cluster:  ~ mare
   Data: ovary.data


Variance/Covariance Components Estimate(s):


 Block: 1
  Structure: identity
  Standard Deviation(s) of Random Effect(s)
 (Intercept)
     2.96746


 Block: 2
  Structure: diagonal
  Standard Deviation(s) of Random Effect(s)
 sin(2 * pi * Time) cos(2 * pi * Time)
            1.27255         0.01309297


 Cluster Residual Variance: 8.63962


 Serial Correlation Structure: ar1
 Serial Correlation Parameter(s): 0.5663131


 Variance Function: power
 Variance Function Parameter(s): 0.0695264


Fixed Effects Estimate(s):
  (Intercept) sin(2 * pi * Time) cos(2 * pi * Time)
    12.16773          -2.968204         -0.8580597


Number of Observations: 308
Number of Clusters: 11