Week Learning Objectives

• Describe, conceptually, what the likelihood function and maximum likelihood estimation are

• Describe the differences between maximum likelihood and restricted maximum likelihood

• Conduct statistical tests for fixed effects, and use the small-sample correction when needed

• Use the likelihood ratio test to test random slopes

• Estimate multilevel models with the Bayesian/Markov Chain Monte Carlo estimator in the brms package

Estimation

Regression: OLS

MLM: Maximum likelihood, Bayesian

The most commonly used methods in MLM are

maximum likelihood (ML) and restricted maximum likelihood (REML)

># Linear mixed model fit by REML ['lmerMod']
># Formula: Reaction ~ Days + (Days | Subject)
>#    Data: sleepstudy
># REML criterion at convergence: 1744
># Random effects:
>#  Groups   Name        Std.Dev. Corr
>#  Subject  (Intercept) 24.74        
>#           Days         5.92    0.07
>#  Residual             25.59        
># Number of obs: 180, groups:  Subject, 18
># Fixed Effects:
># (Intercept)         Days  
>#       251.4         10.5

For MLM

Find $\gamma$s, $\tau$s, and $\sigma$ that maximizes the likelihood function

$$\ell (\gamma ,\tau ,\sigma ;y)=-\frac{1}{2}\left\{\log |V\right(\tau ,\sigma )|+(y-X\gamma {)}^{\top }{V}^{-1}(\tau ,\sigma )(y-X\gamma )\}+K$$

Here's the log-likelihood function for the coefficient of meanses (see code in the provided Rmd):

Numerical Algorithms

># iteration: 1
>#     f(x) = 47022.519159
># iteration: 2
>#     f(x) = 47151.291766
># iteration: 3
>#     f(x) = 47039.480137
># iteration: 4
>#     f(x) = 46974.909593
># iteration: 5
>#     f(x) = 46990.872588
># iteration: 6
>#     f(x) = 46966.453125
># iteration: 7
>#     f(x) = 46961.719993
># iteration: 8
>#     f(x) = 46965.890703
># iteration: 9
>#     f(x) = 46961.367013
># iteration: 10
>#     f(x) = 46961.288830
># iteration: 11
>#     f(x) = 46961.298898
># iteration: 12
>#     f(x) = 46961.284848
># iteration: 13
>#     f(x) = 46961.285238
># iteration: 14
>#     f(x) = 46961.284845
># iteration: 15
>#     f(x) = 46961.284848
># iteration: 16
>#     f(x) = 46961.284845

ML vs. REML

REML has corrected degrees of freedom for the variance component estimates (like dividing by $N-1$ instead of by $N$ in estimating variance)

• REML is generally preferred in smaller samples

• The difference is small when the number of clusters is large

Technically speaking, REML only estimates the variance components¹

Other Estimation Methods

Generalized estimating equations (GEE)

• Robust to some misspecification and non-normality
• Maybe inefficient in small samples (i.e., with lower power)
• See Snijders & Bosker 12.2; the geepack R package

Markov Chain Monte Carlo (MCMC)/Bayesian

• Researchers set prior distributions for the parameters
  • Different from "empirical Bayes": Prior coming from the data
• Does not depend on normality of the sampling distributions
  • More stable in small samples with the use of priors
• Can handle complex models
• See Snijders & Bosker 12.1

Fixed effects $\left(\gamma \right)$

• Usually, the likelihood-based CI/likelihood-ratio (LRT; ${\chi }^2$) test is sufficient
  • Require ML (as fixed effects are not part of the likelihood function in REML)
• Small sample (10--50 clusters): Kenward-Roger approximation of degrees of freedom
• Non-normality: Residual bootstrap¹

Random effects $\left(\tau \right)$

• LRT (with $p$ values divided by 2)

Likelihood Ratio (Deviance) Test

$H_0:\gamma =0$

Likelihood ratio: $\frac{L(\gamma =0)}{L(\gamma =\hat{\gamma })}$

Deviance: $-2\times \log \left(\frac{L(\gamma =0)}{L(\gamma =\hat{\gamma })}\right)$
= $-2LL(\gamma =0)-[-2LL(\gamma =\hat{\gamma }\left)\right]$
= $Deviance{\mid }_{\gamma =0}-Deviance{\mid }_{\gamma =\hat{\gamma }}$

ML (instead of REML) should be used

Example

...
># Linear mixed model fit by maximum likelihood  ['lmerMod']
># Formula: mathach ~ (1 | id)
>#      AIC      BIC   logLik deviance df.resid 
>#    47122    47142   -23558    47116     7182 
...

...
># Linear mixed model fit by maximum likelihood  ['lmerMod']
># Formula: mathach ~ meanses + (1 | id)
>#      AIC      BIC   logLik deviance df.resid 
>#    46967    46995   -23480    46959     7181 
...

pchisq(47115.81 - 46959.11, df = 1, lower.tail = FALSE)

># [1] 5.95e-36

In lme4, you can also use

anova(m_lv2, ran_int)  # Automatically use ML

Problem of LRT in Small Samples

$F$ Test With Small-Sample Correction

It is based on the Wald test (not the LRT):

• $t=\hat{\gamma }/\widehat{se}\left(\hat{\gamma }\right)$,
• Or equivalently, the $F=t^2$ (for a one-parameter test)

The small-sample correction does two things:

• Adjust $\widehat{se}\left(\hat{\gamma }\right)$ as it tends to be underestimated in small samples
• Determine the critical value based on an $F$ distribution, with an approximate denominator degrees of freedom (ddf)

Kenward-Roger (1997) Correction

Generally performs well with < 50 clusters

# Wald
anova(m_contextual, ddf = "lme4")

># Analysis of Variance Table
>#         npar Sum Sq Mean Sq F value
># meanses    1    860     860    26.4
># ses        1   1874    1874    57.5

# K-R 
anova(m_contextual, ddf = "Kenward-Roger")

># Type III Analysis of Variance Table with Kenward-Roger's method
>#         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
># meanses    324     324     1    16    9.96  0.0063 ** 
># ses       1874    1874     1   669   57.53 1.1e-13 ***
># ---
># Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For meanses, the critical value (and the $p$ value) is determined based on an $F(1,15.51)$ distribution, which has a critical value of

qf(.95, df1 = 1, df2 = 15.51)

># [1] 4.52

LRT for Random Slopes

Example: LRT for ${\tau }_1^2$

...
># Formula: mathach ~ meanses + ses_cmc + (ses_cmc | id)
>#    Data: hsball
># REML criterion at convergence: 46558
...

...
># Formula: mathach ~ meanses + ses_cmc + (1 | id)
>#    Data: hsball
># REML criterion at convergence: 46569
...

pchisq(10.92681, df = 2, lower.tail = FALSE)

># [1] 0.00424

Need to divide by 2

Benefits

• Better small sample properties

• Less likely to have convergence issues

• CIs available for any quantities ($R^2$, predictions, etc)

• Support complex models not possible with lme4

  • E.g., 2+ random slopes

• Is getting increasingly popular in recent research

Downsides

• Computationally more intensive

  • But dealing with convergence issues with ML/REML may end up taking more time

• Need to learn some new terminologies

Bayes's Theorem

Posterior $\propto$ Likelihood $\times$ Prior

• So, in addition to the likelihood function (as in ML/REML), we need prior information

• Prior: Belief about a parameter before looking at the data

• For this course, we use default priors set by the brms package

  • Note that the default priors may change when software gets updated, so keep track of the package version when you run analyses

A simulation-based approach to approximate the sampling distribution of fixed and random effects

• Useful for obtaining CIs
• Especially for statistics that are functions of fixed/random effects (e.g., $R^2$)

Parametric, Residual, and Cases bootstrap

In my own work,¹ the residual bootstrap was found to perform best, especially when data are not normally distributed and when the number of clusters is small

See R code for this week

What is Likelihood?

Learning the Parameter From the Sample

Assume that we have scores from 5 representative students

Likelihood

If we assume that $\mu =10$, how likely will we get 5 students with these scores?

># [1] 4.21e-08

If $\mu =13$

># [1] 5.98e-08

Compute the likelihood for a range of $\mu$ values

Curvature and Standard Errors