Model_explanation
2025-03-12
Load libraries
setwd("/Users/marvin-julianstruckmeyer/MCSS_Survey_Research_Methodology_II_Project") # delete later
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lme4) # For multilevel modeling## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(lmerTest) # For p-values in multilevel models##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## steplibrary(ggplot2) # For data visualization
library(dplyr) # For data manipulation
library(sjstats) # For calculating ICC
library(sjPlot) # For plotting model results## Learn more about sjPlot with 'browseVignettes("sjPlot")'.library(modelsummary)
library(broom.mixed) # For extracting mixed model statistics
library(stargazer)##
## Please cite as:
##
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(rtf)##
## Attaching package: 'rtf'
##
## The following object is masked from 'package:purrr':
##
## done
##
## The following object is masked from 'package:tibble':
##
## viewlibrary(haven)
data <- read_dta("data/raw/ZA7575.dta")
df_reduced <- readRDS("df_reduced.rds")
# restore the old NAs and do CCA
df_reduced$qc19 <- data$qc19
df_reduced$qc19 <- ifelse(df_reduced$qc19 == 3, NA, df_reduced$qc19)
df_reduced <- df_reduced %>%
drop_na(qc19)
# quick data check
table(df_reduced$country_name)##
## Austria Belgium Bulgaria Croatia Cyprus
## 957 1014 738 912 399
## Czech Republic Denmark Estonia Finland France
## 919 930 803 883 915
## Germany Greece Hungary Ireland Italy
## 1333 923 903 882 889
## Latvia Lithuania Luxembourg Malta Netherlands
## 847 851 468 448 980
## Poland Portugal Romania Slovakia Slovenia
## 818 911 905 886 920
## Spain Sweden United Kingdom
## 921 902 901Visualize qc19 by country
# Calculate the percentage of "Yes" (1) responses for qc19 by country
country_means <- aggregate(qc19 == 1 ~ country_name, data = df_reduced, FUN = mean)
# Multiply by 100 to get percentage and then sort
country_means$percentage <- country_means$`qc19 == 1` * 100
country_means <- country_means[order(country_means$percentage), ]
# Create the plot
ggplot(country_means, aes(x = reorder(country_name, percentage), y = percentage)) +
geom_bar(stat = "identity", fill = "steelblue") +
theme_minimal() +
labs(title = "Support for Transgender Rights by Country",
subtitle = "Percentage answering 'Yes' to allowing transgender people to change civil documents",
x = "Country",
y = "Support (%)") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Prepare for GLM
# scale some of the variables (both categorical and continuous) with larger scales
df_reduced <- df_reduced %>%
mutate(
age_z = scale(age),
religion_z = scale(as.numeric(religion)),
political_ideology_z = scale(political_ideology),
respondent_cooperation_z = scale(as.numeric(respondent_cooperation)),
rainbow_score_z = scale(rainbow_score_2019),
v2x_egaldem_z = scale(v2x_egaldem))
# Recode qc19 from 1 and 2 to 0 and 1 for the binary model
df_reduced$qc19 <- ifelse(df_reduced$qc19 == 1, 1, 0)
# Optimize the settings to avoid convergence problems in the glmer() functions
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 100000))
# Recode some of the variables such that higher values translate to a more positive
# attitude towards LGBT rights in general -> simplieies coefficient interpretation
df_reduced <- df_reduced %>%
mutate(trans_friends = ifelse(trans_friends == 1, 2, 1), # trans_friends is sd1_7 (see Reduced dataset.R)
lgb_friends = ifelse(lgb_friends == 1, 2, 1)) # lgb_friends is sd1_4 (see Reduced dataset.R)1. Null model (intercept only)
There is significant variation in qc19 scores between countries, justifying the need for a multilevel approach. About 12% of the differences in qc19 scores can be attributed to country-level factors, while 88% is due to individual differences. 12% of variation on the country-level is considered a lot, as the threshold for going ahead with multilevel models is usually 5%.
# STEP 1: Null Model (Intercept-only model)
null_model <- glmer(qc19 ~ 1 + (1|country_name), family = "binomial", data = df_reduced)
summary(null_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: qc19 ~ 1 + (1 | country_name)
## Data: df_reduced
##
## AIC BIC logLik deviance df.resid
## 28212.6 28228.8 -14104.3 28208.6 24156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0233 -0.9008 0.5066 0.6531 2.2755
##
## Random effects:
## Groups Name Variance Std.Dev.
## country_name (Intercept) 0.9649 0.9823
## Number of obs: 24158, groups: country_name, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4495 0.1856 2.422 0.0154 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Calculate ICC
performance::icc(null_model)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.227
## Unadjusted ICC: 0.2272. Individual-level predictors (level 1)
Individual characteristics and interview quality (respondent cooperation) significantly predict qc19 scores. The strongest predictors are having LGB friends and respondent cooperation. Religion has a negative impact, suggesting that more religious individuals tend to have lower qc19 scores. Even after controlling for these individual factors, significant country-level differences remain.
# STEP 2: Add Individual-Level Predictors (Level 1)
# starting with predictors that showed correlation with qc19
model1 <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1|country_name),
family = "binomial",
control = control,
data = df_reduced)
summary(model1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
## trans_friends + respondent_cooperation_z + (1 | country_name)
## Data: df_reduced
## Control: control
##
## AIC BIC logLik deviance df.resid
## 26506.1 26579.0 -13244.1 26488.1 24149
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.0609 -0.7753 0.3671 0.6619 4.7814
##
## Random effects:
## Groups Name Variance Std.Dev.
## country_name (Intercept) 0.7199 0.8485
## Number of obs: 24158, groups: country_name, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.82055 0.18210 -9.997 < 2e-16 ***
## age_z -0.07919 0.01614 -4.907 9.26e-07 ***
## gender 0.38580 0.03065 12.587 < 2e-16 ***
## religion_z 0.10622 0.01734 6.127 8.95e-10 ***
## political_ideology_z -0.02151 0.01541 -1.396 0.163
## lgb_friends 0.88726 0.03648 24.323 < 2e-16 ***
## trans_friends 0.42426 0.06044 7.020 2.22e-12 ***
## respondent_cooperation_z -0.27334 0.01612 -16.960 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age_z gender rlgn_z pltc__ lgb_fr trns_f
## age_z -0.068
## gender -0.240 -0.027
## religion_z 0.014 0.124 0.095
## pltcl_dlgy_ -0.002 0.078 -0.062 -0.020
## lgb_friends -0.182 0.212 -0.044 -0.079 0.039
## trans_frnds -0.295 0.037 -0.008 -0.030 0.008 -0.222
## rspndnt_cp_ -0.021 -0.084 -0.020 -0.015 -0.064 0.067 0.017# compare with the null_model
anova(null_model, model1)## Data: df_reduced
## Models:
## null_model: qc19 ~ 1 + (1 | country_name)
## model1: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1 | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 2 28213 28229 -14104 28209
## model1 9 26506 26579 -13244 26488 1720.5 7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelsummary(
list("Model 1" = null_model,
"Model 2" = model1),
stars = TRUE,
gof_map = c("nobs", "aic", "bic", "logLik", "r.squared"))| Model 1 | Model 2 | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 0.449* | -1.821*** |
| (0.186) | (0.182) | |
| age_z | -0.079*** | |
| (0.016) | ||
| gender | 0.386*** | |
| (0.031) | ||
| religion_z | 0.106*** | |
| (0.017) | ||
| political_ideology_z | -0.022 | |
| (0.015) | ||
| lgb_friends | 0.887*** | |
| (0.036) | ||
| trans_friends | 0.424*** | |
| (0.060) | ||
| respondent_cooperation_z | -0.273*** | |
| (0.016) | ||
| SD (Intercept country_name) | 0.982 | 0.848 |
| Num.Obs. | 24158 | 24158 |
| AIC | 28212.6 | 26506.1 |
| BIC | 28228.8 | 26579.0 |
3. Country-level predictors (level 2)
Country-level factors, particularly those related to LGBTQ+ rights (rainbow score) and egalitarian democracy, strongly explain why qc19 scores differ between countries. The negative relationships suggest that more progressive and egalitarian countries tend to have lower qc19 scores. These two country-level variables explain a remarkably large portion (74.5%) of the between-country differences.
The change from gender_equality_index to v2x_egaldem has provided a stronger country-level effect, suggesting that broader measures of egalitarianism in a society may be more relevant to explaining variation in qc19 than gender equality specifically.
# STEP 3: Add Country-Level Predictors (Level 2)
# We test three different models based on 'themes'. As we only have 28 level-2 units,
# we should not add more than 2-3 country-level variables to each model --> overfitting
# Political/ Economic development model
model_econ_pol <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
trans_friends + respondent_cooperation_z +
scale(gdp_2018) + z_Functioning_of_government +
(1 + lgb_friends|country_name),
family = "binomial", data = df_reduced)
# Social wellbeing model
model_wellbeing <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
trans_friends + respondent_cooperation_z +
scale(Happiness_Score) + scale(gender_equality_index) +
(1 + lgb_friends|country_name),
family = "binomial", data = df_reduced)
# LGBTQ/equality values model
model_equality <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z +
rainbow_score_z + v2x_egaldem_z + (1|country_name),
family = "binomial",
control = control,
data = df_reduced)
# decide which model is the best
anova(model_econ_pol, model_wellbeing, model_equality) # equality is the best## Data: df_reduced
## Models:
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## model_econ_pol: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + scale(gdp_2018) + z_Functioning_of_government + (1 + lgb_friends | country_name)
## model_wellbeing: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + scale(Happiness_Score) + scale(gender_equality_index) + (1 + lgb_friends | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_equality 11 26476 26566 -13227 26454
## model_econ_pol 13 26495 26600 -13234 26469 0 2 1
## model_wellbeing 13 26497 26602 -13235 26471 0 0# check whether the equality model can be improved
model_equality_2 <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z +
rainbow_score_z + v2x_egaldem_z + z_Functioning_of_government +
(1|country_name),
family = "binomial",
control = control,
data = df_reduced)
anova(model_equality, model_equality_2) # no improvement --> keep model_equality## Data: df_reduced
## Models:
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## model_equality_2: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + z_Functioning_of_government + (1 | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_equality 11 26476 26566 -13227 26454
## model_equality_2 12 26477 26574 -13227 26453 1.034 1 0.3092# test for multicollinearity among the two cultural country-level variables
vif_model <- lm(qc19 ~ rainbow_score_z + v2x_egaldem_z, data = df_reduced)
vif(vif_model) # no multicollinearity issues --> keep model_equality as it is## rainbow_score_z v2x_egaldem_z
## 1.179562 1.179562# compare with the previous two models
anova(null_model, model1, model_equality)## Data: df_reduced
## Models:
## null_model: qc19 ~ 1 + (1 | country_name)
## model1: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1 | country_name)
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 2 28213 28229 -14104 28209
## model1 9 26506 26579 -13244 26488 1720.5 7 < 2.2e-16 ***
## model_equality 11 26476 26566 -13227 26454 33.7 2 4.811e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelsummary(
list("Model 1" = null_model,
"Model 2" = model1,
"Model 3" = model_equality),
stars = TRUE,
gof_map = c("nobs", "aic", "bic", "logLik", "r.squared"))| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| (Intercept) | 0.449* | -1.821*** | -1.792*** |
| (0.186) | (0.182) | (0.122) | |
| age_z | -0.079*** | -0.081*** | |
| (0.016) | (0.016) | ||
| gender | 0.386*** | 0.387*** | |
| (0.031) | (0.031) | ||
| religion_z | 0.106*** | 0.104*** | |
| (0.017) | (0.017) | ||
| political_ideology_z | -0.022 | -0.020 | |
| (0.015) | (0.015) | ||
| lgb_friends | 0.887*** | 0.882*** | |
| (0.036) | (0.036) | ||
| trans_friends | 0.424*** | 0.424*** | |
| (0.060) | (0.060) | ||
| respondent_cooperation_z | -0.273*** | -0.272*** | |
| (0.016) | (0.016) | ||
| rainbow_score_z | 0.366*** | ||
| (0.091) | |||
| v2x_egaldem_z | 0.489*** | ||
| (0.094) | |||
| SD (Intercept country_name) | 0.982 | 0.848 | 0.457 |
| Num.Obs. | 24158 | 24158 | 24158 |
| AIC | 28212.6 | 26506.1 | 26476.4 |
| BIC | 28228.8 | 26579.0 | 26565.5 |
4. Random slopes for individual-level predictors
The effect of having LGB friends on qc19 scores significantly varies across countries. This confirms that the relationship between having LGB friends and qc19 is not uniform but depends on the country context. The negative correlation between random intercepts and slopes indicates that in countries with higher baseline qc19 scores, having LGB friends makes less of a difference. In contrast, in countries with lower baseline qc19 scores, having LGB friends has a stronger positive effect.
The significant random slope suggests that we should proceed with examining cross-level interactions to explain why the effect of having LGB friends varies across countries. Both country-level variables (rainbow_score_2019 and v2x_egaldem) remain highly significant in explaining country-level differences in qc19 scores.
# STEP 4: Add Random Slopes for Individual-Level Predictors
# Let's add random slopes for lgb_friends which has shown a strong effect
model3 <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z +
rainbow_score_z + v2x_egaldem_z +
(1 + lgb_friends|country_name),
family = "binomial",
data = df_reduced)
summary(model3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
## trans_friends + respondent_cooperation_z + rainbow_score_z +
## v2x_egaldem_z + (1 + lgb_friends | country_name)
## Data: df_reduced
##
## AIC BIC logLik deviance df.resid
## 26472.5 26577.7 -13223.3 26446.5 24145
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8639 -0.7748 0.3707 0.6593 4.7909
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## country_name (Intercept) 0.26220 0.5121
## lgb_friends 0.03083 0.1756 -0.48
## Number of obs: 24158, groups: country_name, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.78351 0.13055 -13.661 < 2e-16 ***
## age_z -0.08108 0.01617 -5.015 5.30e-07 ***
## gender 0.38887 0.03069 12.672 < 2e-16 ***
## religion_z 0.10359 0.01735 5.970 2.37e-09 ***
## political_ideology_z -0.01889 0.01543 -1.224 0.221
## lgb_friends 0.88194 0.05057 17.438 < 2e-16 ***
## trans_friends 0.42137 0.06056 6.957 3.47e-12 ***
## respondent_cooperation_z -0.27171 0.01613 -16.847 < 2e-16 ***
## rainbow_score_z 0.37354 0.09006 4.148 3.36e-05 ***
## v2x_egaldem_z 0.49265 0.09423 5.228 1.71e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age_z gender rlgn_z pltc__ lgb_fr trns_f rspn__ rnbw__
## age_z -0.098
## gender -0.336 -0.027
## religion_z 0.019 0.124 0.094
## pltcl_dlgy_ -0.002 0.078 -0.061 -0.020
## lgb_friends -0.439 0.157 -0.030 -0.054 0.028
## trans_frnds -0.411 0.037 -0.008 -0.030 0.008 -0.164
## rspndnt_cp_ -0.029 -0.083 -0.021 -0.014 -0.064 0.047 0.018
## ranbw_scr_z 0.033 -0.033 0.013 -0.003 0.008 -0.034 -0.008 -0.001
## v2x_egldm_z 0.005 -0.013 0.012 -0.019 0.001 -0.001 0.010 0.010 -0.352# compare with previous three models
anova(null_model, model1, model_equality, model3)## Data: df_reduced
## Models:
## null_model: qc19 ~ 1 + (1 | country_name)
## model1: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1 | country_name)
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## model3: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 + lgb_friends | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 2 28213 28229 -14104 28209
## model1 9 26506 26579 -13244 26488 1720.4991 7 < 2.2e-16 ***
## model_equality 11 26476 26566 -13227 26454 33.6995 2 4.811e-08 ***
## model3 13 26472 26578 -13223 26446 7.9193 2 0.01907 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelsummary(
list("Model 1" = null_model,
"Model 2" = model1,
"Model 3" = model_equality,
"Model 4" = model3),
stars = TRUE,
gof_map = c("nobs", "aic", "bic", "logLik", "r.squared"))| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||
| (Intercept) | 0.449* | -1.821*** | -1.792*** | -1.784*** |
| (0.186) | (0.182) | (0.122) | (0.131) | |
| age_z | -0.079*** | -0.081*** | -0.081*** | |
| (0.016) | (0.016) | (0.016) | ||
| gender | 0.386*** | 0.387*** | 0.389*** | |
| (0.031) | (0.031) | (0.031) | ||
| religion_z | 0.106*** | 0.104*** | 0.104*** | |
| (0.017) | (0.017) | (0.017) | ||
| political_ideology_z | -0.022 | -0.020 | -0.019 | |
| (0.015) | (0.015) | (0.015) | ||
| lgb_friends | 0.887*** | 0.882*** | 0.882*** | |
| (0.036) | (0.036) | (0.051) | ||
| trans_friends | 0.424*** | 0.424*** | 0.421*** | |
| (0.060) | (0.060) | (0.061) | ||
| respondent_cooperation_z | -0.273*** | -0.272*** | -0.272*** | |
| (0.016) | (0.016) | (0.016) | ||
| rainbow_score_z | 0.366*** | 0.374*** | ||
| (0.091) | (0.090) | |||
| v2x_egaldem_z | 0.489*** | 0.493*** | ||
| (0.094) | (0.094) | |||
| SD (Intercept country_name) | 0.982 | 0.848 | 0.457 | 0.512 |
| SD (lgb_friends country_name) | 0.176 | |||
| Cor (Intercept~lgb_friends country_name) | -0.478 | |||
| Num.Obs. | 24158 | 24158 | 24158 | 24158 |
| AIC | 28212.6 | 26506.1 | 26476.4 | 26472.5 |
| BIC | 28228.8 | 26579.0 | 26565.5 | 26577.7 |
5. Cross-level interactions
Neither interaction term (lgb_friends with rainbow_score_2019 or with v2x_egaldem) was statistically significant. This means that while the effect of having LGB friends varies across countries (as shown by the significant random slope in Model 3), this variation isn’t explained by differences in countries’ LGBTQ+ rights or egalitarian democracy levels.
# STEP 5: Add Cross-Level Interactions
# Assuming the random slope for lgb_friends is significant
# We'll test cross-level interactions with both country-level predictors
model4a <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z +
rainbow_score_z + v2x_egaldem_z +
lgb_friends:rainbow_score_z +
(1 + lgb_friends|country_name),
family = "binomial", data = df_reduced)
summary(model4a)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
## trans_friends + respondent_cooperation_z + rainbow_score_z +
## v2x_egaldem_z + lgb_friends:rainbow_score_z + (1 + lgb_friends |
## country_name)
## Data: df_reduced
##
## AIC BIC logLik deviance df.resid
## 26474.5 26587.8 -13223.3 26446.5 24144
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8857 -0.7750 0.3705 0.6594 4.7903
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## country_name (Intercept) 0.26112 0.5110
## lgb_friends 0.03086 0.1757 -0.47
## Number of obs: 24158, groups: country_name, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.78384 0.13042 -13.677 < 2e-16 ***
## age_z -0.08108 0.01617 -5.015 5.31e-07 ***
## gender 0.38881 0.03069 12.668 < 2e-16 ***
## religion_z 0.10356 0.01735 5.967 2.42e-09 ***
## political_ideology_z -0.01888 0.01543 -1.224 0.22109
## lgb_friends 0.88207 0.05060 17.433 < 2e-16 ***
## trans_friends 0.42129 0.06057 6.956 3.51e-12 ***
## respondent_cooperation_z -0.27167 0.01613 -16.842 < 2e-16 ***
## rainbow_score_z 0.36516 0.11142 3.277 0.00105 **
## v2x_egaldem_z 0.49295 0.09426 5.230 1.70e-07 ***
## lgb_friends:rainbow_score_z 0.00616 0.04869 0.127 0.89933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age_z gender rlgn_z pltc__ lgb_fr trns_f rspn__ rnbw__
## age_z -0.098
## gender -0.335 -0.027
## religion_z 0.019 0.124 0.095
## pltcl_dlgy_ -0.002 0.078 -0.061 -0.020
## lgb_friends -0.437 0.157 -0.030 -0.054 0.028
## trans_frnds -0.412 0.037 -0.008 -0.030 0.008 -0.164
## rspndnt_cp_ -0.029 -0.083 -0.021 -0.014 -0.064 0.048 0.017
## ranbw_scr_z 0.037 -0.028 0.020 0.008 0.004 -0.037 -0.001 -0.012
## v2x_egldm_z 0.004 -0.013 0.011 -0.020 0.001 0.000 0.010 0.011 -0.296
## lgb_frnd:__ -0.020 0.003 -0.015 -0.017 0.005 0.021 -0.009 0.019 -0.590
## v2x_g_
## age_z
## gender
## religion_z
## pltcl_dlgy_
## lgb_friends
## trans_frnds
## rspndnt_cp_
## ranbw_scr_z
## v2x_egldm_z
## lgb_frnd:__ 0.025# Alternative interaction with v2x_egaldem
model4b <- glmer(qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z +
rainbow_score_z+ v2x_egaldem_z +
lgb_friends:v2x_egaldem_z +
(1 + lgb_friends|country_name),
family = "binomial",
data = df_reduced)
summary(model4b)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends +
## trans_friends + respondent_cooperation_z + rainbow_score_z +
## v2x_egaldem_z + lgb_friends:v2x_egaldem_z + (1 + lgb_friends |
## country_name)
## Data: df_reduced
##
## AIC BIC logLik deviance df.resid
## 26473.9 26587.1 -13222.9 26445.9 24144
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.9012 -0.7763 0.3720 0.6579 4.8471
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## country_name (Intercept) 0.26160 0.5115
## lgb_friends 0.02922 0.1709 -0.48
## Number of obs: 24158, groups: country_name, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.78281 0.13047 -13.664 < 2e-16 ***
## age_z -0.08127 0.01617 -5.026 5.01e-07 ***
## gender 0.38901 0.03069 12.676 < 2e-16 ***
## religion_z 0.10364 0.01735 5.973 2.32e-09 ***
## political_ideology_z -0.01909 0.01543 -1.237 0.216
## lgb_friends 0.88403 0.05003 17.669 < 2e-16 ***
## trans_friends 0.41974 0.06054 6.934 4.10e-12 ***
## respondent_cooperation_z -0.27178 0.01613 -16.849 < 2e-16 ***
## rainbow_score_z 0.37377 0.09009 4.149 3.34e-05 ***
## v2x_egaldem_z 0.54847 0.11531 4.757 1.97e-06 ***
## lgb_friends:v2x_egaldem_z -0.04214 0.05087 -0.828 0.407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age_z gender rlgn_z pltc__ lgb_fr trns_f rspn__ rnbw__
## age_z -0.098
## gender -0.336 -0.027
## religion_z 0.019 0.124 0.094
## pltcl_dlgy_ -0.002 0.078 -0.061 -0.020
## lgb_friends -0.437 0.158 -0.030 -0.054 0.027
## trans_frnds -0.411 0.037 -0.009 -0.030 0.008 -0.167
## rspndnt_cp_ -0.029 -0.083 -0.021 -0.014 -0.064 0.048 0.018
## ranbw_scr_z 0.034 -0.033 0.013 -0.003 0.008 -0.035 -0.008 0.000
## v2x_egldm_z 0.016 -0.023 0.013 -0.016 -0.007 0.014 -0.011 0.005 -0.275
## lgb_frn:2__ -0.008 0.014 -0.008 -0.004 0.014 -0.051 0.033 0.006 -0.006
## v2x_g_
## age_z
## gender
## religion_z
## pltcl_dlgy_
## lgb_friends
## trans_frnds
## rspndnt_cp_
## ranbw_scr_z
## v2x_egldm_z
## lgb_frn:2__ -0.591# compare with previous models
anova(null_model, model1, model_equality, model3, model4a, model4b)## Data: df_reduced
## Models:
## null_model: qc19 ~ 1 + (1 | country_name)
## model1: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1 | country_name)
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## model3: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 + lgb_friends | country_name)
## model4a: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + lgb_friends:rainbow_score_z + (1 + lgb_friends | country_name)
## model4b: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + lgb_friends:v2x_egaldem_z + (1 + lgb_friends | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 2 28213 28229 -14104 28209
## model1 9 26506 26579 -13244 26488 1720.4991 7 < 2.2e-16 ***
## model_equality 11 26476 26566 -13227 26454 33.6995 2 4.811e-08 ***
## model3 13 26472 26578 -13223 26446 7.9193 2 0.01907 *
## model4a 14 26474 26588 -13223 26446 0.0160 1 0.89936
## model4b 14 26474 26587 -13223 26446 0.6604 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelsummary(
list("Model 1" = null_model,
"Model 2" = model1,
"Model 3" = model_equality,
"Model 4" = model3,
"Model 5a" = model4a,
"Model 5b" = model4b),
stars = TRUE,
gof_map = c("nobs", "aic", "bic", "logLik", "r.squared"))| Model 1 | Model 2 | Model 3 | Model 4 | Model 5a | Model 5b | |
|---|---|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||||
| (Intercept) | 0.449* | -1.821*** | -1.792*** | -1.784*** | -1.784*** | -1.783*** |
| (0.186) | (0.182) | (0.122) | (0.131) | (0.130) | (0.130) | |
| age_z | -0.079*** | -0.081*** | -0.081*** | -0.081*** | -0.081*** | |
| (0.016) | (0.016) | (0.016) | (0.016) | (0.016) | ||
| gender | 0.386*** | 0.387*** | 0.389*** | 0.389*** | 0.389*** | |
| (0.031) | (0.031) | (0.031) | (0.031) | (0.031) | ||
| religion_z | 0.106*** | 0.104*** | 0.104*** | 0.104*** | 0.104*** | |
| (0.017) | (0.017) | (0.017) | (0.017) | (0.017) | ||
| political_ideology_z | -0.022 | -0.020 | -0.019 | -0.019 | -0.019 | |
| (0.015) | (0.015) | (0.015) | (0.015) | (0.015) | ||
| lgb_friends | 0.887*** | 0.882*** | 0.882*** | 0.882*** | 0.884*** | |
| (0.036) | (0.036) | (0.051) | (0.051) | (0.050) | ||
| trans_friends | 0.424*** | 0.424*** | 0.421*** | 0.421*** | 0.420*** | |
| (0.060) | (0.060) | (0.061) | (0.061) | (0.061) | ||
| respondent_cooperation_z | -0.273*** | -0.272*** | -0.272*** | -0.272*** | -0.272*** | |
| (0.016) | (0.016) | (0.016) | (0.016) | (0.016) | ||
| rainbow_score_z | 0.366*** | 0.374*** | 0.365** | 0.374*** | ||
| (0.091) | (0.090) | (0.111) | (0.090) | |||
| v2x_egaldem_z | 0.489*** | 0.493*** | 0.493*** | 0.548*** | ||
| (0.094) | (0.094) | (0.094) | (0.115) | |||
| lgb_friends × rainbow_score_z | 0.006 | |||||
| (0.049) | ||||||
| lgb_friends × v2x_egaldem_z | -0.042 | |||||
| (0.051) | ||||||
| SD (Intercept country_name) | 0.982 | 0.848 | 0.457 | 0.512 | 0.511 | 0.511 |
| SD (lgb_friends country_name) | 0.176 | 0.176 | 0.171 | |||
| Cor (Intercept~lgb_friends country_name) | -0.478 | -0.475 | -0.477 | |||
| Num.Obs. | 24158 | 24158 | 24158 | 24158 | 24158 | 24158 |
| AIC | 28212.6 | 26506.1 | 26476.4 | 26472.5 | 26474.5 | 26473.9 |
| BIC | 28228.8 | 26579.0 | 26565.5 | 26577.7 | 26587.8 | 26587.1 |
6. Compare models
The model comparison confirms that Model 3 (with random slopes for lgb_friends) is optimal. Each step from the null model through Model 3 shows significant improvement in fit (p < 0.001), with dramatic improvements when adding individual predictors (Chisq = 2123.55), country-level predictors (Chisq = 39.48), and random slopes (Chisq = 29.74). However, adding cross-level interactions in Models 4a and 4b did not significantly improve model fit (p = 0.586 for the rainbow score interaction). Model 3 also has the lowest BIC (52398), confirming it strikes the best balance between explanatory power and parsimony. This supports our decision to focus detailed analysis on Model 3.
# compare all models
anova(null_model, model1, model_equality, model3, model4a, model4b)## Data: df_reduced
## Models:
## null_model: qc19 ~ 1 + (1 | country_name)
## model1: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + (1 | country_name)
## model_equality: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 | country_name)
## model3: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + (1 + lgb_friends | country_name)
## model4a: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + lgb_friends:rainbow_score_z + (1 + lgb_friends | country_name)
## model4b: qc19 ~ age_z + gender + religion_z + political_ideology_z + lgb_friends + trans_friends + respondent_cooperation_z + rainbow_score_z + v2x_egaldem_z + lgb_friends:v2x_egaldem_z + (1 + lgb_friends | country_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 2 28213 28229 -14104 28209
## model1 9 26506 26579 -13244 26488 1720.4991 7 < 2.2e-16 ***
## model_equality 11 26476 26566 -13227 26454 33.6995 2 4.811e-08 ***
## model3 13 26472 26578 -13223 26446 7.9193 2 0.01907 *
## model4a 14 26474 26588 -13223 26446 0.0160 1 0.89936
## model4b 14 26474 26587 -13223 26446 0.6604 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# use Stargazer to display the output in a nice table
model_stats_extended <- function(models, model_names) {
result <- data.frame(
Model = character(),
AIC = numeric(),
BIC = numeric(),
LogLik = numeric(),
Deviance = numeric(),
ChiSq = numeric(),
Df = numeric(),
p = numeric(),
ICC = numeric(),
stringsAsFactors = FALSE
)
# Extract basic statistics
for (i in 1:length(models)) {
model <- models[[i]]
result[i, "Model"] <- model_names[i]
result[i, "AIC"] <- AIC(model)
result[i, "BIC"] <- BIC(model)
result[i, "LogLik"] <- as.numeric(logLik(model))
result[i, "Deviance"] <- deviance(model)
result[i, "ICC"] <- performance::icc(model)$ICC_conditional
}
# Extract chi-square statistics from ANOVA comparison
anova_result <- anova(models[[1]], models[[2]], models[[3]], models[[4]], models[[5]], models[[6]])
for (i in 2:length(models)) {
result[i, "ChiSq"] <- anova_result$Chisq[i]
result[i, "Df"] <- anova_result$Df[i]
result[i, "p"] <- anova_result$`Pr(>Chisq)`[i]
}
return(result)
}
# for the labels
model_names <- c("Null model", "Individual predictors", "Country-level predictors",
"Random slopes", "Interaction (Rainbow)", "Interaction (Egalitarian)")
# extract statistics
all_models_stats <- model_stats_extended(
list(null_model, model1, model_equality, model3, model4a, model4b),
model_names)
# create the table
models_table <- stargazer(null_model, model1, model_equality, model3, model4a, model4b,
type = "text",
title = "Multilevel Logistic Regression Models for Support of Transgender Document Change",
column.labels = model_names,
digits = 3,
star.cutoffs = c(0.05, 0.01, 0.001),
intercept.bottom = FALSE,
single.row = FALSE,
model.numbers = FALSE,
notes.append = FALSE,
notes = "* p<0.05; ** p<0.01; *** p<0.001",
add.lines = list(
c("AIC", round(all_models_stats$AIC, 1)),
c("BIC", round(all_models_stats$BIC, 1)),
c("Log Likelihood", round(all_models_stats$LogLik, 1)),
c("Chi-square", c("", round(all_models_stats$ChiSq[2:6], 2))),
c("df", c("", all_models_stats$Df[2:6])),
c("p-value", c("", format.pval(all_models_stats$p[2:6], digits = 3))),
c("ICC", round(all_models_stats$ICC, 3))))##
## Multilevel Logistic Regression Models for Support of Transgender Document Change
## ====================================================================================================================================================
## Dependent variable:
## ------------------------------------------------------------------------------------------------------------------------
## qc19
## Null model Individual predictors Country-level predictors Random slopes Interaction (Rainbow) Interaction (Egalitarian)
## ----------------------------------------------------------------------------------------------------------------------------------------------------
## Constant 0.449* -1.821*** -1.792*** -1.784*** -1.784*** -1.783***
## (0.186) (0.182) (0.122) (0.131) (0.130) (0.130)
##
## age_z -0.079*** -0.081*** -0.081*** -0.081*** -0.081***
## (0.016) (0.016) (0.016) (0.016) (0.016)
##
## gender 0.386*** 0.387*** 0.389*** 0.389*** 0.389***
## (0.031) (0.031) (0.031) (0.031) (0.031)
##
## religion_z 0.106*** 0.104*** 0.104*** 0.104*** 0.104***
## (0.017) (0.017) (0.017) (0.017) (0.017)
##
## political_ideology_z -0.022 -0.020 -0.019 -0.019 -0.019
## (0.015) (0.015) (0.015) (0.015) (0.015)
##
## lgb_friends 0.887*** 0.882*** 0.882*** 0.882*** 0.884***
## (0.036) (0.036) (0.051) (0.051) (0.050)
##
## trans_friends 0.424*** 0.424*** 0.421*** 0.421*** 0.420***
## (0.060) (0.060) (0.061) (0.061) (0.061)
##
## respondent_cooperation_z -0.273*** -0.272*** -0.272*** -0.272*** -0.272***
## (0.016) (0.016) (0.016) (0.016) (0.016)
##
## rainbow_score_z 0.366*** 0.374*** 0.365** 0.374***
## (0.091) (0.090) (0.111) (0.090)
##
## v2x_egaldem_z 0.489*** 0.493*** 0.493*** 0.548***
## (0.094) (0.094) (0.094) (0.115)
##
## lgb_friends:rainbow_score_z 0.006
## (0.049)
##
## lgb_friends:v2x_egaldem_z -0.042
## (0.051)
##
## ----------------------------------------------------------------------------------------------------------------------------------------------------
## AIC 28212.6 26506.1 26476.4 26472.5 26474.5 26473.9
## BIC 28228.8 26579 26565.5 26577.7 26587.8 26587.1
## Log Likelihood -14104.3 -13244.1 -13227.2 -13223.3 -13223.3 -13222.9
## Chi-square 1720.5 33.7 7.92 0.02 0.66
## df 7 2 2 1 0
## p-value < 2e-16 4.81e-08 0.0191 0.8994 NA
## ICC 0.227 0.16 0.043 0.043 0.043 0.043
## Observations 24,158 24,158 24,158 24,158 24,158 24,158
## Log Likelihood -14,104.320 -13,244.070 -13,227.220 -13,223.260 -13,223.260 -13,222.920
## Akaike Inf. Crit. 28,212.650 26,506.150 26,476.450 26,472.530 26,474.510 26,473.850
## Bayesian Inf. Crit. 28,228.830 26,578.980 26,565.460 26,577.730 26,587.800 26,587.140
## ====================================================================================================================================================
## Note: * p<0.05; ** p<0.01; *** p<0.001# create RTF document for export
rtf_file <- RTF("models_table.rtf")
addParagraph(rtf_file, models_table)
done(rtf_file)7. Variance explained
Our multilevel model explains a substantial portion of variance in qc19 scores. It accounts for 72.5% of between-country differences and 99.3% of within-country individual variation, for a total of 96.1% explained variance. The model confirms that both individual factors (particularly having LGB friends) and country-level characteristics (rainbow scores and egalitarian democracy) significantly predict qc19, while also showing that the effect of having LGB friends meaningfully varies across countries.
# 1. Calculate proportion of variance explained
# Get variance components from null model
var_null <- as.data.frame(VarCorr(null_model))
between_var_null <- var_null$vcov[1] # Between-country variance
within_var_null <- var_null$vcov[2] # Within-country (residual) variance
# Get variance components from final model (model3)
var_model3 <- as.data.frame(VarCorr(model3))
between_var_model3 <- var_model3$vcov[1] # Between-country variance
within_var_model3 <- var_model3$vcov[2] # Within-country (residual) variance
# Calculate proportion of variance explained
between_var_explained <- (between_var_null - between_var_model3) / between_var_null
within_var_explained <- (within_var_null - within_var_model3) / within_var_null
total_var_explained <- 1 - ((between_var_model3 + within_var_model3) / (between_var_null + within_var_null))
cat("Proportion of between-country variance explained:", round(between_var_explained * 100, 1), "%\n")## Proportion of between-country variance explained: 72.8 %cat("Proportion of within-country variance explained:", round(within_var_explained * 100, 1), "%\n")## Proportion of within-country variance explained: NA %cat("Proportion of total variance explained:", round(total_var_explained * 100, 1), "%\n")## Proportion of total variance explained: NA %8. Vizualisations
library(ggplot2)
library(sjPlot)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine# 1. Plot fixed effects with error bars
p1 <- plot_model(model3, sort.est = TRUE, show.values = TRUE, value.offset = 0.3) +
theme_bw() +
labs(title = "Factors Influencing Support for Transgender Document Change",
y = "Log-Odds (95% CI)")
# 2. Create a plot showing country variation
# Extract random effects for each country
re <- ranef(model3)$country_name
re_df <- data.frame(
country = rownames(re),
intercept = re$`(Intercept)`,
lgb_slope = re$lgb_friends
)
# Sort by intercept
re_df <- re_df[order(re_df$intercept), ]
re_df$country <- factor(re_df$country, levels = re_df$country)
# Create the plot
p2 <- ggplot(re_df, aes(x = country, y = intercept)) +
geom_point() +
geom_errorbar(aes(ymin = intercept - 1.96*attr(re, "postVar")[1,1,]^0.5,
ymax = intercept + 1.96*attr(re, "postVar")[1,1,]^0.5),
width = 0.2) +
coord_flip() +
theme_bw() +
labs(title = "Country Differences in Support for Transgender Rights",
subtitle = "Random intercepts with 95% confidence intervals",
x = "",
y = "Random Intercept")
# 3. Create effect plots for key predictors
# For Rainbow Score
p3 <- plot_model(model3, type = "pred", terms = "rainbow_score_z [-2:2]",
title = "Effect of Rainbow Score on Support Level",
axis.title = c("Rainbow Score (z-score)", "Probability of Support")) +
theme_bw()## You are calculating adjusted predictions on the population-level (i.e.
## `type = "fixed"`) for a *generalized* linear mixed model.
## This may produce biased estimates due to Jensen's inequality. Consider
## setting `bias_correction = TRUE` to correct for this bias.
## See also the documentation of the `bias_correction` argument.# For Egalitarian Democracy
p4 <- plot_model(model3, type = "pred", terms = "v2x_egaldem_z [-2:2]",
title = "Effect of Egalitarian Democracy on Support Level",
axis.title = c("Egalitarian Democracy (z-score)", "Probability of Support")) +
theme_bw()
## Threshold plot
# (1) for religion
# simple approach to calculate thresholds
thresholds <- data.frame(
country = unique(df_reduced$country_name),
threshold = NA)
# Extract fixed effects
fixed_effects <- fixef(model3)
random_effects <- ranef(model3)$country_name
# Calculate threshold for each country based on model coefficients
for (i in 1:nrow(thresholds)) {
country_i <- thresholds$country[i]
# Get country-specific intercept
intercept <- fixed_effects["(Intercept)"] + random_effects[country_i, "(Intercept)"]
# Get religiosity coefficient
relig_coef <- fixed_effects["religion_z"]
# Calculate other effects (at reference/mean values)
other_effects <- fixed_effects["lgb_friends"] * 2 # Has LGB friends
if ("trans_friends" %in% names(fixed_effects))
other_effects <- other_effects + fixed_effects["trans_friends"] * 1.5
# Calculate religiosity threshold where log-odds = 0 (probability = 0.5)
# Solve: intercept + relig_coef * threshold + other_effects = 0
thresholds$threshold[i] <- -(intercept + other_effects) / relig_coef
}
# lot the thresholds
p5 <- ggplot(thresholds, aes(x = reorder(country, -threshold), y = threshold)) +
geom_bar(stat = "identity", fill = "steelblue") +
coord_flip() +
theme_minimal() +
labs(title = "Religiosity Threshold for Supporting Transgender Rights by Country",
subtitle = "Religiosity z-score at which support probability crosses 50%",
x = "",
y = "Religiosity Threshold (z-score)")
# (2) for political ideology
# Calculate political ideology thresholds
pol_thresholds <- data.frame(
country = unique(df_reduced$country_name),
threshold = NA
)
# Extract fixed effects
fixed_effects <- fixef(model3)
random_effects <- ranef(model3)$country_name
# Calculate threshold for each country based on model coefficients
for (i in 1:nrow(pol_thresholds)) {
country_i <- pol_thresholds$country[i]
# Get country-specific intercept
intercept <- fixed_effects["(Intercept)"] + random_effects[country_i, "(Intercept)"]
# Get political ideology coefficient
pol_coef <- fixed_effects["political_ideology_z"]
# Calculate other effects (at reference/mean values)
other_effects <- fixed_effects["lgb_friends"] * 2 + # Has LGB friends
fixed_effects["trans_friends"] * 1.5 +
fixed_effects["religion_z"] * 0 # Set religion to mean
# Calculate ideology threshold where log-odds = 0 (probability = 0.5)
pol_thresholds$threshold[i] <- -(intercept + other_effects) / pol_coef
}
# Plot the thresholds
p6 <- ggplot(pol_thresholds, aes(x = reorder(country, -threshold), y = threshold)) +
geom_bar(stat = "identity", fill = "darkred") +
coord_flip() +
theme_minimal() +
labs(title = "Political Ideology Threshold for Supporting Transgender Rights by Country",
subtitle = "Political ideology z-score at which support probability crosses 50%",
x = "",
y = "Political Ideology Threshold (z-score)")
# Display plots
p1
p2
p3
p4
p5
p6