1. Revisit the three models of Table 2 in Green et al. (2020) The p-values of which slope estimates are accurate, and which ones are most likely too small?
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")library(lme4) # Multilevel modeling
library(texreg) # Regression tables
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#
# . 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//St33701.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//St33701.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)
#
# .# . 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
#
# .2. Approximate the df using the \(M-l-1\) rule for the DCE and CLIs.
The analyses are based on 20 countries; \(M
= 20\). Model 1 has c_unall_2011 and
c_gini_2011 as DCEs and a CLI between contact and MIPEX;
\(l = 5\). It follows that \(\text{df} \approx 20 - 5 - 1 = 14\) for
Model 1. Model 2 and Model 3 have three DCEs and no CLI. It follows that
\(\text{df} \approx 20 - 3 - 1 =
16\).
3. Calculate the p-value for everyday contact in Model 1.
# Let's print the model results to get the t-value
summary(Model_1)
# Linear mixed model fit by REML ['lmerMod']
# Formula: 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)
# Data: ESS
#
# REML criterion at convergence: 107999
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -3.663 -0.652 -0.050 0.602 3.913
#
# Random effects:
# Groups Name Variance Std.Dev. Corr
# cntry (Intercept) 0.32979 0.5743
# contact_ij 0.00252 0.0502 -0.56
# Residual 2.95983 1.7204
# Number of obs: 27472, groups: cntry, 20
#
# Fixed effects:
# Estimate Std. Error t value
# (Intercept) 6.76e+00 1.06e+00 6.37
# contact_ij -7.78e-02 1.27e-02 -6.12
# mipex_j -1.04e-02 1.13e-02 -0.92
# agea_ij 1.26e-03 6.58e-04 1.91
# gndr_ijFemale -4.36e-02 2.12e-02 -2.06
# eduyrs_ij -1.01e-01 3.07e-03 -32.92
# immiback_ijYes -4.25e-01 3.29e-02 -12.92
# subj_income_ij -2.70e-01 1.44e-02 -18.74
# relig_ij -4.07e-02 3.81e-03 -10.68
# conserve_ij 3.23e-01 1.41e-02 22.93
# div_nbh_ij -2.50e-02 1.67e-02 -1.50
# c_unall_2011 1.47e-05 1.04e-04 0.14
# c_gini_2011 -1.91e-02 3.86e-02 -0.49
# contact_ij:mipex_j -2.40e-03 1.06e-03 -2.27
# fit warnings:
# Some predictor variables are on very different scales: consider rescaling
# optimizer (nloptwrap) convergence code: 0 (OK)
# Model failed to converge with max|grad| = 0.00269264 (tol = 0.002, component 1)
# Plugging in the t-value and the df we can calculate the p-value.
2 * (1 - pt(q = 6.12, df = 13))
# [1] 3.7e-05#
# . quietly do "./../../assets/ESS7/10-exercise..
# . *// Model results
#
# . mixed symb_threat_ij c.contact_ij##c.mipex_j agea_ij gndr_ij eduyrs_ij immiba
# > ck_ij subj_income_ij relig_ij conserve_ij div_nbh_ij c_unall_2011 c_gini_2011
# > || cntry: contact_ij , reml covariance(unstructured)
#
# Performing EM optimization:
#
# Performing gradient-based optimization:
#
# Iteration 0: log restricted-likelihood = -64386.546
# Iteration 1: log restricted-likelihood = -64386.546
#
# Computing standard errors:
#
# Mixed-effects REML regression Number of obs = 32,461
# Group variable: cntry Number of groups = 20
#
# Obs per group:
# min = 1,095
# avg = 1,623.0
# max = 2,732
#
# Wald chi2(13) = 3398.68
# Log restricted-likelihood = -64386.546 Prob > chi2 = 0.0000
#
# ------------------------------------------------------------------------------
# symb_threa~j | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------
# contact_ij | -.0839535 .0119149 -7.05 0.000 -.1073063 -.0606008
# mipex_j | -.007093 .0106961 -0.66 0.507 -.028057 .013871
# |
# c. |
# contact_ij#|
# c.mipex_j | -.0027575 .000988 -2.79 0.005 -.0046939 -.000821
# |
# agea_ij | .0023729 .0006087 3.90 0.000 .0011799 .0035659
# gndr_ij | -.0510287 .0198686 -2.57 0.010 -.0899704 -.012087
# eduyrs_ij | -.0927992 .0027498 -33.75 0.000 -.0981886 -.0874097
# immiback_ij | -.5710364 .0419352 -13.62 0.000 -.6532279 -.4888449
# subj_incom~j | -.2841113 .0133002 -21.36 0.000 -.3101791 -.2580435
# relig_ij | -.0428884 .0035699 -12.01 0.000 -.0498853 -.0358915
# conserve_ij | .3226175 .0131374 24.56 0.000 .2968686 .3483664
# div_nbh_ij | -.0164396 .0156984 -1.05 0.295 -.0472079 .0143286
# c_unall_2011 | .0000188 .0001034 0.18 0.855 -.0001837 .0002214
# c_gini_2011 | -.0190327 .0381957 -0.50 0.618 -.0938949 .0558294
# _cons | -1.38316 1.049523 -1.32 0.188 -3.440187 .6738673
# ------------------------------------------------------------------------------
#
# ------------------------------------------------------------------------------
# Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
# -----------------------------+------------------------------------------------
# cntry: Unstructured |
# var(contac~j) | .0022144 .0009016 .000997 .0049184
# var(_cons) | .2973545 .1088003 .1451526 .6091501
# cov(contac~j,_cons) | -.0126495 .0079574 -.0282458 .0029468
# -----------------------------+------------------------------------------------
# var(Residual) | 3.071164 .0241247 3.024242 3.118813
# ------------------------------------------------------------------------------
# LR test vs. linear model: chi2(3) = 1654.63 Prob > chi2 = 0.0000
#
# Note: LR test is conservative and provided only for reference.
#
# .
# . # Plugging in the t-value and the df we can calculate the p-value.
# Unknown #command
#
# . dis 2 * (1 - t(13, 7.05))
# 8.671e-06
#
# .4. How does that compare to the result based on the numerical approximation?
library(iimm) # Library for numerical approximation.
# Add numerically approximated df and accurate p-values to model object.
Model_1_t <- lmer_t(Model_1, method = "Satterthwaite")
# Print the results, which are pretty similar but habe 2 df more.
summary(Model_1_t)
# Linear mixed model fit by REML ['lmerMod']
# t-tests use the Satterthwaite method.
# Formula: 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)
# Data: ESS
#
# REML criterion at convergence: 107999
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -3.663 -0.652 -0.050 0.602 3.913
#
# Random effects:
# Groups Name Variance Std.Dev. Corr
# cntry (Intercept) 0.32979 0.5743
# contact_ij 0.00252 0.0502 -0.56
# Residual 2.95983 1.7204
# Number of obs: 27472, groups: cntry, 20
#
# Coefficients:
# Estimate Std.Err t value Df Lower Upper Pr(>|t|)
# (Intercept) 6.762 1.061 6.372 16 4.515 9.008 8.51e-06 ***
# contact_ij -0.078 0.013 -6.117 20 -0.104 -0.051 6.05e-06 ***
# mipex_j -0.010 0.011 -0.924 17 -0.034 0.013 0.3683
# agea_ij 0.001 0.001 1.910 27442 0.000 0.003 0.0561 .
# gndr_ijFemale -0.044 0.021 -2.063 27439 -0.085 -0.002 0.0391 *
# eduyrs_ij -0.101 0.003 -32.922 27453 -0.107 -0.095 < 2e-16 ***
# immiback_ijYes -0.425 0.033 -12.916 27437 -0.490 -0.361 < 2e-16 ***
# subj_income_ij -0.270 0.014 -18.743 27446 -0.298 -0.241 < 2e-16 ***
# relig_ij -0.041 0.004 -10.676 27456 -0.048 -0.033 < 2e-16 ***
# conserve_ij 0.323 0.014 22.935 27453 0.295 0.350 < 2e-16 ***
# div_nbh_ij -0.025 0.017 -1.499 27451 -0.058 0.008 0.1339
# c_unall_2011 0.000 0.000 0.141 16 0.000 0.000 0.8896
# c_gini_2011 -0.019 0.039 -0.494 16 -0.101 0.063 0.6278
# contact_ij:mipex_j -0.002 0.001 -2.268 19 -0.005 0.000 0.0351 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# fit warnings:
# Some predictor variables are on very different scales: consider rescaling
# optimizer (nloptwrap) convergence code: 0 (OK)
# Model failed to converge with max|grad| = 0.00269264 (tol = 0.002, component 1)#
# . quietly do "./../../assets/ESS7/10-exercise..
# . *// Model with numerical df approximation
#
# . mixed symb_threat_ij c.contact_ij##c.mipex_j agea_ij gndr_ij eduyrs_ij immiba
# > ck_ij subj_income_ij relig_ij conserve_ij div_nbh_ij c_unall_2011 c_gini_2011
# > || cntry: contact_ij , reml covariance(unstructured) dfmethod(satterthwaite)
# > dftable(pvalue)
#
# Performing EM optimization:
#
# Performing gradient-based optimization:
#
# Iteration 0: log restricted-likelihood = -64386.546
# Iteration 1: log restricted-likelihood = -64386.546
#
# Computing standard errors:
#
# Computing degrees of freedom:
#
# Mixed-effects REML regression Number of obs = 32,461
# Group variable: cntry Number of groups = 20
#
# Obs per group:
# min = 1,095
# avg = 1,623.0
# max = 2,732
# DF method: Satterthwaite DF: min = 15.97
# avg = 18,540.96
# max = 32,443.69
#
# F(13, 41.49) = 261.44
# Log restricted-likelihood = -64386.546 Prob > F = 0.0000
#
# -----------------------------------------------------------------------------
# symb_threat_ij | Coef. Std. Err. DF t P>|t|
# -----------------------+-----------------------------------------------------
# contact_ij | -.0839535 .0119149 19.6 -7.05 0.000
# mipex_j | -.007093 .0106961 16.6 -0.66 0.516
# |
# c.contact_ij#c.mipex_j | -.0027575 .000988 18.8 -2.79 0.012
# |
# agea_ij | .0023729 .0006087 32421.9 3.90 0.000
# gndr_ij | -.0510287 .0198686 32429.4 -2.57 0.010
# eduyrs_ij | -.0927992 .0027498 32442.5 -33.75 0.000
# immiback_ij | -.5710364 .0419352 32421.4 -13.62 0.000
# subj_income_ij | -.2841113 .0133002 32427.2 -21.36 0.000
# relig_ij | -.0428884 .0035699 32443.7 -12.01 0.000
# conserve_ij | .3226175 .0131374 32443.3 24.56 0.000
# div_nbh_ij | -.0164396 .0156984 32440.6 -1.05 0.295
# c_unall_2011 | .0000188 .0001034 16.0 0.18 0.858
# c_gini_2011 | -.0190327 .0381957 16.1 -0.50 0.625
# _cons | -1.38316 1.049523 16.4 -1.32 0.206
# -----------------------------------------------------------------------------
#
# ------------------------------------------------------------------------------
# Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
# -----------------------------+------------------------------------------------
# cntry: Unstructured |
# var(contac~j) | .0022144 .0009016 .000997 .0049184
# var(_cons) | .2973545 .1088003 .1451526 .6091501
# cov(contac~j,_cons) | -.0126495 .0079574 -.0282458 .0029468
# -----------------------------+------------------------------------------------
# var(Residual) | 3.071164 .0241247 3.024242 3.118813
# ------------------------------------------------------------------------------
# LR test vs. linear model: chi2(3) = 1654.63 Prob > chi2 = 0.0000
#
# Note: LR test is conservative and provided only for reference.
#
# .5. Add the random slope for contact to Model 2 Is the additional random slope a significant improvement of the model?
# Estimate model with additional random slope.
Model_2b <- 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 + contact_ij | cntry))
# Compare the two nested models.
# BIC suggests the added random slope is no improvement,
# but the likelihood ratio tests is significant. Since the more complex model
# doesn't even converge, better leave the random slope out.
anova(Model_2, Model_2b)
# Data: ESS
# Models:
# Model_2: 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_2b: 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 + contact_ij | cntry)
# npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
# Model_2 15 103251 103374 -51610 103221
# Model_2b 17 103247 103387 -51607 103213 7.93 2 0.019 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#
# . quietly do "./../../assets/ESS7/10-exercise..
# . *// Estimate original model
#
# . quietly mixed real_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: , reml
#
# . estimates store Model_2
#
# . *// Estimate model with additional random slope.
#
# . quietly mixed real_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_2b
#
# .
# . *// Compare the two nested models.
#
# . *// BIC suggests the added random slope is no improvement,
#
# . *// and with a p-value of 0.23 the likelihood ratio test is insignificant, to
# > o.
#
# . lrtest Model_2 Model_2b, stats
#
# Likelihood-ratio test LR chi2(2) = 2.01
# (Assumption: Model_2 nested in Model_2b) Prob > chi2 = 0.3664
#
# Note: The reported degrees of freedom assumes the null hypothesis is not on
# the boundary of the parameter space. If this is not true, then the
# reported test is conservative.
# Note: LR tests based on REML are valid only when the fixed-effects
# specification is identical for both models.
#
# Akaike's information criterion and Bayesian information criterion
#
# -----------------------------------------------------------------------------
# Model | Obs ll(null) ll(model) df AIC BIC
# -------------+---------------------------------------------------------------
# Model_2 | 32,461 . -61794.86 16 123621.7 123755.9
# Model_2b | 32,461 . -61793.86 18 123623.7 123774.7
# -----------------------------------------------------------------------------
# Note: N=Obs used in calculating BIC; see [R] BIC note.
#
# .