1. Replicate Model 2 from Schlueter, Masso, and Davidov (2019).
library(tidyverse) # Data manipulation
library(haven) # Read Stata dta files
library(essurvey) # API to ESS

ESS <- import_rounds(rounds = 7, ess_email = YOUR_EMAIL) %>%
  recode_missings() %>%
  # redefine some variables as factors
  mutate(
    # Educational attainment
    edu_attain = case_when(
      edulvlb < 200 ~ 1,
      edulvlb > 199 & edulvlb < 300 ~ 2,
      edulvlb > 299 & edulvlb < 400 ~ 3,
      edulvlb > 399 & edulvlb < 500 ~ 4,
      edulvlb > 499 & edulvlb < 600 ~ 5,
      edulvlb > 599 & edulvlb < 700 ~ 6,
      edulvlb > 699 & edulvlb < 801 ~ 7,
      TRUE ~ as.numeric(NA)
    ),
    # Watching TV, but not news on TV
    TV_not_news = tvtot - tvpol,
    # Immigrant friends
    dfegcf = max(dfegcf, na.rm = TRUE) - dfegcf,
    brncntr = as_factor(brncntr), 
    ctzcntr = as_factor(ctzcntr),
    rlgblg = as_factor(rlgblg),
    rlgdnm = as_factor(rlgdnm),
    uempla = as_factor(uempla),
    rlgdnm = case_when( # Add atheists to the religion variable
      rlgblg == "No" ~ "Atheist",
      rlgdnm=="Not applicable"|rlgdnm=="Refusal"|rlgdnm=="No answer" ~ as.character(NA),
      TRUE ~ as.character(rlgdnm)
    ) %>% factor(),
    gndr = as_factor(gndr),
    gndr = case_when(
      gndr == "No answer" ~ as.character(NA),
      TRUE ~ as.character(gndr)
    ) %>% factor()
  ) %>% zap_labels() %>%
  # keep only those born and with citizenship in country of interview,
  # non-Muslim and drop Israel
  filter(brncntr == "Yes" & ctzcntr == "Yes" & rlgdnm != "Islamic" & cntry != "IL") %>%
  inner_join(., read_dta("./../../assets/ESS7/Schlueter.dta") %>%
               mutate(
                 z_muslim_j = (muslim_j - mean(muslim_j, na.rm= TRUE)) / sd(muslim_j, na.rm = TRUE)
               ), by = "cntry") %>%
  # Keep only the following variables
  select(cntry, almuslv, agea, gndr, edu_attain, TV_not_news,
         uempla, hincfel, rlgdnm, dfegcf, tvtot, tvpol, 
         muslim_j, z_muslim_j, mipex_j, relstate_j, claims_j) %>%
  drop_na() # Casewise deletion
library(lme4) # Multilevel modeling

# Model 2
Model_2 <- lmer(
  almuslv ~ gndr + agea + edu_attain + uempla + hincfel +
    dfegcf + rlgdnm + TV_not_news + muslim_j + mipex_j + 
    relstate_j + claims_j + (1 | cntry), data  = ESS
)

summary(Model_2)
# Linear mixed model fit by REML ['lmerMod']
# Formula: almuslv ~ gndr + agea + edu_attain + uempla + hincfel + dfegcf +  
#     rlgdnm + TV_not_news + muslim_j + mipex_j + relstate_j +      claims_j + (1 | cntry)
#    Data: ESS
# 
# REML criterion at convergence: 73799
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.557 -0.694  0.009  0.743  3.020 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  cntry    (Intercept) 0.0344   0.186   
#  Residual             0.6774   0.823   
# Number of obs: 30064, groups:  cntry, 20
# 
# Fixed effects:
#                                      Estimate Std. Error t value
# (Intercept)                          3.529811   0.250494   14.09
# gndrMale                            -0.001773   0.009618   -0.18
# agea                                 0.005507   0.000272   20.23
# edu_attain                          -0.090369   0.002825  -31.99
# uemplaMarked                        -0.026586   0.025022   -1.06
# hincfel                              0.101535   0.006686   15.19
# dfegcf                              -0.234470   0.007446  -31.49
# rlgdnmEastern Orthodox               0.000852   0.054319    0.02
# rlgdnmEastern religions             -0.092259   0.103340   -0.89
# rlgdnmJewish                        -0.022459   0.171995   -0.13
# rlgdnmOther Christian denomination  -0.077413   0.050117   -1.54
# rlgdnmOther non-Christian religions -0.221022   0.100889   -2.19
# rlgdnmProtestant                     0.032939   0.015282    2.16
# rlgdnmRoman Catholic                 0.057076   0.013169    4.33
# TV_not_news                          0.028947   0.002832   10.22
# muslim_j                            -0.088233   0.018832   -4.69
# mipex_j                             -0.010364   0.003905   -2.65
# relstate_j                          -0.017433   0.014196   -1.23
# claims_j                             0.122111   0.096072    1.27

. use "./../../assets/ESS7/ESS7e02_2.dta", cl. quietly do "./../../assets/ESS7/ESS7e02_2_formats_unicode.do"
. 
. * Micro data preparation

. ** Keep only those born in and with citizenship of country of interview, 

. ** Non-Muslims, and drop Israel

. keep if brncntr == 1 & ctzcntr == 1 & rlgdnm != 6 & cntry != "IL"
(6,789 observations deleted)

. 
. ** Add atheists to the religion variable

. replace rlgdnm = 0 if rlgblg == 2
(15,411 real changes made)

. 
. ** Educational attainment

. gen edu_attain = edulvlb
(93 missing values generated)

. recode edu_attain (0 113 129 = 1) (212 213 221 222 223 229 = 2) (311 312 313 
> 321 322 323 = 3) (412 413 421 422 423 = 4) (510 520 = 5) (610 620 = 6) (710 7
> 20 800 = 7) (else = .)
(edu_attain: 33396 changes made)

. label var edu_attain "Education"

. 
. ** Watching TV, but not news on TV

. gen TV_not_news = tvtot - tvpol
(1,506 missing values generated)

. label var TV_not_news "TV exposure"

. 
. ** Immigrant friends

. sum dfegcf

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
      dfegcf |     33,278    2.426979    .6874448          1          3

. replace dfegcf = `r(max)' - dfegcf
(33,396 real changes made, 118 to missing)

. lab var dfegcf "Friendships with immigrants"

. 
. ** Some more labels

. label define uempla 0 "Not unemployed" 1 "Unemployed", replace

. label variable hincfel "Economic deprivation"

. 
. ** Use the following variables for the analysis

. keep cntry almuslv agea gndr edu_attain uempla hincfel rlgblg rlgdnm dfegcf T
> V_not_news

. ** Listwise deletion of missing values

. foreach var in cntry almuslv agea gndr edu_attain uempla hincfel rlgblg rlgdn
> m dfegcf TV_not_news {
  2. drop if missing(`var')
  3. }
(0 observations deleted)
(1,217 observations deleted)
(54 observations deleted)
(20 observations deleted)
(160 observations deleted)
(0 observations deleted)
(242 observations deleted)
(114 observations deleted)
(65 observations deleted)
(75 observations deleted)
(1,385 observations deleted)

. 
. ** Save memory

. compress
  variable edu_attain was float now byte
  variable TV_not_news was float now byte
  (180,384 bytes saved)

. ** Make temporary file and save it

. tempfile micro_data

. save "`micro_data'"
file /var/folders/2_/0vc2cx65781cm0qm9tw72fg80000gn/T//St33330.000001 saved

. 
. * Macro data

. use "./../../assets/ESS7/Schlueter.dta", clear

. ** z-standardize macro predictor

. egen z_muslim_j = std(muslim_j)
(20 missing values generated)

. 
. * Conjoin Schlueter et al. macro data

. merge 1:m cntry using "`micro_data'", keep(match) nogen

    Result                           # of obs.
    -----------------------------------------
    not matched                             0
    matched                            30,064  
    -----------------------------------------

. 
. * Model 2

. mixed almuslv i.gndr agea edu_attain i.uempla hincfel dfegcf i.rlgdnm TV_not_
> news muslim_j mipex_j relstate_j claims_j || cntry: , reml

Performing EM optimization: 

Performing gradient-based optimization: 

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

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =     30,064
Group variable: cntry                           Number of groups  =         20

                                                Obs per group:
                                                              min =        953
                                                              avg =    1,503.2
                                                              max =      2,478

                                                Wald chi2(18)     =    4600.41
Log restricted-likelihood =  -36899.45          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
     almuslv |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        gndr |
     Female  |   .0017727   .0096184     0.18   0.854     -.017079    .0206244
        agea |   .0055072   .0002722    20.23   0.000     .0049737    .0060406
  edu_attain |  -.0903693   .0028245   -31.99   0.000    -.0959053   -.0848333
             |
      uempla |
 Unemployed  |   -.026586   .0250219    -1.06   0.288    -.0756281    .0224561
     hincfel |   .1015351   .0066859    15.19   0.000      .088431    .1146391
      dfegcf |  -.2344702   .0074458   -31.49   0.000    -.2490638   -.2198767
             |
      rlgdnm |
Roman Cat..  |   .0570756   .0131686     4.33   0.000     .0312655    .0828856
 Protestant  |    .032939   .0152821     2.16   0.031     .0029866    .0628914
Eastern O..  |   .0008523   .0543188     0.02   0.987    -.1056106    .1073151
Other Chr..  |  -.0774134    .050117    -1.54   0.122    -.1756409    .0208141
     Jewish  |  -.0224594   .1719952    -0.13   0.896    -.3595639     .314645
Eastern r..  |  -.0922595   .1033398    -0.89   0.372    -.2948018    .1102828
Other non..  |  -.2210224   .1008889    -2.19   0.028    -.4187609   -.0232839
             |
 TV_not_news |    .028947   .0028321    10.22   0.000     .0233962    .0344979
    muslim_j |  -.0882334   .0188324    -4.69   0.000    -.1251441   -.0513226
     mipex_j |  -.0103644   .0039052    -2.65   0.008    -.0180185   -.0027104
  relstate_j |  -.0174327   .0141957    -1.23   0.219    -.0452557    .0103903
    claims_j |   .1221106   .0960718     1.27   0.204    -.0661867    .3104079
       _cons |   3.528038   .2504508    14.09   0.000     3.037164    4.018913
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cntry: Identity              |
                  var(_cons) |   .0344445   .0127685      .0166564    .0712291
-----------------------------+------------------------------------------------
               var(Residual) |   .6774208   .0055284      .6666716    .6883433
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 935.93        Prob >= chibar2 = 0.0000

. 
  1. Add a random slope for “TV exposure”; interpret the results by drawing how you iagine the random slope to look like by hand.
ESS <- ESS %>%
  mutate(
    z_almuslv = (almuslv - mean(almuslv, na.rm = TRUE)) / sd(almuslv, na.rm = TRUE),
    z_TV_not_news = (TV_not_news - mean(TV_not_news, na.rm = TRUE)) / sd(TV_not_news, na.rm = TRUE))

# Model 4
Model_3 <- lmer(
  z_almuslv ~ gndr + agea + edu_attain + uempla + hincfel +
    dfegcf + rlgdnm + z_TV_not_news + muslim_j + mipex_j + 
    relstate_j + claims_j + (1 + z_TV_not_news | cntry), data  = ESS
)

summary(Model_3)
# Linear mixed model fit by REML ['lmerMod']
# Formula: z_almuslv ~ gndr + agea + edu_attain + uempla + hincfel + dfegcf +  
#     rlgdnm + z_TV_not_news + muslim_j + mipex_j + relstate_j +      claims_j + (1 + z_TV_not_news | cntry)
#    Data: ESS
# 
# REML criterion at convergence: 75542
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.561 -0.695  0.008  0.742  3.103 
# 
# Random effects:
#  Groups   Name          Variance Std.Dev. Corr 
#  cntry    (Intercept)   0.03771  0.1942        
#           z_TV_not_news 0.00148  0.0384   -0.22
#  Residual               0.71730  0.8469        
# Number of obs: 30064, groups:  cntry, 20
# 
# Fixed effects:
#                                      Estimate Std. Error t value
# (Intercept)                          0.925346   0.258058    3.59
# gndrMale                            -0.001809   0.009901   -0.18
# agea                                 0.005674   0.000281   20.23
# edu_attain                          -0.092874   0.002912  -31.90
# uemplaMarked                        -0.027001   0.025754   -1.05
# hincfel                              0.104613   0.006885   15.19
# dfegcf                              -0.240715   0.007669  -31.39
# rlgdnmEastern Orthodox               0.001690   0.055905    0.03
# rlgdnmEastern religions             -0.097854   0.106378   -0.92
# rlgdnmJewish                        -0.020189   0.177027   -0.11
# rlgdnmOther Christian denomination  -0.076974   0.051589   -1.49
# rlgdnmOther non-Christian religions -0.226086   0.103827   -2.18
# rlgdnmProtestant                     0.034904   0.015736    2.22
# rlgdnmRoman Catholic                 0.058991   0.013561    4.35
# z_TV_not_news                        0.052573   0.010113    5.20
# muslim_j                            -0.087053   0.019369   -4.49
# mipex_j                             -0.010417   0.004029   -2.59
# relstate_j                          -0.019605   0.014609   -1.34
# claims_j                             0.121544   0.098840    1.23


. quietly do "./../../assets/ESS7/9-exercise.. 
. * Z-standardization

. egen z_TV_not_news = std(TV_not_news)

. egen z_almuslv = std(almuslv)

. 
. ** Model 3

. mixed z_almuslv i.gndr agea edu_attain i.uempla hincfel dfegcf i.rlgdnm z_TV_
> not_news muslim_j mipex_j relstate_j claims_j || cntry: z_TV_not_news, reml

Performing EM optimization: 

Performing gradient-based optimization: 

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

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =     30,064
Group variable: cntry                           Number of groups  =         20

                                                Obs per group:
                                                              min =        953
                                                              avg =    1,503.2
                                                              max =      2,478

                                                Wald chi2(18)     =    4222.24
Log restricted-likelihood = -37771.034          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
   z_almuslv |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        gndr |
     Female  |   .0018031   .0099012     0.18   0.856     -.017603    .0212091
        agea |   .0056736   .0002805    20.22   0.000     .0051237    .0062234
  edu_attain |  -.0928767   .0029118   -31.90   0.000    -.0985837   -.0871696
             |
      uempla |
 Unemployed  |    -.02694   .0257548    -1.05   0.296    -.0774185    .0235385
     hincfel |    .104622   .0068853    15.20   0.000     .0911271    .1181169
      dfegcf |  -.2407303   .0076692   -31.39   0.000    -.2557615    -.225699
             |
      rlgdnm |
Roman Cat..  |   .0588494   .0135636     4.34   0.000     .0322653    .0854335
 Protestant  |   .0347839   .0157364     2.21   0.027     .0039411    .0656268
Eastern O..  |   .0019552   .0559046     0.03   0.972    -.1076158    .1115261
Other Chr..  |  -.0770574   .0515893    -1.49   0.135    -.1781706    .0240557
     Jewish  |  -.0196228   .1770271    -0.11   0.912    -.3665896    .3273441
Eastern r..  |  -.0978406   .1063807    -0.92   0.358     -.306343    .1106617
Other non..  |  -.2260544   .1038275    -2.18   0.029    -.4295526   -.0225562
             |
z_TV_not_n~s |   .0523385   .0101408     5.16   0.000     .0324629    .0722141
    muslim_j |  -.0918123   .0197136    -4.66   0.000    -.1304502   -.0531745
     mipex_j |  -.0106394    .004088    -2.60   0.009    -.0186517    -.002627
  relstate_j |  -.0166677   .0148585    -1.12   0.262    -.0457899    .0124546
    claims_j |   .1287731   .1005631     1.28   0.200     -.068327    .3258732
       _cons |   .9210309   .2619374     3.52   0.000      .407643    1.434419
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cntry: Independent           |
               var(z_TV_n~s) |   .0014856   .0006267      .0006498    .0033963
                  var(_cons) |    .037738    .013992      .0182466    .0780506
-----------------------------+------------------------------------------------
               var(Residual) |    .717295   .0058555      .7059099    .7288638
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 972.91                Prob > chi2 = 0.0000

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

. 

  1. Can the % Muslims in a country explain part of the random slope? Explain your results.
# Model 4
Model_4 <- lmer(
  z_almuslv ~ gndr + agea + edu_attain + uempla + hincfel +
    dfegcf + rlgdnm + z_TV_not_news*z_muslim_j + mipex_j + 
    relstate_j + claims_j + (1 + z_TV_not_news | cntry), data  = ESS
)

summary(Model_4)
# Linear mixed model fit by REML ['lmerMod']
# Formula: z_almuslv ~ gndr + agea + edu_attain + uempla + hincfel + dfegcf +  
#     rlgdnm + z_TV_not_news * z_muslim_j + mipex_j + relstate_j +      claims_j + (1 + z_TV_not_news | cntry)
#    Data: ESS
# 
# REML criterion at convergence: 75545
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.555 -0.695  0.007  0.743  3.106 
# 
# Random effects:
#  Groups   Name          Variance Std.Dev. Corr 
#  cntry    (Intercept)   0.03761  0.1939        
#           z_TV_not_news 0.00134  0.0366   -0.22
#  Residual               0.71730  0.8469        
# Number of obs: 30064, groups:  cntry, 20
# 
# Fixed effects:
#                                      Estimate Std. Error t value
# (Intercept)                          0.663256   0.268713    2.47
# gndrMale                            -0.001866   0.009901   -0.19
# agea                                 0.005677   0.000281   20.24
# edu_attain                          -0.092806   0.002912  -31.87
# uemplaMarked                        -0.026783   0.025754   -1.04
# hincfel                              0.104567   0.006885   15.19
# dfegcf                              -0.240643   0.007669  -31.38
# rlgdnmEastern Orthodox               0.001761   0.055905    0.03
# rlgdnmEastern religions             -0.098644   0.106378   -0.93
# rlgdnmJewish                        -0.018014   0.177032   -0.10
# rlgdnmOther Christian denomination  -0.076695   0.051589   -1.49
# rlgdnmOther non-Christian religions -0.226189   0.103827   -2.18
# rlgdnmProtestant                     0.035037   0.015736    2.23
# rlgdnmRoman Catholic                 0.059160   0.013561    4.36
# z_TV_not_news                        0.052603   0.009764    5.39
# z_muslim_j                          -0.231661   0.049479   -4.68
# mipex_j                             -0.010465   0.004028   -2.60
# relstate_j                          -0.019461   0.014604   -1.33
# claims_j                             0.122244   0.098800    1.24
# z_TV_not_news:z_muslim_j             0.015001   0.009830    1.53
# optimizer (nloptwrap) convergence code: 0 (OK)
# Model failed to converge with max|grad| = 0.0024852 (tol = 0.002, component 1)

. quietly do "./../../assets/ESS7/9-exercise.. 
. * Z-standardization

. egen z_TV_not_news = std(TV_not_news)

. egen z_almuslv = std(almuslv)

. 
. ** Model 4

. mixed z_almuslv i.gndr agea edu_attain i.uempla hincfel dfegcf i.rlgdnm c.z_T
> V_not_news##c.muslim_j mipex_j relstate_j claims_j || cntry: z_TV_not_news, r
> eml

Performing EM optimization: 

Performing gradient-based optimization: 

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

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =     30,064
Group variable: cntry                           Number of groups  =         20

                                                Obs per group:
                                                              min =        953
                                                              avg =    1,503.2
                                                              max =      2,478

                                                Wald chi2(19)     =    4233.82
Log restricted-likelihood = -37774.531          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
   z_almuslv |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        gndr |
     Female  |   .0018604   .0099012     0.19   0.851    -.0175456    .0212664
        agea |   .0056768   .0002805    20.24   0.000     .0051269    .0062266
  edu_attain |   -.092809   .0029121   -31.87   0.000    -.0985166   -.0871015
             |
      uempla |
 Unemployed  |  -.0267224   .0257551    -1.04   0.299    -.0772014    .0237566
     hincfel |   .1045781   .0068852    15.19   0.000     .0910833    .1180728
      dfegcf |  -.2406598   .0076692   -31.38   0.000    -.2556912   -.2256284
             |
      rlgdnm |
Roman Cat..  |   .0590203   .0135638     4.35   0.000     .0324358    .0856048
 Protestant  |   .0349144   .0157365     2.22   0.027     .0040715    .0657573
Eastern O..  |   .0020329   .0559045     0.04   0.971    -.1075379    .1116036
Other Chr..  |  -.0767766   .0515894    -1.49   0.137    -.1778899    .0243367
     Jewish  |  -.0174538   .1770322    -0.10   0.921    -.3644304    .3295229
Eastern r..  |  -.0986338   .1063808    -0.93   0.354    -.3071363    .1098687
Other non..  |  -.2261565   .1038273    -2.18   0.029    -.4296543   -.0226588
             |
z_TV_not_n~s |   .0344196   .0153439     2.24   0.025      .004346    .0644931
    muslim_j |  -.0917756   .0197202    -4.65   0.000    -.1304265   -.0531247
             |
          c. |
z_TV_not_n~s#|
  c.muslim_j |   .0059433   .0039112     1.52   0.129    -.0017225    .0136091
             |
     mipex_j |  -.0106688   .0040894    -2.61   0.009    -.0186839   -.0026537
  relstate_j |  -.0165834   .0148638    -1.12   0.265    -.0457158    .0125491
    claims_j |   .1291909   .1005977     1.28   0.199    -.0679769    .3263588
       _cons |   .9212703   .2620248     3.52   0.000     .4077111    1.434829
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cntry: Independent           |
               var(z_TV_n~s) |   .0013476   .0006016      .0005618    .0032326
                  var(_cons) |   .0377645   .0140022       .018259    .0781071
-----------------------------+------------------------------------------------
               var(Residual) |   .7172967   .0058555      .7059115    .7288656
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 972.43                Prob > chi2 = 0.0000

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

. 

References

Schlueter, Elmar, Anu Masso, and Eldad Davidov. 2019. “What Factors Explain Anti-Muslim Prejudice? An Assessment of the Effects of Muslim Population Size, Institutional Characteristics and Immigration-Related Media Claims.” Journal of Ethnic and Migration Studies 0(0):1–16. doi: 10.1080/1369183X.2018.1550160.