1. Add the context variables MIPEX (from Schlueter.dta), the Gini and unemployment rate from year 2011 (from the ESS country data) to your dataset.
library(tidyverse) # Data manipulation
library(haven) # Read Stata dta files
library(essurvey) # API to ESS
# Import the ESS round 7 data via the API
ESS <- import_rounds(rounds = 7, ess_email = YOUR_EMAIL) %>%
recode_missings() %>%
mutate(
# Threat & contact
symb_threat_ij = 10 - ((zap_labels(imueclt) + zap_labels(imwbcnt) + zap_labels(rlgueim)) / 3),
real_threat_ij = 10 - ((zap_labels(imtcjob) + zap_labels(imbleco) + zap_labels(imwbcrm)) / 3),
contact_ij = zap_labels(dfegcon), # Everyday contact
## Controls
subj_income_ij = 4 - zap_labels(hincfel),
eduyrs_ij = case_when(
eduyrs > mean(eduyrs, na.rm = TRUE) + 3*sd(eduyrs, na.rm = TRUE) ~
mean(eduyrs, na.rm = TRUE) + 3*sd(eduyrs, na.rm = TRUE),
TRUE ~ zap_labels(eduyrs)),
agea_ij = zap_labels(agea),
relig_ij = zap_labels(rlgdgr),
div_nbh_ij = zap_labels(acetalv),
immiback_ij = case_when(
as_factor(facntr) == "No" | as_factor(mocntr) == "No" ~ "Yes",
as_factor(facntr) == "Yes" | as_factor(mocntr) == "Yes" ~ "No",
TRUE ~ as.character(NA)) %>% factor(),
conserve_ij = 6 - ((zap_labels(impsafe) + zap_labels(ipstrgv) + zap_labels(ipfrule) + zap_labels(ipmodst) + zap_labels(ipbhprp) + zap_labels(imptrad)) / 6),
# Categorical variables
cntry = as_factor(cntry) %>% fct_drop(),
ctzcntr = as_factor(ctzcntr) %>% fct_drop(),
gndr_ij = as_factor(gndr) %>% fct_drop()
) %>%
filter(cntry != "IL" & agea >= 18 & ctzcntr == "Yes") %>%
select(symb_threat_ij, real_threat_ij, contact_ij, eduyrs_ij, agea_ij, gndr_ij,
subj_income_ij, relig_ij, div_nbh_ij, immiback_ij, conserve_ij, cntry) %>%
drop_na() %>%
inner_join(.,
read_dta("./../../assets/ESS7/Schlueter.dta") %>%
mutate(mipex_j = mipex_j - mean(mipex_j, na.rm = TRUE)),
by = "cntry") %>%
inner_join(.,
read_dta("./../../assets/ESS7/ESSMD-2014-cntry_F1.dta"),
by = "cntry")#
# . use "./../../assets/ESS7/ESS7e02_2.dta", cl. quietly do "./../../assets/ESS7/ESS7e02_2_formats_unicode.do"
# . keep if cntry != "IL" & agea >= 18 & ctzcntr == 1
# (5,757 observations deleted)
#
# .
# . *// Threat and contact
#
# . alpha imueclt imwbcnt rlgueim, gen(symb_threat_ij) reverse(imueclt imwbcnt rl
# > gueim)
#
# Test scale = mean(unstandardized items)
# Reversed items: imueclt imwbcnt rlgueim
#
# Average interitem covariance: 2.983469
# Number of items in the scale: 3
# Scale reliability coefficient: 0.7940
#
# . alpha imtcjob imbleco imwbcrm, gen(real_threat_ij) reverse(imtcjob imbleco im
# > wbcrm)
#
# Test scale = mean(unstandardized items)
# Reversed items: imtcjob imbleco imwbcrm
#
# Average interitem covariance: 2.132521
# Number of items in the scale: 3
# Scale reliability coefficient: 0.7153
#
# . rename dfegcon contact_ij
#
# .
# . *// Controls
#
# . gen subj_income_ij = 4 - hincfel
# (241 missing values generated)
#
# . rename eduyrs eduyrs_ij
#
# . sum eduyrs_ij
#
# Variable | Obs Mean Std. Dev. Min Max
# -------------+---------------------------------------------------------
# eduyrs_ij | 34,130 12.98186 3.955893 0 50
#
# . replace eduyrs_ij = `r(mean)' + 3*`r(sd)' if eduyrs_ij > `r(mean)' + 3*`r(sd)
# > '
# variable eduyrs_ij was byte now float
# (482 real changes made)
#
# . rename agea agea_ij
#
# . rename rlgdgr relig_ij
#
# . rename acetalv div_nbh_ij
#
# . gen immiback_ij = 1 if facntr == 2 | mocntr == 2
# (30,226 missing values generated)
#
# . replace immiback_ij = 0 if facntr == 1 | mocntr == 1
# (32,265 real changes made)
#
# . rename gndr gndr_ij
#
# . alpha impsafe ipstrgv ipfrule ipmodst ipbhprp imptrad, gen(conserve_ij) rever
# > se(impsafe ipstrgv ipfrule ipmodst ipbhprp imptrad)
#
# Test scale = mean(unstandardized items)
# Reversed items: impsafe ipstrgv ipfrule ipmodst ipbhprp imptrad
#
# Average interitem covariance: .4675392
# Number of items in the scale: 6
# Scale reliability coefficient: 0.7072
#
# .
# . *// Keep variables and only complete obserbvations
#
# . keep cntry symb_threat_ij real_threat_ij contact_ij subj_income_ij eduyrs_ij
# > agea_ij relig_ij div_nbh_ij immiback_ij gndr_ij conserve_ij
#
# . foreach var in cntry symb_threat_ij real_threat_ij contact_ij subj_income_ij
# > eduyrs_ij agea_ij relig_ij div_nbh_ij immiback_ij gndr_ij conserve_ij {
# 2. drop if missing(`var')
# 3. }
# (0 observations deleted)
# (250 observations deleted)
# (175 observations deleted)
# (286 observations deleted)
# (215 observations deleted)
# (0 observations deleted)
# (50 observations deleted)
# (172 observations deleted)
# (249 observations deleted)
# (40 observations deleted)
# (17 observations deleted)
# (513 observations deleted)
#
# .
# . *// Save temporary file
#
# . tempfile micro_data
#
# . save `micro_data'
# file /var/folders/2_/0vc2cx65781cm0qm9tw72fg80000gn/T//St33560.000001 saved
#
# .
# . *// Conjoin the Schlueter data
#
# . use "./../../assets/ESS7/Schlueter.dta", clear
#
# . *// !! grand-mean center mipex !!
#
# . sum mipex_j
#
# Variable | Obs Mean Std. Dev. Min Max
# -------------+---------------------------------------------------------
# mipex_j | 20 56.9 12.25561 38 80
#
# . replace mipex_j = mipex_j - `r(mean)'
# variable mipex_j was byte now float
# (20 real changes made)
#
# . merge 1:m cntry using `micro_data' , keep(match) nogen
#
# Result # of obs.
# -----------------------------------------
# not matched 0
# matched 32,461
# -----------------------------------------
#
# .
# . *// Save temporary file
#
# . tempfile micro_data
#
# . save `micro_data'
# file /var/folders/2_/0vc2cx65781cm0qm9tw72fg80000gn/T//St33560.000002 saved
#
# .
# . *// The ESS country data
#
# . use cntry c_gini_2011 c_unall_2011 using "./../../assets/ESS7/ESSMD-2014-cntr
# > y_F1", clear
#
# . merge 1:m cntry using `micro_data' , keep(match) nogen
#
# Result # of obs.
# -----------------------------------------
# not matched 0
# matched 32,461
# -----------------------------------------
#
# . compress
# variable subj_income_ij was float now byte
# variable immiback_ij was float now byte
# (194,766 bytes saved)
#
# .2. Replicate all three models of Table 2 in Green et al. (2020), although we lack the change in the immigration rate.
library(lme4)
library(texreg)
Model_1 <- lmer(data = ESS, symb_threat_ij ~ contact_ij*mipex_j + agea_ij + gndr_ij + eduyrs_ij + immiback_ij + subj_income_ij + relig_ij + conserve_ij + div_nbh_ij + c_unall_2011 + c_gini_2011 + (1 + contact_ij | cntry))
Model_2 <- lmer(data = ESS, real_threat_ij ~ contact_ij + mipex_j + agea_ij + gndr_ij + eduyrs_ij + immiback_ij + subj_income_ij + relig_ij + conserve_ij + div_nbh_ij + c_unall_2011 + c_gini_2011 + (1 | cntry))
Model_3 <- lmer(data = ESS, contact_ij ~ mipex_j + agea_ij + gndr_ij + eduyrs_ij + immiback_ij + subj_income_ij + relig_ij + conserve_ij + div_nbh_ij + c_unall_2011 + c_gini_2011 + (1 | cntry))
# Regression table of correct and mis-specified model
screenreg(list(Model_1, Model_2, Model_3), single.row = FALSE, digits = 3)
#
# =================================================================================
# Model 1 Model 2 Model 3
# ---------------------------------------------------------------------------------
# (Intercept) 6.762 *** 8.164 *** 4.769 ***
# (1.061) (0.825) (1.201)
# contact_ij -0.078 *** -0.046 ***
# (0.013) (0.005)
# mipex_j -0.010 -0.009 0.032 **
# (0.011) (0.007) (0.011)
# agea_ij 0.001 -0.003 *** -0.024 ***
# (0.001) (0.001) (0.001)
# gndr_ijFemale -0.044 * 0.064 *** -0.002
# (0.021) (0.019) (0.022)
# eduyrs_ij -0.101 *** -0.062 *** 0.065 ***
# (0.003) (0.003) (0.003)
# immiback_ijYes -0.425 *** -0.389 *** 0.243 ***
# (0.033) (0.030) (0.034)
# subj_income_ij -0.270 *** -0.267 *** 0.091 ***
# (0.014) (0.013) (0.015)
# relig_ij -0.041 *** -0.035 *** -0.005
# (0.004) (0.003) (0.004)
# conserve_ij 0.323 *** 0.217 *** -0.110 ***
# (0.014) (0.013) (0.015)
# div_nbh_ij -0.025 -0.010 0.842 ***
# (0.017) (0.015) (0.016)
# c_unall_2011 0.000 0.000 0.000
# (0.000) (0.000) (0.000)
# c_gini_2011 -0.019 -0.050 -0.044
# (0.039) (0.030) (0.044)
# contact_ij:mipex_j -0.002 *
# (0.001)
# ---------------------------------------------------------------------------------
# AIC 108034.862 103358.093 109977.012
# BIC 108182.839 103481.407 110092.105
# Log Likelihood -53999.431 -51664.047 -54974.506
# Num. obs. 27472 27472 27472
# Num. groups: cntry 20 20 20
# Var: cntry (Intercept) 0.330 0.138 0.297
# Var: cntry contact_ij 0.003
# Cov: cntry (Intercept) contact_ij -0.016
# Var: Residual 2.960 2.501 3.183
# =================================================================================
# *** p < 0.001; ** p < 0.01; * p < 0.05# . quietly do "./../../assets/ESS7/10-exercise.do"
# .
# . quietly mixed symb_threat_ij c.contact_ij##c.mipex_j agea_ij gndr_ij eduyrs_i
# > j immiback_ij subj_income_ij relig_ij conserve_ij div_nbh_ij c_unall_2011 c_g
# > ini_2011 || cntry: contact_ij , reml covariance(unstructured)
#
# . estimates store Model_1
#
# .
# . quietly mixed real_threat_ij contact_ij mipex_j agea_ij gndr_ij eduyrs_ij imm
# > iback_ij subj_income_ij relig_ij conserve_ij div_nbh_ij c_unall_2011 c_gini_2
# > 011 || cntry: , reml
#
# . estimates store Model_2
#
# .
# . quietly mixed contact_ij c.mipex_j agea_ij gndr_ij eduyrs_ij immiback_ij subj
# > _income_ij relig_ij conserve_ij div_nbh_ij c_unall_2011 c_gini_2011 || cntry:
# > , reml
#
# . estimates store Model_3
#
# .
# . * Make regression table (change names of vars in fields (e.g. "contact_ij"))
#
# . esttab Model_1 Model_2 Model_3, b(3) se(3) nobaselevels mtitles("Model 1" "Mo
# > del 2" "Model 3") nonumbers transform(ln*: exp(@)^2 exp(@)^2 at*: tanh(@) (1-
# > tanh(@)^2)) eqlabels("" "var(contact_ij)" "var(_cons)" "corr(contact_ij, _con
# > s)" "var(Residual)", none) varlabels(, elist(weight:_cons "{break}{hline @wid
# > th}")) compress
#
# -------------------------------------------------
# Model 1 Model 2 Model 3
# -------------------------------------------------
# contact_ij -0.084*** -0.051***
# (0.012) (0.005)
# mipex_j -0.007 -0.008 0.032**
# (0.011) (0.007) (0.011)
# c.contac~j -0.003**
# (0.001)
# agea_ij 0.002*** -0.001* -0.025***
# (0.001) (0.001) (0.001)
# gndr_ij -0.051* 0.051** -0.016
# (0.020) (0.018) (0.020)
# eduyrs_ij -0.093*** -0.059*** 0.061***
# (0.003) (0.003) (0.003)
# immiback~j -0.571*** -0.550*** 0.248***
# (0.042) (0.039) (0.043)
# subj_inc~j -0.284*** -0.280*** 0.099***
# (0.013) (0.012) (0.014)
# relig_ij -0.043*** -0.035*** -0.006
# (0.004) (0.003) (0.004)
# conserve~j 0.323*** 0.221*** -0.120***
# (0.013) (0.012) (0.013)
# div_nbh_ij -0.016 0.004 0.881***
# (0.016) (0.014) (0.015)
# c_una~2011 0.000 0.000 0.000
# (0.000) (0.000) (0.000)
# c_gin~2011 -0.019 -0.046 -0.044
# (0.038) (0.030) (0.044)
# _cons -1.383 -0.750 4.121***
# (1.050) (0.813) (1.218)
# var(cont~) 0.002*** 0.135*** 0.307***
# (0.000) (0.024) (0.055)
# var(_cons) 0.297***
# (0.054)
# corr(con~) -0.493
# (0.213)
# var(Resi~) 3.071*** 2.621*** 3.214***
# (0.012) (0.010) (0.013)
# -------------------------------------------------
# N 32461 32461 32461
# -------------------------------------------------
# Standard errors in parentheses
# * p<0.05, ** p<0.01, *** p<0.001
#
# .