MARP exploration - Team 40b

install.packages(c("tidyverse", "psych", "skimr", "ggridges", "tidytext", "lme4", "performance", "sjPlot", "here", "emo", "broom.mixed"))

library(tidyverse)
library(psych)
library(skimr)
library(ggridges)
library(tidytext)
library(lme4)
library(performance)
library(sjPlot)
library(broom.mixed)

theme_set(theme_light())

Main questions

In this blind analysis, we are going to address 2 research questions:

1. Do religious people report higher well-being?
2. Does the relation between religiosity and well-being depend on how important people consider religion to be in their country (i.e., perceived cultural norms of religion)?

Blinded data

The point of blinded data is to let analysts explore the data without the danger of p-hacking. Therefore, only the relationship of the outcomes and predctor variables are destroyed by shuffling, and the remainder is kept intact. This means that:

• Data reduction techinques (PCA, EFA) will yield valid results.
• Outiers, missing data can be observed and treated.
• Confounders can be explored, meaning that not only univariate distributions but correlations are also kept intact.
• Country level means should remain intact

Analysis issues to resolve

In order to best address the hypotheses, we need to decide a few things.

• Operationalization of variables: How should we conceptualize key variables?
• Outliers how should we handle them?
• Choose statistical model (multilevel?)
  • Confounders (which ones to use?)
  • Moderators
  • Lumping levels for nominal variables

Read and pre-process data

At this point, we only exclude participants who did not pass the attention check.

# Read raw data
marp_raw <- 
  # read_csv(here::here("data/MARP_data_blinded.csv"))
  read_csv(here::here("data/MARP_data.csv"))

## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   country = col_character(),
##   gender = col_character(),
##   ethnicity = col_character(),
##   denomination = col_character(),
##   sample_type = col_character(),
##   compensation = col_character()
## )
## i Use `spec()` for the full column specifications.

marp_proc <-
  marp_raw %>% 
  filter(attention_check == 1) %>% 
  mutate( gender = recode(gender, 
                          "man" = "Male", 
                          "woman" = "Female",
                          "other" = "Other"))

Descriptives and potential moderators

marp_proc %>% 
        select(country, 
               gender, ses, education, ethnicity, denomination,
               sample_type, compensation, 
               ends_with("_mean")) %>% 
        skim()

Data summary

Name	Piped data
Number of rows	10195
Number of columns	12
_______________________
Column type frequency:
character	6
numeric	6
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
country	0	1.00	2	11	0	24	0
gender	0	1.00	4	6	0	3	0
ethnicity	405	0.96	5	22	0	17	0
denomination	5567	0.45	4	41	0	20	0
sample_type	0	1.00	5	14	0	4	0
compensation	0	1.00	6	31	0	5	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ses	5	1	6.10	1.77	1.00	5.00	6.00	7.00	10	▁▃▇▇▁
education	0	1	4.64	1.26	1.00	4.00	5.00	5.00	7	▁▃▆▇▆
wb_overall_mean	0	1	3.67	0.61	1.22	3.28	3.78	4.11	5	▁▁▃▇▂
wb_phys_mean	0	1	3.84	0.66	1.00	3.43	4.00	4.29	5	▁▁▃▇▆
wb_psych_mean	0	1	3.51	0.72	1.00	3.00	3.67	4.00	5	▁▂▅▇▂
wb_soc_mean	0	1	3.56	0.87	1.00	3.00	3.67	4.17	5	▁▂▇▇▆

Demographics by country

N of participants

marp_proc %>% 
        mutate(country = fct_infreq(country) %>% fct_rev()) %>% 
        ggplot() +
        aes(y = country) +
        geom_bar() +
        scale_x_continuous(breaks = seq(0, 1500, 250)) +
        labs(title = "Number of participants by country",
             y = NULL)

Age

Age has a few obvious outliers (e.g. age of 0 or 1), that we could remove. However, the sample size makes it unlikely that these outliers have a large influence on the model, therefore we won’t remove them.

marp_proc %>% 
  drop_na(age) %>% 
  mutate(country = fct_reorder(country, age)) %>%
  ggplot() +
  aes(y = country, x = age, fill = country) +
  geom_boxplot(alpha = .5, width = .2, 
               outlier.alpha = .2, show.legend = FALSE) +
  geom_density_ridges(alpha = .5, show.legend = FALSE) +
  labs(title = "Age distribution by country",
       y = NULL)

## Picking joint bandwidth of 2.45

Gender

marp_proc %>% 
        group_by(country) %>% 
        summarise(pct_female = mean(gender == "Female", 
                                    na.rm = TRUE)) %>% 
        mutate(country = fct_reorder(country, pct_female)) %>%
        ggplot() +
        aes(y = country, x = pct_female) +
        geom_point() +
        scale_x_continuous(limits = c(0,1),
                           labels = scales::percent_format()) +
        labs(title = "Proportion of females by country",
             x = NULL, y = NULL)

## `summarise()` ungrouping output (override with `.groups` argument)

Education

Education is a common confounder when the relationship between religiosity and well-being is considered. There is a considerable spread of education within and between countries.

marp_proc %>% 
  mutate(country = fct_reorder(country, education, median)) %>%
  ggplot() +
  aes(y = country, x = education, fill = country) +
  geom_density_ridges(show.legend = FALSE) +
  labs(title = "Education by country",
       y = NULL, x = NULL)

## Picking joint bandwidth of 0.288

SES

Socio-economic status is a common confounder when the relationship between religiosity and well-being is considered. There is a considerable spread of SES within and between countries.

marp_proc %>% 
  drop_na(ses) %>% 
  mutate(country = fct_reorder(country, ses)) %>%
  ggplot() +
  aes(y = country, x = ses, fill = country) +
  geom_density_ridges(show.legend = FALSE) +
  labs(title = "Socio-economic status by country",
       y = NULL, x = NULL)

## Picking joint bandwidth of 0.403

Denomination

Religious denomination is often used when the connection between religiosity and well being is considered.

marp_proc %>% 
  mutate(country = fct_infreq(country) %>% fct_rev()) %>%
  count(country, denomination) %>%
  ggplot() +
  aes(y = country, x = n, fill = denomination) +
  geom_col(position = "stack") +
  labs(title = "Number/Proportion of denominations by country",
       y = NULL,
       fill = "Denomination")

However, in its current form, the variable contains categories that are ratther sparse, and this may interfere with the statistical model that we want to use. Therefore we choose to lump together levels that constitute less than 1% of the categories into “Other”.

lumped_denom <-
  marp_proc %>% 
  transmute(subject, 
            denom_lump = case_when(str_detect(denomination, "Muslim") ~ "Muslim",
                                   str_detect(denomination, "Christian|Evangelical") ~ "Christian",
                                   str_detect(denomination, "Other") ~ "Other",
                                   is.na(denomination) ~ "No denomination",
                                   TRUE ~ denomination) %>% 
                          fct_lump_prop(.01, other_level = "Other"))

marp_proc %>% 
  left_join(lumped_denom, by = "subject") %>% 
  mutate(country = fct_infreq(country) %>% fct_rev()) %>%
  count(country, denom_lump) %>%
  ggplot() +
  aes(y = country, x = n, fill = denom_lump, label = denom_lump) +
  geom_col(position = "stack") +
  labs(title = "Number/Proportion of denominations by country",
       subtitle = "Denominations that were infrequent (<1%) were lumped together",
       y = NULL,
       fill = "Denomination")

Ethnicity

Well-being might be influenced by ethicity, but not directly. For e.g. minority status might be associated with WB, but there were no questions about this. Therefore, in my opinion, raw ethnicity should not be added to the model.

marp_proc %>% 
  mutate(ethnicity = fct_infreq(ethnicity) %>% fct_rev()) %>% 
  ggplot() +
  aes(y = ethnicity, fill = ethnicity) +
  geom_bar(show.legend = FALSE) +
  labs(y = NULL)

Importance of religion

cnorm_questions <- 
  tibble(name = c("cnorm_1", "cnorm_2"),
         question = c("Importance of religious lifestyle for average person in country",
                      "Importance of belief in God/Gods for average person in country"
                    ))

marp_proc %>% 
  select(subject, country, cnorm_1, cnorm_2) %>% 
  pivot_longer(c("cnorm_1", "cnorm_2")) %>% 
  group_by(country, name) %>% 
  summarise(avg = mean(value)) %>% 
  ungroup() %>% 
  left_join(cnorm_questions, by = "name") %>%
  mutate(country = reorder_within(country, avg, question)) %>% 
  ggplot() +
  aes(x = avg, y = country) +
  geom_point() +
  scale_y_reordered(NULL) +
  facet_wrap(~question, ncol = 2, scales = "free_y") +
  labs(title = "Normativity of religion by country in two questions")  

## `summarise()` regrouping output by 'country' (override with `.groups` argument)

GDP per capita

marp_proc %>% 
  group_by(country) %>% 
  summarise(gdp = mean(gdp)) %>% 
  ungroup() %>% 
  mutate(country = fct_reorder(country, gdp)) %>% 
  ggplot() +
  aes(x = gdp, y = country) +
  geom_point() +
  scale_x_continuous(labels = scales::dollar_format()) +
  labs(title = "GDP per capita by country")  

## `summarise()` ungrouping output (override with `.groups` argument)

Sample type

marp_proc %>% 
  mutate(country = fct_infreq(country) %>% fct_rev()) %>%
  count(country, sample_type) %>% 
  ggplot() +
  aes(y = country, x = n, fill = sample_type) +
  geom_col(position = "stack") +
  labs(title = "Proportion of sample type by country",
       x = NULL, y = NULL, fill = "Sample type")

Compensation

Compensation is very closely associated with the sample type, therefore using it in the model would be redundant.

marp_proc %>% 
  mutate(country = fct_infreq(country) %>% fct_rev()) %>%
  count(country, compensation) %>% 
  ggplot() +
  aes(y = country, x = n, fill = compensation) +
  geom_col(position = "stack") +
  labs(title = "Proportion of compensation by country",
       x = NULL, y = NULL, fill = "Compensation")

Operationalization of variables

Religiosity and well-being has multiple items that we can use. According to previous studies of similar topics, norms about religion should be aggregated to the country level.

Religiosity

Although religiosity is an elusive concept, and no one-size-fits-all metric is available. We don’t feel competent to choose just one question, so We try to use as much information from all available questions as possible. I’m also not feeling confident to relevel specific questions (e.g. rel_3). Therefore, We choose to use PCA to extract an aggregated variable.

marp_proc %>% 
  select(starts_with("rel_")) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(value)) +
  geom_histogram() +
  facet_wrap(~name) +
  labs(title = "Distribution of rel_ variables to be used in PCA")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Number of components
marp_proc %>% 
        select(starts_with("rel_")) %>% 
        fa.parallel()

## Parallel analysis suggests that the number of factors =  4  and the number of components =  1

# Parallel analysis suggests to use a single component.
rel_pca <-
        marp_proc %>% 
        select(starts_with("rel_")) %>% 
        pca(nfactors = 1)

# Correlation of religiosity items and the PCA component
marp_proc %>% 
  select(starts_with("rel")) %>% 
  mutate(rel_pca = rel_pca$scores[,1]) %>% 
  cor.plot()

The religiosity values seems to vary considerably by country.

  marp_proc %>%
  mutate(rel_pca = rel_pca$scores[,1]) %>% 
  mutate(country = fct_reorder(country, rel_pca, median)) %>%
  ggplot() +
  aes(y = country, x = rel_pca, fill = country) +
  geom_boxplot(outlier.alpha = .2, show.legend = FALSE) +
  labs(title = "Religiosity by country",
       subtitle = "Aggregated religiosity component",
       x = NULL,
       y = NULL)

marp_proc %>%
  bind_cols(rel_pca = rel_pca$scores[,1]) %>% 
  left_join(lumped_denom, by = "subject") %>% 
  mutate(denom_lump = fct_reorder(denom_lump, rel_pca, median)) %>%
  ggplot() +
  aes(y = denom_lump, x = rel_pca, fill = denom_lump) +
  geom_boxplot(alpha = .5,
               width = .2,
               outlier.alpha = .2, show.legend = FALSE) +
  geom_density_ridges(alpha = .5) +
  labs(title = "Religiosity by denomination",
       subtitle = "Aggregated religiosity component",
       fill = "Denomination",
       x = NULL,
       y = NULL)

## Picking joint bandwidth of 0.173

We compared the PCA operationalization with the self-admitted single item religiosity. The difference on the country level seems quite subtle.

marp_proc %>% 
  mutate(rel_pca = rel_pca$scores[,1]) %>% 
  group_by(country) %>% 
  summarise(rel_item = mean(rel_3 == 1),
                  rel_pca = mean(rel_pca),
                  n = n()) %>% 
  mutate(across(starts_with("rel_"), ~scale(.x) %>% as.numeric())) %>% 
  pivot_longer(cols = c("rel_item", "rel_pca")) %>% 
  mutate(country = fct_reorder(country, value)) %>% 
  ggplot() +
  aes(x = value, y = country, color = name) +
  geom_point(size = 2) +
  labs(title = "Different operationalizations of religiosity lead to similar country-wise values (r = .91)",
       subtitle = "One item religiosity (rel_3) vs. 9-item PCA religiosity component (rel_pca) values",
       y = NULL, x = NULL, color = "Operationalization")

## `summarise()` ungrouping output (override with `.groups` argument)

Well-being

The well-being questions have calculated values for subscales and overall

marp_proc %>% 
  select(starts_with("wb")) %>% 
  cor.plot()

marp_proc %>% 
  mutate(country = fct_reorder(country, wb_overall_mean)) %>% 
  ggplot() +
  aes(x = wb_overall_mean, y = country, fill = country) +
  geom_boxplot(alpha = .5,
           width = .2,
           outlier.alpha = .2, show.legend = FALSE) +
  geom_density_ridges(alpha = .5, show.legend = FALSE) +
  labs(title = "Overall well-being by country",
       y = NULL)

## Picking joint bandwidth of 0.153

Cultural norms about religiosity into one variable.

According to previous research, norms should be handled on the country level. We use the mean of the two variables that are otherwise strongly correlated, since it makes no sense to use a more complicated method (PCA or EFA). Then we calculate the country level mean from this variable.

# Country norm of religiosity
country_norms <-
  marp_proc %>% 
  # Create new variable with the average of the 2 cnorm vairalbes
  mutate(cnorm_avg = ((cnorm_1 + cnorm_2)/2)) %>% 
  # Calculate country level norms
  group_by(country) %>% 
  summarise(cnorm_mean = mean(cnorm_avg)) %>% 
  ungroup() %>% 
  # Standardize the variable
  mutate(cnorm_mean = scale(cnorm_mean) %>% as.numeric())

## `summarise()` ungrouping output (override with `.groups` argument)

country_norms %>% 
  mutate(country = fct_reorder(country, cnorm_mean)) %>% 
  ggplot() +
  aes(y = country, x = cnorm_mean) +
  geom_point() +
  labs(title = "Average of cultural norms about religion by country",
       x = NULL, y = NULL)

Construct final dataset

Using all information from the exploratory analysis, we create a dataset for modeling. This dataset still doesn’t contain potential problems that may emerge during model diagnostics.

We add the religiosity component, the country-wise norms, the lumped denomination data, and set baselines for categorical variables. We also drop participants with missing values in varibles that we want to use in the statistical models, as those can cause difficulties when comparing models. This means dropping 25 participants.

marp_nodiag <-
  marp_proc %>% 
  # Add religiosity scores from PCA
  bind_cols(religiosity = rel_pca$scores[,1]) %>%
  # Add country level norms
  left_join(country_norms, by = "country") %>% 
  # Merge different branches of the same religion, lump levels < 1%
  left_join(lumped_denom, by = "subject") %>% 
  # Set baselines
  mutate(sample_type = fct_relevel(sample_type, "online panel"),
         denom_lump = fct_relevel(denom_lump, "No denomination"),
         gender = fct_relevel(gender, "Female")) %>% 
  # Drop participants with missing variables
  drop_na(age,  gender, ses, education, denom_lump, sample_type)

Investigating model assumptions

Before creating the final dataset and models, we investigate if there is anything strange in the model diagnostics that would necessitate further changes in the dataset. Therefore we create a model that contains all the terms that we want to include in the analysis, and we check all assumptions.

model_diag <- 
  lmer(wb_overall_mean ~ religiosity * cnorm_mean + 
         # personal level confounders
                         age + gender + ses + education + denom_lump +
         # country and sample level confounders
                         gdp_scaled + sample_type + 
         # random intercept and slope model
                         (religiosity|country), data = marp_nodiag)

check_model(model_diag)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 10170 rows containing missing values (geom_text_repel).

## `geom_smooth()` using formula 'y ~ x'

Model diagnostics show:

• ✅ No multicollinearity,
• ✅ Normally distributed residuals
• ✅ No influential cases
• ✅ Normally distributed random effects
• ❌ Homoskedasticity

Apart from heteroscedasticity, it seems like there is a strange separation in the fitted values. All residuals on the left hand side come from the Japanese sample. As the separation is complete and the difference is huge, we should handle the Japanese data with extra care. Further, there is very small variability in the Japanese fitted values.
Taken together, we decided to remove the Japanese data.

augment(model_diag) %>% 
  mutate(country = fct_reorder(country, .fitted)) %>% 
  ggplot() +
  aes(x = .fitted, y = country, fill = country) +
  geom_boxplot(show.legend = FALSE) +
  labs(title = "The fitted Japanese values are much lower than for any other country",
       y = NULL)

Correcting the final dataset

In the final dataset we remove the Japanese answers.

marp <- 
  marp_nodiag %>% 
  filter(country != "Japan")

Building models

1) Do religious people report higher well-being?

h1 <- 
  lmer(wb_overall_mean ~ religiosity + 
         # personal level confounders
                         age + gender + ses + education + denom_lump +
         # country and sample level confounders
                         gdp_scaled + sample_type +
         # random intercept and slope model
                         (religiosity|country), 
       data = marp)

# Create a null model for comparisons that does not contain the main predictor
h0 <- update(h1, . ~ . -religiosity)

check_model(h1)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 9747 rows containing missing values (geom_text_repel).

## `geom_smooth()` using formula 'y ~ x'

We can handle heteroscedasticity by using cluster robust standard errors (CR2), using the clubSandwich package. https://strengejacke.github.io/sjPlot/articles/tab_model_robust.html

tab_model(h1, 
          show.aic = TRUE, 
          show.reflvl = TRUE, 
          string.ci = "95% CI", 
          vcov.fun = "CR", 
          vcov.type = "CR2", 
          vcov.args = list(cluster = h1@frame$country))

	wb_overall_mean
Predictors	Estimates	95% CI	p
(Intercept)	2.68	2.54 – 2.82	<0.001
age	0.00	-0.00 – 0.00	0.076
Female	Reference
Male	0.04	0.01 – 0.08	0.011
Other	-0.29	-0.44 – -0.13	<0.001
No denomination	Reference
Buddhist	0.02	-0.05 – 0.10	0.567
Christian	-0.02	-0.07 – 0.03	0.422
Hindu	-0.08	-0.16 – -0.01	0.027
education	0.04	0.02 – 0.05	<0.001
gdp_scaled	0.03	-0.01 – 0.07	0.204
Jewish	-0.08	-0.15 – -0.00	0.049
Muslim	-0.09	-0.17 – -0.01	0.035
Other	-0.09	-0.19 – 0.01	0.072
religiosity	0.08	0.06 – 0.10	<0.001
online panel	Reference
general public	0.07	0.04 – 0.11	<0.001
mixed	0.09	-0.00 – 0.19	0.051
students	0.14	-0.10 – 0.38	0.248
ses	0.12	0.10 – 0.13	<0.001
Random Effects
σ2	0.28
τ00 country	0.01
τ11 country.religiosity	0.00
ρ01 country	-0.27
ICC	0.03
N country	23
Observations	9747
Marginal R2 / Conditional R2	0.171 / 0.200
AIC	15522.052

2) Does the relation between religiosity and well-being depend on how important people consider religion to be in their country (i.e., perceived cultural norms of religion)?

h2 <- 
  lmer(wb_overall_mean ~ religiosity * cnorm_mean + 
         # personal level confounders
                         age + gender + ses + education + denom_lump +
         # country level confounders
                         gdp_scaled + sample_type + 
         # random intercept and slope model
                         (religiosity|country), 
       data = marp)

check_model(h2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 9747 rows containing missing values (geom_text_repel).

## `geom_smooth()` using formula 'y ~ x'

Model diagnostics show heteroscedasticity, therefore cluster robust standard errors are calculated.

tab_model(h2, 
          show.aic = TRUE, 
          show.reflvl = TRUE, 
          string.ci = "95% CI", 
          vcov.fun = "CR", 
          vcov.type = "CR2", 
          vcov.args = list(cluster = h1@frame$country))

	wb_overall_mean
Predictors	Estimates	95% CI	p
(Intercept)	2.67	2.53 – 2.81	<0.001
age	0.00	-0.00 – 0.00	0.067
cnorm_mean	-0.03	-0.07 – 0.01	0.133
Female	Reference
Male	0.04	0.01 – 0.08	0.011
Other	-0.29	-0.44 – -0.13	<0.001
No denomination	Reference
Buddhist	0.03	-0.04 – 0.10	0.465
Christian	-0.02	-0.07 – 0.03	0.483
Hindu	-0.09	-0.17 – -0.01	0.030
education	0.04	0.03 – 0.05	<0.001
gdp_scaled	0.02	-0.03 – 0.07	0.526
Jewish	-0.08	-0.15 – -0.00	0.039
Muslim	-0.10	-0.17 – -0.04	0.002
Other	-0.09	-0.19 – 0.00	0.063
religiosity	0.08	0.06 – 0.09	<0.001
religiosity:cnorm_mean	0.03	0.02 – 0.04	<0.001
online panel	Reference
general public	0.07	0.04 – 0.11	<0.001
mixed	0.09	-0.00 – 0.19	0.053
students	0.15	-0.08 – 0.39	0.207
ses	0.12	0.10 – 0.13	<0.001
Random Effects
σ2	0.28
τ00 country	0.01
τ11 country.religiosity	0.00
ρ01 country	0.09
ICC	0.03
N country	23
Observations	9747
Marginal R2 / Conditional R2	0.173 / 0.200
AIC	15520.787

Model comparisons and Bayes Factors

anova(h0, h1)

## refitting model(s) with ML (instead of REML)

## Data: marp
## Models:
## h0: wb_overall_mean ~ age + gender + ses + education + denom_lump + 
## h0:     gdp_scaled + sample_type + (religiosity | country)
## h1: wb_overall_mean ~ religiosity + age + gender + ses + education + 
## h1:     denom_lump + gdp_scaled + sample_type + (religiosity | country)
##    npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## h0   20 15448 15592 -7704.2    15408                         
## h1   21 15418 15569 -7688.0    15376 32.335  1  1.297e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(h1, h2)

## refitting model(s) with ML (instead of REML)

## Data: marp
## Models:
## h1: wb_overall_mean ~ religiosity + age + gender + ses + education + 
## h1:     denom_lump + gdp_scaled + sample_type + (religiosity | country)
## h2: wb_overall_mean ~ religiosity * cnorm_mean + age + gender + ses + 
## h2:     education + denom_lump + gdp_scaled + sample_type + (religiosity | 
## h2:     country)
##    npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## h1   21 15418 15569  -7688    15376                         
## h2   23 15402 15567  -7678    15356 19.971  2  4.606e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(h0, h2)

## refitting model(s) with ML (instead of REML)

## Data: marp
## Models:
## h0: wb_overall_mean ~ age + gender + ses + education + denom_lump + 
## h0:     gdp_scaled + sample_type + (religiosity | country)
## h2: wb_overall_mean ~ religiosity * cnorm_mean + age + gender + ses + 
## h2:     education + denom_lump + gdp_scaled + sample_type + (religiosity | 
## h2:     country)
##    npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## h0   20 15448 15592 -7704.2    15408                         
## h2   23 15402 15567 -7678.0    15356 52.306  3  2.577e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Calculate BIC based Bayes factors for 
# H0 vs H1
exp((BIC(h0) - BIC(h1))/2)

## [1] 2127.347

# H1 vs H2
exp((BIC(h1) - BIC(h2))/2)

## [1] 0.00142744

# H0 vs H2
exp((BIC(h0) - BIC(h2))/2)

## [1] 3.03666