Maximum likelihood estimation (MLE) and restricted maximum likelihood estimation (REML) are both methods used to estimate the parameters of a statistical model. The main difference between MLE and REML lies in the way they estimate the variance-covariance matrix of the random effects in mixed-effects models.
In MLE, the variance-covariance matrix of the random effects is estimated using the full likelihood, which includes both the fixed and random effects. However, the maximum likelihood estimate of the variance-covariance matrix can be biased downward due to the inclusion of the fixed effects in the likelihood function. This bias is referred to as the “degrees of freedom” problem, and it can lead to an overestimation of the significance of the fixed effects.
In REML, the variance-covariance matrix of the random effects is estimated using a likelihood function that is based only on the random effects. The fixed effects are eliminated from the likelihood function by taking conditional distributions of the data given the random effects. By eliminating the fixed effects from the likelihood function, REML avoids the degrees of freedom problem and produces an unbiased estimate of the variance-covariance matrix of the random effects.
In summary, MLE estimates both the fixed and random effects using the full likelihood function, while REML estimates only the variance-covariance matrix of the random effects using a likelihood function that eliminates the fixed effects. REML is preferred over MLE when the goal is to estimate the variance-covariance matrix of the random effects without bias. However, MLE is sometimes preferred when the goal is to estimate the fixed effects as well as the variance-covariance matrix of the random effects.
One famous example of using REML is in the analysis of variance (ANOVA) of a mixed-effects model. In this case, the goal is to estimate the variance components of the random effects, which can be used to make inferences about the variability among the different levels of the grouping factor.
For example, consider a study where the effect of a treatment is being tested on a response variable measured in different subjects over time. The subjects are considered as random effects, and the time and treatment effects are considered fixed effects. A mixed-effects model can be used to model this data, where the response variable is modeled as a function of time, treatment, and subject:
Y = Xβ + Zb + ε
where Y is the response variable, X and Z are design matrices for the fixed and random effects, β and b are the corresponding parameter vectors, and ε is the residual error.
In this case, the REML method can be used to estimate the variance components of the random effects, which represent the variability among the different subjects. The fixed effects, including the treatment and time effects, can be estimated using either REML or MLE.
REML is preferred over MLE in this case because the goal is to estimate the variability among the subjects without bias. By using REML, the degrees of freedom problem associated with MLE is avoided, which can lead to biased estimates of the variance components. The estimates obtained using REML can then be used to test the significance of the treatment effect and to make inferences about the variability among the different subjects.
The likelihood function for a mixed-effects model can be written as:
L(θ|y) = (2π)-n/2|V(θ)|^-1/2 exp[-1/2(y-μ(θ))’V(θ)^-1(y-μ(θ))]
where θ is the vector of unknown parameters, y is the vector of observed data, V(θ) is the variance-covariance matrix of the random effects, and μ(θ) is the vector of expected values for the response variable.
The REML estimator is obtained by maximizing the likelihood function with respect to the parameters after integrating out the random effects. This involves maximizing a modified likelihood function:
L’(θ|y) = (2π)-n/2|W(θ)|^-1/2 exp[-1/2(y-Xβ(θ))’W(θ)^-1(y-Xβ(θ)))]
where W(θ) = V(θ) - V(θ)X(X’V(θ)-1X)-1X’V(θ) is a matrix that depends on both the variance-covariance matrix of the random effects and the design matrix X.
To show that the REML estimator is unbiased, we need to show that the expected value of the estimator equals the true value of the parameter:
E(θREML) = θ
We can do this by taking the derivative of the log-likelihood function with respect to the parameters and evaluating it at the true value of the parameter:
d/dθ log L’(θ|y)|θ=θtrue = 0
This gives us the score equations, which can be used to solve for the maximum likelihood estimates of the parameters.
We can then take the expected value of the score equations and evaluate them at the true value of the parameter:
E[d/dθ log L’(θ|y)|θ=θtrue] = 0
This shows that the expected value of the score equations equals zero, which implies that the REML estimator is unbiased.
Sure, here is an example of how REML can be used to estimate the coefficients of a mixed-effects model:
Suppose we have a dataset of 100 patients who are measured for their blood pressure (BP) at four different time points (t=0, 1, 2, 3) before and after receiving a treatment. We want to estimate the effect of the treatment on BP while taking into account the correlation among the repeated measurements within each patient.
We can model this data using a mixed-effects model with patient-specific intercepts and slopes over time:
BPij = β0i + β1itij + β2Txij + εij
where BPij is the blood pressure measurement for patient i at time j, β0i and β1i are the patient-specific intercept and slope for time, β2 is the treatment effect, Txij is the treatment indicator variable (1 if patient i received the treatment at time j, 0 otherwise), and εij is the residual error.
To estimate the coefficients of this model using REML, we need to first estimate the variance-covariance matrix of the random effects (patient-specific intercepts and slopes). This can be done using the lme() function in R with the method argument set to “REML”:
#library(nlme)
#
## Fit mixed-effects model using REML
#model <- lme(BP ~ t*Tx, random = ~t|patient, data = data, method = "REML")
#
## Extract coefficient estimates
#summary(model)$tTable
#This code fits the mixed-effects model to the data using the lme() function from the nlme package in R. The random argument specifies that the intercept and slope for time are patient-specific random effects. The method argument is set to “REML” to use the REML method to estimate the variance-covariance matrix of the random effects.
The output of this code gives us the estimated coefficients for the fixed effects, including the treatment effect (β2), as well as the estimated variance-covariance matrix of the random effects. The estimated coefficient for the treatment effect tells us the expected change in blood pressure associated with receiving the treatment, while controlling for the effect of time and the correlation among the repeated measurements within each patient.
Note that the variance-covariance matrix of the random effects estimated using REML can be used to test the significance of the patient-level variability, as well as to make inferences about the correlation among the repeated measurements within each patient.