# MLE and REML Estimates

- 7 mins### The Restricted Maximum Likelihood Estimates (REML)

#### Estimate ,, and the ﬁxed eﬀects 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$.

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 by replacing the variance components $\theta$ by the REML estimators of variance components .
- We call

the REML estimator of V.

- Then we plug in into

- So we obtain the variance of REML estimates of fixed effects.

- Last we take the diagonals of and the square root of that are the standard errors for the REML estimate of fixed effects. The ’s standard error is the j
^{th}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 .
- Let REML estimates of variance components be

.

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)

- Simulate a sample from this multivariate normal distribution
- Fit a Linear Mixed Model to that sample to find REML estimates of .
- Let be the REML estimates of respectively from the sample.
- The bootstrap estimates are

where

is the mean of : b=1,2,3,…,B

is the mean of : 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 are much more higher than that of for both methods.
- The standard errors I obtained for by two methods are more similar compared to the s.e. I obtained for by two methods.
- That is because of different degree of freedom of . The standard error of is much more smaller since the degree of freedom is 99. The degree of freedom of 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 |

#### ML estimates of fixed effects

#### ACM method for standard errors of Variance Component

- Replace the with in . 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
_{REML}and_{REML}with_{ML}and_{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')
```