MLE and REML Estimates

- 6 mins
Statistics Code

The Restricted Maximum Likelihood Estimates (REML)

Estimate \(\sigma^2\),\(\tau^2\), and the fixed effects \(\beta\) via REML

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 fixed effects

\[(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}\] \[\hat{V}_{\text{REML}} = V(\hat{\theta}_{\text{REML}})\]

the REML estimator of V.

\[S = Var(\hat{\beta}) = (X'\hat{V}^{-1}X)^{-1}\] \[\hat{S} = Var(\hat{\beta}) = (X'\hat{V}^{-1}X)^{-1}\]
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

\(\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}\]
$$ 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} $$
$$ I^{-1}(\hat{\theta}_{\text{REML}}) = \Big( \begin{matrix} a^{\ast} & b^{\ast} \\ c^{\ast} & d^{\ast} \end{matrix} \Big) $$
\[\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

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

from the original data.

  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
\[\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.

\[\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.

\[\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.

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

$$ 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} $$
$$ I^{-1}(\hat{\theta}_{\text{ML}}) = \Big( \begin{matrix} a^{\ast} & b^{\ast} \\ c^{\ast} & d^{\ast} \end{matrix} \Big) $$
\[\text{s.e.}(\hat{\sigma}^{2}_{\text{ML}}) = \sqrt{a^{\ast}}\] \[\text{s.e.}(\hat{\tau}^{2}_{\text{ML}}) = \sqrt{d^{\ast}}\] \[\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

\[\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')

Related Posts

Ying Zhang

Ying Zhang

A statistician who gets lost in analysis.

comments powered by Disqus
rss facebook twitter github youtube mail spotify instagram linkedin google google-plus pinterest medium vimeo stackoverflow reddit quora