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

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

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

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

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

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

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

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

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

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

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

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

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

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?

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}}}\]

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})\]

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?

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

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

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)}\]

dotplot(ranef(modMLR, condVar= T))

## $group

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

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}\]

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

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

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

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

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

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

## $subject

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

## $subject

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

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)

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

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

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

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

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

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) \]

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

Example