MLE and REML Estimates
- 6 minsThe Restricted Maximum Likelihood Estimates (REML)
Estimate \(\sigma^2\),\(\tau^2\), and the fixed effects \(\beta\) via REML
- I use lmer function in lme4 package to fit the model to obtain the REML estimates of the variance components $\theta = (\sigma^2, \tau^2)$ and the REML estimates of the fixed effects $\beta$. \(\hat{\theta}_{\text{REML}} = (\hat{\sigma}^{2}_{\text{REML}},\hat{\tau}^{2}_{\text{REML}}) = (42106.61, 5850.373)\)
Intercept | year1 | AGE | SHARES | REV | INC |
---|---|---|---|---|---|
513.56951705 | 78.93746000 | -0.72576169 | 0.05438211 | 13.90654655 | 68.48501834 |
Obtain standard error for the REML for the fixed effects
- We know that
- Hence, first of all, we can obtain the REML estimator of \(V(\theta)\) by replacing the variance components $\theta$ by the REML estimators of variance components \(\hat{\theta}_{\text{REML}}\).
- We call
the REML estimator of V.
- Then we plug in \(\hat{V}_{\text{REML}}\) into
- So we obtain the variance of REML estimates of fixed effects.
- Last we take the diagonals of \(\hat{S}\) and the square root of that are the standard errors for the REML estimate of fixed effects. The \(\beta_{j}\)’s standard error is the jth element of the square root of the diagonal of S.
Intercept | year1 | AGE | SHARES | REV | INC |
---|---|---|---|---|---|
282.79806838 | 15.29754655 | 4.93312803 | 0.02941859 | 6.09109415 | 130.97186340 |
Obtain the standard errors for the variance-component estimates.
Note that the R output does not include standard errors for the REML estimates of the variance components.
There are two methods to obtain the s.e. for the variance-component estimates.
Asymptotic Covariance Matrix (ACM) method
- The REML estimators are consistent and aymptoticallly normal under certain conditions. And asymptotic covariance matrix is equal to the inverse of the restricted Fisher Information matrix. So under regularity conditions and assuming V is twice differentiable with respect to \(\theta\).
- Let REML estimates of variance components be
\(\hat{\theta}_{\text{REML}} = (\hat{\sigma}^{2}_{\text{REML}},\hat{\tau}^{2}_{\text{REML}})\).
And we define
\[\hat{P} = \hat{V}^{-1} - \hat{V}^{-1}X(X'\hat{V}^{-1}X)^{-1}X'\hat{V}^{-1}\]- The REML estimate of the inverse of the restricted Fisher Information matrix is
$$ I^{-1}(\hat{\theta}_{\text{REML}}) = \frac{1}{2}\Big( \begin{matrix} \text{tr}(\hat{P}ZZ'\hat{P}ZZ') & \text{tr}(\hat{P}^{2}ZZ') \\ \text{tr}(\hat{P}^2ZZ') & \text{tr}(\hat{P}^{2}) \end{matrix} \Big) ^{-1} $$
- Then we can obtain
$$ I^{-1}(\hat{\theta}_{\text{REML}}) = \Big( \begin{matrix} a^{\ast} & b^{\ast} \\ c^{\ast} & d^{\ast} \end{matrix} \Big) $$
- Hence
Parametric Bootstrap (PB) method
- Also let REML estimates of variance components be
from the original data.
- For b in 1:B (repeat the following two steps for B times)
- Simulate a sample from this multivariate normal distribution \(N(X\hat{\beta},\hat{V(\theta)})\)
- Fit a Linear Mixed Model to that sample to find REML estimates of \(\theta = (\sigma^{2}, \tau^{2})\).
- Let \((\hat{\theta}^{(b)}_{\sigma^{2}},\hat{\theta}^{(b)}_{\tau^{2}})\) be the REML estimates of \(\theta = (\sigma^{2},\tau^{2})\) respectively from the sample.
- The bootstrap estimates are
where
\(\bar{\theta}_{\sigma^{2}}\) is the mean of \(\hat{\theta}^{(b)}_{\sigma^{2}}\) : b=1,2,3,…,B
\(\bar{\theta}_{\tau^{2}}\) is the mean of \(\hat{\theta}^{(b)}_{\tau^{2}}\) : b=1,2,3,…,B
Compute the standard errors of REML estimates via ACM method.
- The standard errors for the REML estimates of the variance components by ACM method are
Compute the standard errors in REML estimates via PB method.
- The standard errors for the REML estimates of the variance components by PB method are
Compare the results of ACM and PB Comment on what we observe.
- The standard errors of \(\hat{\sigma}^{2}_{\text{\text{REML}}}\) are much more higher than that of \(\hat{\tau}^{2}_{\text{\text{REML}}}\) for both methods.
- The standard errors I obtained for \(\hat{\tau}^{2}_{\text{REML}}\) by two methods are more similar compared to the s.e. I obtained for \(\hat{\sigma}^{2}_{\text{REML}}\) by two methods.
- That is because of different degree of freedom of \(\hat{\tau}^{2}_{\text{REML}}$ and $\hat{\sigma}^{2}_{\text{REML}}\). The standard error of \(\hat{\tau}^{2}_{\text{REML}}\) is much more smaller since the degree of freedom is 99. The degree of freedom of \(\hat{\sigma}^{2}_{\text{REML}}\) is 49.
Repeat previous parts using MLE instead of REML.
variance Component Estimates
\[\hat{\theta}_{\text{ML}} = (\hat{\sigma}^{2}_{\text{ML}},\hat{\tau}^{2}_{\text{ML}}) = (37661.93, 5733.366)\]Intercept | year1 | AGE | SHARES | REV | INC |
---|---|---|---|---|---|
513.56951705 | 78.93746000 | -0.72576169 | 0.05438211 | 13.90654655 | 68.48501834 |
ML estimates of fixed effects \(\beta\)
\[\hat{\beta}_{\text{ML}} = (513.56951705,78.93746000,-0.72576169,0.05438211,13.90654655,68.48501834)'\]ACM method for standard errors of Variance Component
- Replace the \(\hat{P}\) with \(\hat{V}^{-1}\) in \(\hat{I}^{-1}(\hat{\theta})\). i.e.
$$ I^{-1}(\hat{\theta}_{\text{ML}}) = \frac{1}{2}\Big( \begin{matrix} \text{tr}(\hat{V}^{-1}ZZ'\hat{V}^{-1}ZZ') & \text{tr}(\hat{V}^{-2}ZZ') \\ \text{tr}(\hat{V}^{-2}ZZ') & \text{tr}(\hat{V}^{-2}) \end{matrix}\Big) ^{-1} $$
- Then we can obtain
$$ I^{-1}(\hat{\theta}_{\text{ML}}) = \Big( \begin{matrix} a^{\ast} & b^{\ast} \\ c^{\ast} & d^{\ast} \end{matrix} \Big) $$
- Hence
- The standard errors for the ML estimates of the variance components by ACM method are
PB method for standard errors of Variance Component
- replace the \(\sigma^{2}\)REML and \(\tau^{2}\)REML with \(\sigma^{2}\)ML and \(\tau^{2}\)ML
- The standard errors for the ML estimates of the variance components by PB method are
#R Code for STA232B-Project I
#==========================
#Organize Data
corps <- read.table("corps.txt", header = TRUE)
corps$$chairman <- rep(1:50,each = 2)
corps$$year <- corps$$year - 83
corps$$year <- as.factor(corps$$year)
n <- length(corps$$year)
#a
library(lme4)
library(HLMdiag)
fit_lmer = function(data, method = 'REML') {
mod = lmer(formula = y ~ year + AGE + SHARES + REV + INC + (1 | chairman),
data = data,
REML = (method == 'REML'))
list(mod = mod,
beta.hat = fixef(mod),
sigma2.hat = varcomp.mer(mod)['D11'],
tau2.hat = varcomp.mer(mod)['sigma2'])
}
REML = fit_lmer(corps)
# REML est. of \sigma^2
( REML$$sigma2.hat)
## D11
## 42106.61
# REML est. of \tau^2
( REML$$tau2.hat)
## sigma2
## 5850.373
# REML est. of \beta (fixed effects)
(beta.hat.reml = REML$$beta.hat)
## (Intercept) year1 AGE SHARES REV
## 513.56951705 78.93746000 -0.72576169 0.05438211 13.90654655
## INC
## 68.48501834
#b
summary(REML$$mod)$$coefficients[,2]
## (Intercept) year1 AGE SHARES REV
## 282.79806317 15.29754681 4.93312794 0.02941859 6.09109403
## INC
## 130.97186099
Z = as.matrix(getME(REML$$mod, "Z"))
X = as.matrix(getME(REML$$mod, "X"))
ZZ = Z %% t(Z)
V.hat = REML$$sigma2.hatZZ+REML$$tau2.hat diag(1,n)
Varbeta.hat = solve( t(X) %% solve(V.hat) %% X )
(s.e.beta.hat = sqrt(diag(Varbeta.hat)))
## (Intercept) year1 AGE SHARES REV
## 282.79806317 15.29754681 4.93312794 0.02941859 6.09109403
## INC
## 130.97186099
#d ACM method for REML
P.hat = solve(V.hat) - solve(V.hat) %% X %% solve(t(X) %% solve(V.hat) %% X) %% t(X) %% solve(V.hat)
tr = function(x) sum(diag(x))
se_asym = function(P, ZZ){
I.REML = matrix(NA, 2, 2)
I.REML[1, 1] = tr(P %% ZZ %% P %% ZZ)
I.REML[1, 2] = I.REML[2, 1] = tr(P %% P %% ZZ)
I.REML[2, 2] = tr(P %% P)
I.inv.REML = solve(0.5 I.REML)
sqrt(diag(I.inv.REML))
}
se.reml = se_asym(P.hat, ZZ)
(sigma2REML.hat = se.reml[1])
(tau2REML.hat = se.reml[2])
#e Bootstrap Method for REML
library(mvtnorm)
se_boot = function(data, beta, V, method, B = 100){
var.boot = matrix(NA, 2, B)
Xb = X %% beta
for(i in 1:B) {
data$$y = t(rmvnorm(1, Xb, V))
mod = fit_lmer(data, method)
var.boot[1, i] = mod$$sigma2.hat
var.boot[2, i] = mod$$tau2.hat
}
# based on bootstrap formula, we can use sd()
c(sigma2 = sd(var.boot[1,]), tau2 = sd(var.boot[2,]))
}
# REML
se_boot(corps, REML$$beta.hat, V.hat, 'REML')
#g ML METHOD
##Variance Components
ML = fit_lmer(corps,method = "ML")
(ML$$sigma2.hat)
(ML$$tau2.hat)
##Fixed Effects
(ML$$beta.hat)
#s.e. of variance components
##ACM
V.hat.ml = ML$$sigma2.hatZZ+ML$$tau2.hat diag(1,n)
se.ml = se_asym(solve(V.hat.ml), ZZ)
(sigma2ML.se = se.ml[1])
(tau2ML.se = se.ml[2])
##BOOTSTRAP
se_boot(corps, ML$$beta.hat, V.hat.ml, 'ML')