Applied Multilevel Analysis

class: center, middle, inverse, title-slide

# Applied Multilevel Analysis
## Session 7b: Model Assumptions
### Merlin Schaeffer<br> Department of Sociology
### 2022-03-21

---

# The outcome we want to explain

---
layout: false
class: clear
# A random intercept model .font60[Stata]

.push-left[
.bg-washed-green.b--transparent.ba.bw2.br3.shadow-5.ph4.mt2[
**Small recap:**
- What is a micro predictor?
- Which variables potentially correct for composition differences?
- What is a macro predictor?
- Why do the df vary for the different predictors?
- What is a cross-level interaction?
- What do we learn from this model substantially about what predicts anti-Muslim prejudice?
]]

.push-right[
.font60[

```
. quietly do "./../../assets/ESS7/12-Assumptions.do"
. mixed z_anti_muslim_ij agea dfegcf TV_not_news_ij hincfel c_gini_2011 c.z_edu
> yrs_ij##c.z_geothreat_j || cntry: z_eduyrs_ij, reml covariance(unstructured) 
> dfmethod(satterthwaite) dftable(pvalue) vsquish noheader

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -34761.715  
Iteration 1:   log restricted-likelihood = -34761.715

Computing standard errors:

Computing degrees of freedom:

----------------------------------------------------------------------
z_anti_muslim~j |      Coef.   Std. Err.           DF       t    P>|t|
----------------+-----------------------------------------------------
           agea |   .0045544   .0002896       27640.9    15.73   0.000
         dfegcf |   .2443825   .0080135       27764.7    30.50   0.000
 TV_not_news_ij |   .0309307   .0030103       27762.9    10.28   0.000
        hincfel |   .1132938   .0069974       26499.3    16.19   0.000
    c_gini_2011 |  -.0024174   .0149423          16.5    -0.16   0.873
    z_eduyrs_ij |  -.1797526   .0130006          18.8   -13.83   0.000
  z_geothreat_j |   .1827897   .0690602          17.4     2.65   0.017
  c.z_eduyrs_ij#|
c.z_geothreat_j |    .042974   .0133167          18.7     3.23   0.004
          _cons |  -1.040551    .433216          16.9    -2.40   0.028
----------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cntry: Unstructured          |
               var(z_eduy~j) |   .0026373   .0011292      .0011394     .006104
                  var(_cons) |   .0886524   .0304357      .0452334    .1737484
         cov(z_eduy~j,_cons) |   .0107446    .004829      .0012799    .0202093
-----------------------------+------------------------------------------------
               var(Residual) |   .7100294    .006029      .6983106     .721945
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 3022.26               Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. 
```
]]

---
# Model assumptions

.push-left[
Like OLS, multilevel models only predict **accurate conditional means**, if all random effects have an expected value of zero after all  predictors have been taken into account (e.g., `$\mathbb{E}(U_{0j}|x_{ij},z_{j})=0$`):
1. **No confounder bias** the slope estimates. .alert[Impossible to test, you will need to theorize carefully.]
2. There are no **outliers**.
  + Micro-level outliers are rare in multilevel studies, because the micro-level sample size tends to be large.
  + .alert[Macro-level outlier are a serious concern, particularly since the macro-level sample size tends to be small (Van der Meer, Te Grotenhuis, and Pelzer, 2010)!]
2. The estimated associations are **linear**.
]

.push-right[
REML **standard errors** for the slope estimates are only accurate if:
3. No **(auto-)correlation**: All cluster-driven (auto-)correlation is accurately captured by random effects.
  + .alert[Difficult to test, you will need to think carefully about: At which levels do your predictors vary?]
4. **homoscedasticity**:
  + `$\text{Var}(R_{ij}) = \sigma^2$`, `$\sigma_{j}^2= \sigma^2$`,
  + `$\text{Var}(U_{0j}) = \tau_{0}^2$`,
  + ...
5. **Normally distributed errors**:
  + Matters little for `$R_{ij}$` because micro-level samples tend to be large.
  + Also matters little for `$U_{0j}$` and `$U_{1j}$`, although macro-level samples tend to be small .
]

---
layout: true
# No confounder bias

.push-left[
- "A variable `$X$` us a *cause* of a variable `$Y$` if `$Y$` in any way relies on `$X$` for its value. [...] think of causation as a form of listening; `$X$` is a cause of `$Y$` if `$Y$` listens to `$X$` and decisdes its value in response to what it hears" .
- A confounder is a common cause of both `$X_{ij}$` (or `$Z_{j}$`) and `$Y_{ij}$`. It makes it appear as if `$Y_{ij}$` listens to `$X_{ij}$` (or `$Z_{j}$`), but in reality both listen to the confounder.
- In the absence of an experiment (i.e. randmized controlled trial), you need to control for all confounders to get proper causal estimate.
]

---

.push-right[

<blockquote class="twitter-tweet"><p lang="en" dir="ltr">It&#39;s going in the teaching slides <a href="https://t.co/bspSFNxgam">pic.twitter.com/bspSFNxgam</a></p>&mdash; Pietro Biroli (@pietrobiroli) <a href="https://twitter.com/pietrobiroli/status/1309805942180118528?ref_src=twsrc%5Etfw">September 26, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

]

---

.push-right[
We can compare superficial contact to close friendships, with the latter being hypothesized to be the most effective. But: 
- We cannot model the longitudinal process.
- Supportive conditions are not necessarily reflected.
- .alert[Across all measures, selection bias is a problem:]
<img src="7b-Model-Assumptions_files/figure-html/DAG1-1.png" width="90%" style="display: block; margin: auto;" />
]

---
layout: false
# Post-estimation .font70[for model diagnostics]

.panelset[
.panel[.panel-name[Stata]

```stata
// Estimate the model
mixed z_anti_muslim_ij agea dfegcf TV_not_news_ij hincfel c_gini_2011 ///
  c.z_eduyrs_ij##c.z_geothreat_j || cntry: z_eduyrs_ij, reml covariance(unstructured)

// Calculate prediction based on fixed random effects
predict y_hat_ij, fitted

// Calculate R_ij, i.e., y_ij - y_hat_ij
predict r_ij, residuals

// Predict U_0j and U_1j, 
// "*" ensures that both are predicted and called eblup1 and eblup2
predict eblup*, reffects
```
]
.panel[.panel-name[R]

```r
library(modelr) # Useful tools for postestimation.
ESS <- ESS %>% # Make sure there are no unused factor levelsleft (i.e. here in cntry).
  select(z_anti_muslim_ij, agea, TV_not_news_ij, hincfel, z_eduyrs_ij, z_geothreat_j, cntry) %>%
  drop_na() %>%
  mutate(cntry = fct_drop(cntry))

# Estimate the model
ris_mod <- lmer(z_anti_muslim_ij ~ agea + TV_not_news_ij + hincfel + z_eduyrs_ij*z_geothreat_j + 
                  (1 + z_eduyrs_ij | cntry), data = ESS)

# Calculate prediction based on fixed and random effects
ESS <- add_predictions(model = ris_mod, data = ESS)

# Calculate R_ij (i.e. y_ij - y_hat_ij)
ESS <- add_residuals(model = ris_mod, data = ESS)

# Predict U_0j and U_1j
Blups <- ranef(ris_mod) %>% as_tibble()
```
]]

---
# `$\sigma_{j}^2= \sigma^2$` & `$R_{ij} ~ \mathcal{N}(0, \sigma^2)$`

.panelset[
.panel[.panel-name[Stata]
.left-column[

```stata
// Plot R_ij by cluster
hist r_ij, by(cntry) xline(0)
```

- Overall, `$R_{ij}$` seems to look pretty similar across the countries.
- But the sparks suggest that our model does not seem to work equally well in all countries.
  + Missing random slopes?
  + Missing micro-level predictors?
]

.right-column[

<img src="media/histogram3.svg" width="80%" style="display: block; margin: auto;" />
]]
.panel[.panel-name[R]

.left-column[
- Overall, `$R_{ij}$` seems to look pretty similar across the countries.
- But our model does not seem to work well in Czech Republic and Hungary.
  + Missing random slopes?
  + Missing micro-level predictors?
]

.right-column[

```r
# Plot R_ij by cluster
ggplot(ESS, aes(resid, fill = cntry, colour = cntry)) +
  geom_density(alpha = 0.1) + theme_minimal()
```

<img src="7b-Model-Assumptions_files/figure-html/unnamed-chunk-7-1.png" width="80%" style="display: block; margin: auto;" />
]]]

---
# Linearity & `$\text{Var}(R_{ij}) = \sigma^2$`

.panelset[
.panel[.panel-name[Stata]
.push-left[
.center[**Plot Residuals against predictions**]

```stata
// Plot r_ij by y_hat_ij including a non-parametric smooth
twoway scatter r_ij y_hat_ij, mcolor(%1) ///
  || lowess r_ij y_hat_ij
```

<img src="media/linearity1.svg" width="80%" style="display: block; margin: auto;" />
]

.push-right[
.center[**Also test against the single categorical predictors**]

```stata
// Plot r_ij by single predictors, too
twoway scatter r_ij agea, mcolor(%1) //
  || lowess r_ij agea
```

<img src="media/linearity2.svg" width="80%" style="display: block; margin: auto;" />
]]

.panel[.panel-name[R]
.push-left[
.center[**Plot Residuals against predictions**]

```r
library(ggpointdensity)
# Plot r_ij by y_hat_ij including a non-parametric smooth
ggplot(ESS, aes(y = resid, x = pred)) +
  geom_pointdensity() +
  geom_smooth(se = FALSE, color = "red") +
  theme_minimal()
```

<img src="7b-Model-Assumptions_files/figure-html/unnamed-chunk-13-1.png" width="80%" style="display: block; margin: auto;" />
]

.push-right[
.center[**Also test against the single categorical predictors**]

```r
# Plot r_ij by single predictors, too
ggplot(ESS, aes(y = resid, x = agea)) +
  geom_pointdensity() +
  geom_smooth(se = FALSE, color = "red") +
  theme_minimal()
```

<img src="7b-Model-Assumptions_files/figure-html/unnamed-chunk-14-1.png" width="80%" style="display: block; margin: auto;" />
]]]

---
class: clear
# Visually inspect normality &  outliers among BLUP

.push-left[

```stata
// Plot the predicted U_0j (random intercept)
graph dot (mean) eblup1, over(cntry, sort(eblup1)) ///
  yline(0) vertical ytitle("U{sub:0j}")
```

<img src="media/blup1.svg" width="100%" style="display: block; margin: auto;" />
]

.push-right[

```stata
// Plot the predicted U_1j (random slope)
graph dot (mean) eblup2, over(cntry, sort(eblup2)) ///
  yline(0) vertical ytitle("U{sub:0j}")
```

<img src="media/blup2.svg" width="100%" style="display: block; margin: auto;" />
]

---
# Systematic macro outlier analysis

.push-left[
We tend to use two statistics for systematic outlier analysis: Cook's D and DFBETA. Cook's D evaluates the overall model. But most of the time, we are interested in a specific slope estimate. Thus we will focus on DFBETAs:

`$$\text{DFBETA}_{jk} = \frac{\hat{\beta}_{k} - \hat{\beta}_{(-j)k}}{\text{SE}(\hat{\beta}_{(-j)k})},$$`
where `$\hat{\beta}_{k}$` is the full-sample slope estimate for predictor `$k$`, and `$\hat{\beta}_{(-j)k}$` is the slope estimate for predictor `$k$` excluding cluster `$j$`. The denominator standardizes for uncertainty because with larger uncertainty, considerable differences are expected.
]
.push-right[
The accepted cut-off for identifying outliers is `$\text{DFBETA} > \frac{2}{\sqrt{N}}$`, where `$N$` is the number of clusters.

`$\rightarrow$` In our case the cut-off is: `$\frac{2}{\sqrt{20}} = 0.447$`.
]

---
# Systematic macro outlier analysis

.panelset[
.panel[.panel-name[Stata code]
.left-column[
The `mlt` ado contains the command `mltcooksd`, which allows you to perform macro outlier analysis in Stata. To install it, type `ssc install mlt`.

- Note, that `mltcooksd` requires you to specify your model via the `xtmixed` command.
- Also, `mltcooksd` cannot handle `##` interaction terms. Thus you will need to generate the interaction variable by hand.
- Finally, `mltcooksd` takes quite some time because it estimates N models.
]

.right-column[
.font70[

```
. quietly do "./../../assets/ESS7/12-Assumptions.do"
. gen cli = c.z_eduyrs_ij*z_geothreat
(254 missing values generated)

. quietly xtmixed z_anti_muslim_ij agea dfegcf TV_not_news hincfel c.z_eduyrs_i
> j z_geothreat cli || cntry: c.z_eduyrs_ij, reml cov(unstructured)

. mltcooksd, graph
Level 2 variable is cntry 
 
Calculating DFBETAs for the fixed effects of 
 agea dfegcf TV_not_news_ij hincfel z_eduyrs_ij z_geothreat_j cli _cons 
 
Cutoff value for DFBETAs is
0.4472
Cutoff value for Cook's D is
0.2000
 
Level-two units with Cook's D above the cut off value:

+-----------------+
  | L2ID     CooksD |
  |-----------------|
  |   SI   .5586839 |
  |   EE   .4662971 |
  |   HU   .3918715 |
  |   CZ   .3610041 |
  |   AT   .2778937 |
  |   SE   .2777313 |
  |   DE   .2551393 |
  |   GB   .2357559 |
  +-----------------+
   Legend: CooksD   = overall Cook's D
 
 
Level-two units with DFBETAs above cut off value:

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   SI |   0.2827 |     0.2357 |    -0.3090 |     0.4786 |    -0.4446 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.1455       |        0.5644       |        -0.0615       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   EE |   2.2442 |    -0.7551 |     0.1211 |     0.2622 |     0.1127 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.2171       |       -0.1330       |        -0.1309       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   HU |  -0.5861 |     0.2904 |     0.2635 |    -0.7298 |     0.2398 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.5142       |        0.5850       |         0.4178       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   CZ |  -0.5973 |     0.0406 |    -0.7703 |    -0.7519 |     0.2145 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.4726       |        0.4808       |         0.7150       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   AT |  -0.8592 |     0.6626 |     1.2596 |     0.1987 |    -0.2186 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.1792       |       -0.3105       |        -0.2878       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   SE |   0.0210 |    -0.5426 |    -0.1192 |     0.2157 |    -0.2285 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.2632       |        0.0694       |        -0.4034       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   DE |  -0.1469 |     0.3062 |     0.7211 |     0.7193 |    -0.2348 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.7287       |       -0.4293       |        -0.8562       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   GB |   0.2869 |    -1.1386 |    -0.2764 |     0.4578 |    -0.2688 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.2815       |       -0.2786       |        -0.0380       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   ES |  -0.4525 |     0.6181 |     0.2115 |    -0.3444 |     0.1010 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.0183       |        0.0304       |        -0.0712       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   IE |  -0.3597 |     0.4846 |     0.4828 |     0.3581 |    -0.0126 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.0779       |        0.0378       |        -0.0454       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   LT |  -0.2633 |    -0.2458 |    -0.6066 |    -0.5911 |     0.0905 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.1116       |        0.0680       |         0.5281       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   NL |  -0.8510 |    -0.6864 |     0.0543 |     0.1715 |     0.0978 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.0701       |       -0.2483       |         0.3778       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   FI |   0.9920 |     0.1890 |    -0.1008 |    -0.1363 |     0.1429 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.0337       |       -0.0428       |        -0.1442       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   NO |   0.5085 |    -0.1035 |    -0.2153 |    -0.3758 |    -0.1131 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |           0.0908       |        0.1252       |        -0.0487       |
  +---------------------------------------------------------------------+

+---------------------------------------------------------------------+
  | L2ID | DFB_agea | DFB_dfegcf | DFB_TV_n~j | DFB_hinc~l | DFB_z_ed~j |
  |   BE |  -0.1991 |     0.6008 |    -0.1125 |     0.0609 |    -0.1148 |
  |---------------------------------------------------------------------|
  |       DFB_z_ge~j       |       DFB_cli       |       DFB_cons       |
  |          -0.0326       |       -0.0573       |        -0.1921       |
  +---------------------------------------------------------------------+

. graph export "./media/DFBETA.svg", replace
(file ./media/DFBETA.svg written in SVG format)

. 
```
]]]

.panel[.panel-name[Stata result]
<img src="media/DFBETA.svg" width="55%" style="display: block; margin: auto;" />
]

.panel[.panel-name[R code]
.left-column[
The `influence.ME` packages allows you to perform macro outlier analysis in R. `influence()` will estimate N models, afterwards `dfbetas()` calculates the DFBETAS.
]

.right-column[
.font80[

```r
library(influence.ME)
# Estimate -j models
one_out_models <- influence(model = ris_mod, group = "cntry")
# Calculate the DFBETAs and turn into tibble
(DFBETAs <- dfbetas(one_out_models) %>% 
    data.frame() %>% 
    rownames_to_column() %>% 
    as_tibble() %>%
    select(rowname, z_eduyrs_ij, z_geothreat_j, 
           z_eduyrs_ij.z_geothreat_j))
# # A tibble: 20 × 4
#    rowname z_eduyrs_ij z_geothreat_j z_eduyrs_ij.z_geothreat_j
#    <chr>         <dbl>         <dbl>                     <dbl>
#  1 AT         -0.327        -0.171                     -0.362 
#  2 BE         -0.116        -0.0272                    -0.0477
#  3 CH         -0.00573       0.00731                   -0.0194
#  4 CZ          0.201         0.483                      0.452 
#  5 DE         -0.225        -0.667                     -0.418 
#  6 DK          0.299         0.0178                    -0.116 
#  7 EE          0.133         0.170                     -0.123 
#  8 ES          0.0910       -0.0171                     0.0417
#  9 FI          0.124        -0.0427                    -0.0374
# 10 FR          0.00699       0.0604                     0.0144
# 11 GB         -0.206        -0.326                     -0.281 
# 12 HU          0.293         0.510                      0.734 
# 13 IE         -0.0471       -0.0793                     0.0455
# 14 LT          0.0402        0.106                      0.0509
# 15 NL          0.131        -0.0772                    -0.261 
# 16 NO         -0.0595        0.0807                     0.0753
# 17 PL          0.116         0.0485                     0.0457
# 18 PT          0.219        -0.191                     -0.291 
# 19 SE         -0.218         0.266                      0.0829
# 20 SI         -0.406        -0.136                      0.483
```
]
]]

.panel[.panel-name[R result]
.push-left[

```r
# Reshape to long
DFBETAs <- pivot_longer(
  DFBETAs, 
  cols = c("z_eduyrs_ij", "z_geothreat_j",
           "z_eduyrs_ij.z_geothreat_j"), 
  names_to = "Variable")
```

```r
ggplot(data = DFBETAs, 
       aes(y = value, x = Variable, label = rowname)) +
  geom_violin(color = "gray") +
  geom_text(position=position_jitter(width=0.2,height=0)) +
  geom_hline(yintercept = 2/sqrt(20), color = "#901A1E") +
  geom_hline(yintercept = -2/sqrt(20), color = "#901A1E") +
  theme_minimal() +
  labs(y = "DFBETA")
```
]
.push-right[
<img src="7b-Model-Assumptions_files/figure-html/unnamed-chunk-24-1.png" width="100%" style="display: block; margin: auto;" />
]
]]

---
class: clear
# Re-estimate model without outliers .font60[Report both results]
.panelset[
.panel[.panel-name[Stata]
.font70[

```
. quietly do "./../../assets/ESS7/12-Assumptions.do"
. drop if cntry == "HU" | cntry == "DE" | cntry == "CZ" | cntry == "SI"
(6,986 observations deleted)

. mixed z_anti_muslim_ij agea dfegcf TV_not_news hincfel c.z_eduyrs_ij##c.z_geo
> threat || cntry: z_eduyrs_ij, reml covariance(unstructured) dfmethod(satterth
> waite) dftable(pvalue) vsquish noheader

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -27384.286  
Iteration 1:   log restricted-likelihood = -27384.285

Computing standard errors:

Computing degrees of freedom:

----------------------------------------------------------------------
z_anti_muslim~j |      Coef.   Std. Err.           DF       t    P>|t|
----------------+-----------------------------------------------------
           agea |   .0049224    .000335       20894.5    14.69   0.000
         dfegcf |   .2359193   .0090343       21510.5    26.11   0.000
 TV_not_news_ij |   .0313767    .003516       21458.6     8.92   0.000
        hincfel |   .1161512   .0080579       19309.3    14.41   0.000
    z_eduyrs_ij |  -.1788924   .0108054          15.4   -16.56   0.000
  z_geothreat_j |   .1619543    .072061          14.0     2.25   0.041
  c.z_eduyrs_ij#|
c.z_geothreat_j |   .0220288   .0125698          15.1     1.75   0.100
          _cons |   -1.13353   .0685701          21.4   -16.53   0.000
----------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cntry: Unstructured          |
               var(z_eduy~j) |   .0010738   .0006376      .0003353    .0034387
                  var(_cons) |   .0575167   .0219738      .0272019    .1216153
         cov(z_eduy~j,_cons) |   .0060299   .0031594     -.0001623    .0122221
-----------------------------+------------------------------------------------
               var(Residual) |   .7388641   .0071239      .7250327    .7529594
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 1278.41               Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. 
```
]]

.panel[.panel-name[R]
.left-column[
.font80[

```r
ESS_NoOutlier <- ESS %>% # Drop outliers
  filter(cntry != "HU" & cntry != "DE" & 
           cntry != "CZ" & cntry != "SI")

# Estimate model without outliers
ris_mod_noOutlier <- lmer(
  z_anti_muslim_ij ~ agea + TV_not_news_ij + 
    hincfel + z_eduyrs_ij*z_geothreat_j + 
    (1 + z_eduyrs_ij | cntry), 
  data = ESS_NoOutlier) %>%
  # Add approx. df and accurate p-values
  lmer_t(., method = "Satterthwaite")
```
]]

.right-column[
.font70[

```r
summary(ris_mod_noOutlier)
# Linear mixed model fit by REML ['lmerMod']
#    t-tests use the Satterthwaite method.
# Formula: z_anti_muslim_ij ~ agea + TV_not_news_ij + hincfel + z_eduyrs_ij *  
#     z_geothreat_j + (1 + z_eduyrs_ij | cntry)
#    Data: ESS_NoOutlier
# 
# REML criterion at convergence: 55521
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.214 -0.711  0.006  0.758  2.832 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr
#  cntry    (Intercept) 0.0708   0.2661       
#           z_eduyrs_ij 0.0015   0.0387   0.72
#  Residual             0.7619   0.8729       
# Number of obs: 21585, groups:  cntry, 16
# 
# Coefficients:
#                           Estimate Std.Err t value    Df    Lower  Upper  Pr(>|t|)    
# (Intercept)               -0.648   0.072      -8.94    17.7 -0.800 -0.495 5.69e-08 ***
# agea                       0.007   0.000      19.79 21135.6  0.006  0.007 < 2e-16  ***
# TV_not_news_ij             0.037   0.004      10.29 21535.2  0.030  0.044 < 2e-16  ***
# hincfel                    0.109   0.008      13.33 20267.9  0.093  0.125 < 2e-16  ***
# z_eduyrs_ij               -0.200   0.012     -16.58    15.0 -0.226 -0.175 4.83e-11 ***
# z_geothreat_j              0.175   0.080       2.19    14.0  0.003  0.346 0.0462   *  
# z_eduyrs_ij:z_geothreat_j  0.019   0.014       1.36    14.9 -0.011  0.049 0.1926      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]]]]

---
layout: false
# References

.font80[
Van der Meer, T., M. Te Grotenhuis, and B. Pelzer (2010). "Influential Cases in Multilevel Modeling: A
Methodological Comment". In: _American Sociological Review_ 75.1, pp. 173-178.
]