July 2016

Slides

Multilevel-modeling

Luke. (2004). Multilevel Modeling. (p. 5)

Multilevel-modeling

Moerbeek & Teerenstra. (2016). Power analysis of trials with multilevel data. (p. 7)

Multilevel-modeling

Examples:

  • Games Choices in dictator game, participant (sex, age) (level1), cultural background of participant (level 2)
  • Educational Psychology Students (level 1), classes (level 2), schools (Level 3), school districts (level 4)
  • Clinical Psychology Patients (level 1), therapists (level 2)
  • Metanalysis Effects, type of independent variable (level 1): Impact of type of fear induction (subliminal priming, supraliminal priming, …) on risky choice

Multilevel models…

  • …can account for this hierarchy.
  • …allow for testing interesting (context-)questions (e.g., relation between variables of different levels; relative importance of differents levels/contexts).

Intercept + Error

\[Y_{i} = \beta_{0} + e_{i}\]

Intercept + Error

\[Y_{i} = \beta_{0} + e_{i} \rightarrow Y_{i} = \beta_{0} + 0 \times X_{i} + e_{i}\]

Intercept + Error

\[Y_{i} = \beta_{0} + e_{i} \rightarrow Y_{i} = \beta_{0} + 0 \times X_{i} + e_{i}\]

Intercept + Error

Residuals

\[e \sim N(0,\sigma)\]

Intercept + Error

Residuals

Intercept + Error

\[Y_{i} = \beta_{0} + e_{i}\]

Intercept + Error

\[Y_{i} = \beta_{0j} + e_{i}\]

Fixed and Random Intercept + Error

\[Y_{ij} = {\color{DarkMagenta}\beta_{0j}} + e_{i} = {\color{DarkMagenta}\beta_{00} + u_{0j}} + e_{ij}\]

Fixed and Random Intercept + Error

\[Y_{ij} = \beta_{00} + {\color{red} u_{0j}} + e_{ij}\]

Fixed and Random Intercept + Error

\[Y_{ij} = \beta_{00} + {\color{red} u_{0j}} + {\color{blue} e_{ij}}\]

Fixed and Random Intercept + Error

\[e \sim N(0,\sigma)\]

Fixed and Random Intercept + Error

Fixed and Random Intercept + Error

Assumptions

\[ Y_{ij} = \beta_{00} + {\color{red}u_{0j}} + e_{ij}\]

Distribution of residuals:

\[e \sim N(0,\sigma_e)\]

\[ {\color{red}u_0} \sim N(0,\sigma_0)\]

Covariation of residuals:

\[cov(e_{ij},u_{j}) = 0\]

Variance Partitioning Coefficient

Random Intercept Model

\[\rho = \frac{{\color{blue}\sigma_u^2}}{{\color{red}\sigma_u^2 + \sigma_e^2}}\]

The perentage of variance of the criterion explained by the groups is the estimated variance of random intercepts divided by the total variance.

\(\rho\) can vary between 0 and 1.

Variance Partitioning Coefficient

Alternative?

Linear regression with fixed effects for groups

Fixed effects model (i.e., dummy coding and continuous variable(s)):

  • Many groups means many free parameters.
  • For groups with few participants, estimates of group effects may be unreliable.
  • No inference beyond the observed groups are possible.
  • Due to linear dependency, the effect of group variables cannot be estimated when all dummies for groups are included.
  • Statistical power is higher.
  • Having missing values does not mean you have to exclude participants.
  • Heteroscedasticity is no problem.

If we have nested data, why don't we use multilevel models?

Residuals and Estimating Random Intercepts

Total variation left after accounting for fixed effect

\[ \text{resid}_{ij} = {\color{red} u_{0j}} + {\color{blue} e_{ij}}\]

\[ \sigma_{ \text{resid}_{ij}} = { \color{red} \sigma_{u_{0j}}} + {\color{blue} \sigma_{e_{ij}}}\]

Residuals and Estimating Random Intercepts

The mean residual is the difference between the the fixed intercept and the mean of group j:

\[{\color{Black}\bar{r}_{0j}} = {\color{red} \bar{y_j}} - {\color{blue} \beta_0}\]

The estimated residual \(u_j\) for group \(j\) is a weighted mean residual:

\[u_{0j} = {\color{green} k} \times \bar{r}_{0j}\]

Shrinkage factor:

\[ {\color{green}k} = \frac{\sigma^2_u}{\sigma^2_u+({\color{DarkOrchid}\sigma^2_e} / {\color{orange}n_j})} \]

Consequence: Few cases with high error variance lead to high shrinkage (i.e., estimates regressing to the fixed intercept).

Calculating individual residuals is easy:

\[e_{ij} = \hat{y}_{ij} - (\beta_0 + u_{0j})\]

Single Linear Regressions

Comparing Intercepts, 25 groups with 100 participants each

Single Linear Regressions

Comparing Intercepts, 25 groups with 100 participants each

Questions

  • When we decrease the ratio \(\sigma_u / \sigma_e\), do random effects change?
  • When we decrease the number of participants in each group, do random effects change?

Single Linear Regressions, \(n_j = 100\)

Single Linear Regressions, \(n_j = 10\)

Significance Tests

Significance Testing of Fixed and Random Effects

  • T-test with corrected degrees of freedom (e.g., Satterthwaite)
  • Likelihood Ratio Test
  • Profiling likelihood confidence intervals
  • (Bootstrapping)
  • (Bayesian Analysis)

Satterthwaite: Testing Fixed Effects

It is basically a t-test with corrected degrees of freadom:

\[t = \frac{b_i}{SE(b_i)}\]

coef(summary(modMLR))
##             Estimate Std. Error       df   t value  Pr(>|t|)
## (Intercept) 5.019234   13.71936 2.000642 0.3658503 0.7495396

Significance Testing of Random Effects

Estimate of the effect is defined as:

\[ \hat{u_j} = (\bar{y_j}-\beta_0) \times \frac{\sigma^2_u}{\sigma^2_u+(\sigma^2_e /n_j)} \]

Variance of the estimate is defined as:

\[var(\hat{u_j} - u_j) = \frac{\color{blue}\sigma^2_u \color{red}\sigma^2_e}{\color{orange}n_j\color{blue}\sigma^2_u+\color{red}\sigma^2_e}\]

95% Confidence intervals are defined as:

\[\hat{u_j} \pm 1.96 \times \sqrt{var(\hat{u_j} - u_j)}\]

Significance Testing of Random Effects

dotplot(ranef(modMLR, condVar= T))
## $group

Significance Testing of Fixed and Random Effects

Model Comparison

\[T_{LR} = 2 \times log\left(\frac{{\color{blue}L_{g}}}{{\color{red}L_{n}}}\right) = 2 \times (log(L_g) - log(L_n))\]

The likelihood of the general model is compared to the likelihood of the nested model. The test statistic \(T_{LR}\) ist (asymptotically) \(\chi^2\) distributed with degrees of freedom \(df = n_{p_g} - n_{p_{n}}\).

Number of Parameters

What is the addtional parameter in the more general model?

General model:

\[Y_{ij} = \beta_{00} + u_{0j} + e_{ij}\]

Nested model:

\[Y_{i} = \beta_{0} + e_{i}\]

Significance Testing of Fixed and Random Effects

Model Comparison

anova(modMLR,model)
## Data: NULL
## Models:
## ..1: Yij ~ 1
## object: Yij ~ 1 + (1 | group)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## ..1     2 799.49 804.49 -397.74   795.49                             
## object  3 580.21 587.71 -287.10   574.21 221.28      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example

In our example, \(T_{LR} = 2 \times (-287.1 + 397.74) = 221.28\) and \(df = 3 - 2 = 1\).

Estimating Confidence Intervals

Likelihood Profile

\[L_1(\theta) = max L(\theta,\hat{\delta})\]

\[2\times (log(L(\hat{\theta},\hat{\delta}))-log(L_1(\theta))) < \chi^2_{1-\alpha}(1)\]

confint(modMLR,level=.95)
##                  2.5 %    97.5 %
## .sig01       10.256922 57.290125
## .sigma        4.616839  6.219009
## (Intercept) -26.254657 36.293093

More Complex Regression Models

Random Intercept and Fixed Slope

\[ Y_{ij} = \beta_{00} + u_{0j} + \beta_{1} \times X_{ij} + e_{ij}\]

Random Intercept and Random Slope

\[ Y_{ij} = (\beta_{00} + u_{0j}) + (\beta_{01} + u_{1j}) \times X_{ij} + e_{ij}\]

Random Intercept and Random Slope

Assumptions

\[Y_{ij} = (\beta_{00} + u_{0j}) + (\beta_{01} + u_{1j}) \times X_{ij} + e_{ij}\]

Distribution of residuals:

\[e \sim N(0,\sigma_e)\]

\[u \sim MVN(0,\Sigma)\]

\[\Sigma \sim \left(\begin{array}{ccc} \sigma^2_{u_{0j}} & - \\ \color{red}{\sigma_{01}} & \sigma^2_{u_{1j}} \end{array} \right) \]

Random effects can covariate.

Covariation of residuals:

\[cov(e_{ij},\Sigma) = 0\]

Example

exampleMLRplot.R

You can find the code here.

Relation Between Random Effects

## $subject

Variance Partitioning Coefficient

Variance Partitioning Coefficient

Variance due to random effects

\[ var(u_{0j} + u_{1j}x_{ij}) = var(u_{0j}) + 2 \times x_{ij} \times cov(u_{0j},u_{1j}) + x_{ij}^2 \times var(u_{1j}) \]

\[ var(u_{0j} + u_{1j}x_{ij}) = \sigma^2_{u_0} + 2\sigma_{u01} x_{ij} + x^2_{ij} \sigma^2_{u_1} \]

Variance Partitioning Coefficient

\[\frac{{\color{green}\sigma^2_{u_0} + 2\sigma_{u01} x_{ij} + x^2_{ij} \sigma^2_{u_1}}}{{\color{green}\sigma^2_{u_0} + 2\sigma_{u01} x_{ij} + x^2_{ij} \sigma^2_{u_1}} + {\color{red}\sigma^2_e}} \]

Conclusion: The importance of context is dependent on the size of the predictor (i.e., heteroscedasticity is ok!)

Variance Partitioning Coefficient

Variance Partitioning Coefficient

Significance Testing of Estimates

## $subject

Significance Testing of Estimates

confint(modMLRFile,condVar=T)
##                  2.5 %     97.5 %
## .sig01       1.1754879  2.9133000
## .sig02       0.7341677  0.9851946
## .sig03       0.7432186  1.8898389
## .sigma       1.0153326  1.1981148
## (Intercept) 17.8491616 20.2427882
## pred         8.2388327  9.7990099

Model Diagnostics: Distribution of Residuals

residuals <- resid(modMLRFile)
summary(residuals)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.37400 -0.77870  0.03502  0.00000  0.73990  2.95400
hist(residuals)

Model Diagnostics: Q-Q Plots for Residuals

residuals <- resid(modMLRFile)
qqnorm(residuals)
qqline(residuals)

Model Diagnostics: Q-Q Plots for Random Intercepts

qqnorm(ranef(modMLRFile)$subject$`(Intercept)`, 
       main="Q-Q plot for the random intercept")
qqline(ranef(modMLRFile)$subject$`(Intercept)`)

Model Diagnostics: Q-Q Plots for Random Slopes

qqnorm(ranef(modMLRFile)$subject$`pred`, 
       main="Q-Q plot for the random slope")
qqline(ranef(modMLRFile)$subject$`pred`)

Model Diagnostics: Prediction versus Residual

residuals <- resid(modMLRFile)
xyplot(resid(modMLRFile) ~ fitted(modMLRFile), pch = 16)

Further Topics

  • (Adding level 2 explanatory variables)
  • Complex variance-covariance specifications
  • (Multilevel models with level > 2)
  • Power analysis in multilevel-models
  • (Multilevel version of logistic regression)

Complex Variance-Covariance Specifications

Uncorrelated residuals with equal variances

\[\left(\begin{array}{ccc} \sigma^2_{e} & - & - \\ 0 & \sigma^2_{e} &- \\ 0 & 0 & \sigma^2_{e}\\ \end{array} \right) \]

Most complex error-model

\[\left(\begin{array}{ccc} \sigma^2_{e_{11}} & - & - \\ \sigma_{e_{21}} & \sigma^2_{e_{22}} &- \\ \sigma_{e_{31}} & \sigma_{e_{23}} & \sigma^2_{e_{33}}\\ \end{array} \right) \]

Autoregressive residuals

\[\sigma_e^2 \times \left(\begin{array}{ccc} 1 & \rho & \rho^2 \\ \rho & 1 & \rho \\ \rho^2 & \rho & 1\\ \end{array} \right) \]

Power Analysis in Multilevel-Models

A-priori power analysis

\(\alpha\)-Level

\(\beta\)-Level (resp. power \(1-\beta\))

effect size

n = ???

Power Analysis in Multilevel-Models

MLPowSim

  • MLPowSim is based on a Monte-Carlo simulation
  • MLPowSim provides an interface for specifying your model
  • MLPowSim generates R-Code

Example