Bayesian hereditarianism
21 February, 2025
Init
library(kirkegaard)## Loading required package: 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.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── 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 errors
## Loading required package: magrittr
##
##
## Attaching package: 'magrittr'
##
##
## The following object is masked from 'package:purrr':
##
## set_names
##
##
## The following object is masked from 'package:tidyr':
##
## extract
##
##
## Loading required package: weights
##
## Loading required package: Hmisc
##
##
## Attaching package: 'Hmisc'
##
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
##
## The following objects are masked from 'package:base':
##
## format.pval, units
##
##
## Loading required package: assertthat
##
##
## Attaching package: 'assertthat'
##
##
## The following object is masked from 'package:tibble':
##
## has_name
##
##
## Loading required package: psych
##
##
## Attaching package: 'psych'
##
##
## The following object is masked from 'package:Hmisc':
##
## describe
##
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
##
##
## Attaching package: 'kirkegaard'
##
##
## The following object is masked from 'package:psych':
##
## rescale
##
##
## The following object is masked from 'package:assertthat':
##
## are_equal
##
##
## The following object is masked from 'package:purrr':
##
## is_logical
##
##
## The following object is masked from 'package:base':
##
## +load_packages(
)
theme_set(theme_bw())
options(
digits = 3
)Functions
#plot distribution
plot_beta = function(shape1, shape2) {
dd = tibble(
x = seq(-1, 1, 0.001),
y = dbeta(x, shape1, shape2)
)
ggplot(dd, aes(x, y)) +
geom_density(stat = "identity") +
labs(
title = paste0("Beta(", shape1, ", ", shape2, ")"),
x = "x",
y = "Density"
)
}Analysis
#PGS models
pgs_range = seq(-1, 1, 0.001)
pgs_range_dx = diff(pgs_range)[1]
pgs_h_prior = pgs_h_prior <- dbeta((pgs_range + 1)/2, 8, 1.5) * 0.5 # Beta(2, 1) on [-1, 1]
pgs_e_prior = dnorm(pgs_range, 0, 0.1)
pgs_a_prior = dunif(pgs_range, -1, 1)
#normalize to 1
pgs_h_prior <- pgs_h_prior / (sum(pgs_h_prior) * pgs_range_dx)
pgs_e_prior <- pgs_e_prior / (sum(pgs_e_prior) * pgs_range_dx)
pgs_a_prior <- pgs_a_prior / (sum(pgs_a_prior) * pgs_range_dx)
#plot together
tibble(
x = pgs_range,
y = pgs_h_prior,
model = "Hereditarianism"
) %>%
bind_rows(
tibble(
x = pgs_range,
y = pgs_e_prior,
model = "Egalitarianism"
),
tibble(
x = pgs_range,
y = pgs_a_prior,
model = "Agnostic"
)
) %>%
ggplot(aes(x, y, color = model)) +
geom_line() +
labs(
title = "Prior beliefs about group-level PGS correlation",
x = "PGS",
y = "Density"
)
#save
GG_save("figs/pgs_priors.png")
#observations
pgs_obs = c(
0.89, #Piffer 2019, 1kg, figure 2, https://www.mdpi.com/2624-8611/1/1/5
0.98, #Piffer 2019, gnomad, figure 9, https://www.mdpi.com/2624-8611/1/1/5
0.90 #Fuerst 2023, abcd, figure 8a, 10.46469/mq.2023.63.4.2
)
# Likelihood function (scalar output for a given x)
likelihood <- function(x, obs) {
prod(dnorm(obs, mean = x, sd = 0.1))
}
# Compute likelihood across the grid
lik <- sapply(pgs_range, likelihood, obs = pgs_obs)
# Compute unnormalized posteriors
pgs_h_post_unnorm <- pgs_h_prior * lik
pgs_e_post_unnorm <- pgs_e_prior * lik
pgs_a_post_unnorm <- pgs_a_prior * lik
# Compute marginal likelihoods (evidence)
pgs_h_marg <- sum(pgs_h_prior * lik) * pgs_range_dx
pgs_e_marg <- sum(pgs_e_prior * lik) * pgs_range_dx
pgs_a_marg <- sum(pgs_a_prior * lik) * pgs_range_dx
# Normalize posteriors
pgs_h_post <- pgs_h_post_unnorm / pgs_h_marg
pgs_e_post <- pgs_e_post_unnorm / pgs_e_marg
pgs_a_post <- pgs_a_post_unnorm / pgs_a_marg
# Compute Bayes Factors (scalar ratios)
bf_eh <- pgs_e_marg / pgs_h_marg # Egalitarianism vs Hereditarianism
bf_ah <- pgs_a_marg / pgs_h_marg # Agnostic vs Hereditarianism
bf_ae <- pgs_a_marg / pgs_e_marg # Agnostic vs Egalitarianism
# Print marginal likelihoods and Bayes factors
cat("Marginal Likelihoods:\n")## Marginal Likelihoods:cat("Hereditarianism:", pgs_h_marg, "\n")## Hereditarianism: 12.2cat("Egalitarianism:", pgs_e_marg, "\n")## Egalitarianism: 3.25e-13cat("Agnostic:", pgs_a_marg, "\n")## Agnostic: 3.27cat("\nBayes Factors:\n")##
## Bayes Factors:cat("BF (Hereditarianism vs. Egalitarianism):", round(1/bf_eh, 2), "\n")## BF (Hereditarianism vs. Egalitarianism): 3.75e+13cat("BF (Hereditarianism vs. Agnostic):", round(1/bf_ah, 2), "\n")## BF (Hereditarianism vs. Agnostic): 3.72cat("BF (Agnostic vs Egalitarianism):", round(bf_ae, 2), "\n")## BF (Agnostic vs Egalitarianism): 1.01e+13# Posterior model probabilities (assuming equal prior model probabilities)
prior_model_prob <- 1/3
total_evidence <- pgs_h_marg + pgs_e_marg + pgs_a_marg
pgs_h_post_prob <- (pgs_h_marg * prior_model_prob) / total_evidence
pgs_e_post_prob <- (pgs_e_marg * prior_model_prob) / total_evidence
pgs_a_post_prob <- (pgs_a_marg * prior_model_prob) / total_evidence
cat("\nPosterior Model Probabilities (assuming equal prior probs):\n")##
## Posterior Model Probabilities (assuming equal prior probs):cat("Hereditarianism:", pgs_h_post_prob, "\n")## Hereditarianism: 0.263cat("Egalitarianism:", pgs_e_post_prob, "\n")## Egalitarianism: 7.01e-15cat("Agnostic:", pgs_a_post_prob, "\n")## Agnostic: 0.0707#model probabilities before and after
tibble(
model = c("Hereditarianism", "Egalitarianism", "Agnostic"),
prior = c(1/3, 1/3, 1/3),
posterior = c(pgs_h_post_prob, pgs_e_post_prob, pgs_a_post_prob) / sum(c(pgs_h_post_prob, pgs_e_post_prob, pgs_a_post_prob))
) %>%
pivot_longer(-model) %>%
mutate(
name = fct_relevel(name, "prior", "posterior")
) %>%
ggplot(aes(name, value, fill = model)) +
geom_col(position = "dodge") +
labs(
title = "Model probabilities",
x = "Before and after seeing evidence",
y = "Probability"
) +
scale_y_continuous(labels = scales::percent)
GG_save("figs/model_probs.png")
#plot posteriors together
bind_rows(
tibble(
x = pgs_range,
y = pgs_h_post,
model = "Hereditarianism"
),
tibble(
x = pgs_range,
y = pgs_e_post,
model = "Egalitarianism"
),
tibble(
x = pgs_range,
y = pgs_a_post,
model = "Agnostic"
)
) %>%
ggplot(aes(x, y, color = model)) +
geom_line() +
labs(
title = "Posterior beliefs about group-level PGS correlation",
x = "PGS",
y = "Density"
)
#save
GG_save("figs/pgs_posteriors.png")Meta
#versions
write_sessioninfo()## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 21.1
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_DK.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_DK.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Brussels
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kirkegaard_2025-01-08 psych_2.4.6.26 assertthat_0.2.1
## [4] weights_1.0.4 Hmisc_5.1-3 magrittr_2.0.3
## [7] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [10] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [13] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [16] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.2.0 digest_0.6.37
## [5] rpart_4.1.24 timechange_0.3.0 lifecycle_1.0.4 cluster_2.1.8
## [9] survival_3.8-3 gdata_3.0.0 compiler_4.4.2 rlang_1.1.4
## [13] sass_0.4.9 tools_4.4.2 utf8_1.2.4 yaml_2.3.10
## [17] data.table_1.16.0 knitr_1.48 labeling_0.4.3 htmlwidgets_1.6.4
## [21] mnormt_2.1.1 withr_3.0.1 foreign_0.8-88 nnet_7.3-20
## [25] grid_4.4.2 fansi_1.0.6 jomo_2.7-6 colorspace_2.1-1
## [29] mice_3.16.0 scales_1.3.0 gtools_3.9.5 iterators_1.0.14
## [33] MASS_7.3-64 cli_3.6.3 rmarkdown_2.28 ragg_1.3.2
## [37] generics_0.1.3 rstudioapi_0.16.0 tzdb_0.4.0 minqa_1.2.8
## [41] cachem_1.1.0 splines_4.4.2 parallel_4.4.2 base64enc_0.1-3
## [45] vctrs_0.6.5 boot_1.3-31 glmnet_4.1-8 Matrix_1.7-2
## [49] jsonlite_1.8.8 hms_1.1.3 mitml_0.4-5 Formula_1.2-5
## [53] htmlTable_2.4.3 systemfonts_1.1.0 foreach_1.5.2 jquerylib_0.1.4
## [57] glue_1.7.0 nloptr_2.1.1 pan_1.9 codetools_0.2-19
## [61] stringi_1.8.4 gtable_0.3.5 shape_1.4.6.1 lme4_1.1-35.5
## [65] munsell_0.5.1 pillar_1.9.0 htmltools_0.5.8.1 R6_2.5.1
## [69] textshaping_0.4.0 evaluate_0.24.0 lattice_0.22-5 highr_0.11
## [73] backports_1.5.0 broom_1.0.6 bslib_0.8.0 Rcpp_1.0.13
## [77] gridExtra_2.3 nlme_3.1-167 checkmate_2.3.2 xfun_0.47
## [81] pkgconfig_2.0.3# #write data to file for reuse
# d %>% write_rds("data/data_for_reuse.rds")
#
# #OSF
# if (F) {
# library(osfr)
#
# #login
# osf_auth(readr::read_lines("~/.config/osf_token"))
#
# #the project we will use
# osf_proj = osf_retrieve_node("https://osf.io/XXX/")
#
# #upload all files in project
# #overwrite existing (versioning)
# osf_upload(
# osf_proj,
# path = c("data", "figures", "papers", "notebook.Rmd", "notebook.html", "sessions_info.txt"),
# conflicts = "overwrite"
# )
# }