MLE and REML Estimates

- 7 mins

The Restricted Maximum Likelihood Estimates (REML)

Estimate $\sigma^2$,$\tau^2$, and the ﬁxed eﬀects $\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 ﬁxed eﬀects

• 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

• 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)
1. Simulate a sample from this multivariate normal distribution $N(X\hat{\beta},\hat{V(\theta)})$
2. Fit a Linear Mixed Model to that sample to find REML estimates of $\theta = (\sigma^{2}, \tau^{2})$.
3. 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

Intercept year1 AGE SHARES REV INC
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$$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) (tau2REML.hat = se.reml) #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)
(tau2ML.se = se.ml)

##BOOTSTRAP
se_boot(corps, MLbeta.hat, V.hat.ml, 'ML')