HTML-slides
http://coherence-based-reasoning-and-rationality.de/amsterdam/multilevelReg.html
PDF-print
http://coherence-based-reasoning-and-rationality.de/amsterdam/slides.pdf
July 2016
HTML-slides
http://coherence-based-reasoning-and-rationality.de/amsterdam/multilevelReg.html
PDF-print
http://coherence-based-reasoning-and-rationality.de/amsterdam/slides.pdf
Luke. (2004). Multilevel Modeling. (p. 5)
Moerbeek & Teerenstra. (2016). Power analysis of trials with multilevel data. (p. 7)
Examples:
Multilevel models…
\[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)):
If we have nested data, why don't we use multilevel models?
LEMMA: University of Bristol; http://www.bristol.ac.uk/cmm/learning/online-course/index.html
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
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\]
You can find the code here.
## $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)
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) \]
A-priori power analysis
\(\alpha\)-Level
\(\beta\)-Level (resp. power \(1-\beta\))
effect size
n = ???
Example