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())In this blind analysis, we are going to address 2 research questions:
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:
In order to best address the hypotheses, we need to decide a few things.
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"))marp_proc %>%
select(country,
gender, ses, education, ethnicity, denomination,
sample_type, compensation,
ends_with("_mean")) %>%
skim()| 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 | ▁▂▇▇▆ |
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 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
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 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
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
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")
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)
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)
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)
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 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")
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.
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)
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
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)
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)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:
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)
In the final dataset we remove the Japanese answers.
marp <-
marp_nodiag %>%
filter(country != "Japan")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 | ||
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 | ||
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 ' ' 1anova(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 ' ' 1anova(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