MLE and REML Estimates

- 7 mins

The Restricted Maximum Likelihood Estimates (REML)

Estimate ,, and the fixed effects via REML

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

the REML estimator of V.

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

.

And we define

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

Parametric Bootstrap (PB) method

from the original data.

  1. Simulate a sample from this multivariate normal distribution
  2. Fit a Linear Mixed Model to that sample to find REML estimates of .
  3. 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.

Compute the standard errors in REML estimates via PB method.

Compare the results of ACM and PB Comment on what we observe.

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

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

PB method for standard errors of Variance Component


  #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

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