MLE and REML Estimates

Thursday. January 28, 2016 - 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

$\hat{\beta}_{\text{REML}} = (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

$(X'\hat{V}^{-1}X)^{\frac{1}{2}}(\hat{\beta} - \beta) \dot{\sim} N(0, I_{p})$ $V = V(\theta)=\sigma^{2} Z Z′+\tau^{2} I_{n}$

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

$\hat{V}_{\text{REML}} = V(\hat{\theta}_{\text{REML}})$

the REML estimator of V.

Then we plug in $\hat{V}_{\text{REML}}$ into

$S = Var(\hat{\beta}) = (X'\hat{V}^{-1}X)^{-1}$

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

$\hat{S} = Var(\hat{\beta}) = (X'\hat{V}^{-1}X)^{-1}$

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 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 $\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

$\text{s.e.}(\hat{\sigma}^{2}_{\text{REML}}) = \sqrt{a^{\ast}}$ $\text{s.e.}(\hat{\tau}^{2}_{\text{REML}}) = \sqrt{d^{\ast}}$

Parametric Bootstrap (PB) method

Also let REML estimates of variance components be

$\hat{\theta}_{\text{REML}} = (\hat{\sigma}^{2}_{\text{REML}},\hat{\tau}^{2}_{\text{REML}})$

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 be the REML estimates of respectively from the sample.
- The bootstrap estimates are

$\hat{\sigma}^{2}_{\text{REML}} = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}(\hat{\theta}^{(b)}_{\sigma^{2}}-\bar{\theta}_{\sigma^{2}})^{2}}$ $\hat{\tau}^{2}_{\text{REML}} = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}(\hat{\theta}^{(b)}_{\tau^{2}}-\bar{\theta}_{\tau^{2}})^{2}}$

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

$\text{s.e.}(\hat{\sigma}^{2}_{\text{REML}}) = 9511.913$ $\text{s.e.}(\hat{\tau}^{2}_{\text{REML}}) = 1181.954$

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

$\text{s.e.}(\hat{\sigma}^{2}_{\text{REML}}) = 9197.191$ $\text{s.e.}(\hat{\tau}^{2}_{\text{REML}}) = 1184.55$

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

$\text{s.e.}(\hat{\sigma}^{2}_{\text{ML}}) = \sqrt{a^{\ast}}$ $\text{s.e.}(\hat{\tau}^{2}_{\text{ML}}) = \sqrt{d^{\ast}}$

The standard errors for the ML estimates of the variance components by ACM method are

$\text{s.e.}(\hat{\sigma}^{2}_{\text{ML}}) = 8125.975$ $\text{s.e.}(\hat{\tau}^{2}_{\text{ML}}) = 1146.673$

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

$\text{s.e.}(\hat{\sigma}^{2}_{\text{ML}}) = 8081.138$ $\text{s.e.}(\hat{\tau}^{2}_{\text{ML}}) = 1083.170$

  #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')

Ying Zhang

A statistician who gets lost in analysis.