Probability of a girl birth given placenta previa (BDA3 p. 37)

Illustrate the effect of prior and compare posterior distributions with different parameter values for the beta prior distribution.

ggplot2 is used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
library(dplyr)

Observed data: 437 girls and 543 boys

a <- 437
b <- 543

Evaluate densities at evenly spaced points between 0.375 and 0.525

df1 <- data.frame(theta = seq(0.375, 0.525, 0.001))

Posterior with Beta(1,1), ie. uniform prior

df1$pu <- dbeta(df1$theta, a+1, b+1)

3 different choices for priors

n <- c(2, 20, 200) # prior counts
apr <- 0.488 # prior ratio of success

# helperf returns for given number of prior observations, prior ratio
# of successes, number of observed successes and failures and a data
# frame with values of theta, a new data frame with prior and posterior
# values evaluated at points theta.
helperf <- function(n, apr, a, b, df)
  cbind(df, pr = dbeta(df$theta, n*apr, n*(1-apr)), po = dbeta(df$theta, n*apr + a, n*(1-apr) + b), n = n)
# lapply function over prior counts n and pivot results into key-value pairs.
df2 <- lapply(n, helperf, apr, a, b, df1) %>%
  do.call(rbind, args = .) %>%
  pivot_longer(!c(theta, n), names_to = "grp", values_to = "p") %>%
  mutate(grp = factor(grp, labels=c('Posterior','Prior','Posterior with unif prior')))
# add correct labels for plotting
df2$title <- factor(paste0('alpha/(alpha+beta)=0.488, alpha+beta=',df2$n))

Plot distributions

ggplot(data = df2) +
  geom_line(aes(theta, p, color = grp)) +
  geom_vline(xintercept = 0.488, linetype = 'dotted') +
  facet_wrap(~title, ncol = 1) +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL) +
  theme(legend.position = 'bottom', legend.title = element_blank())

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDIuMiIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIFByb2JhYmlsaXR5IG9mIGEgZ2lybCBiaXJ0aCBnaXZlbiBwbGFjZW50YSBwcmV2aWEgKEJEQTMgcC4gMzcpCgpJbGx1c3RyYXRlIHRoZSBlZmZlY3Qgb2YgcHJpb3IgYW5kIGNvbXBhcmUgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbnMKd2l0aCBkaWZmZXJlbnQgcGFyYW1ldGVyIHZhbHVlcyBmb3IgdGhlIGJldGEgcHJpb3IgZGlzdHJpYnV0aW9uLgoKZ2dwbG90MiBpcyB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yIG1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkoZHBseXIpCmBgYAoKT2JzZXJ2ZWQgZGF0YTogNDM3IGdpcmxzIGFuZCA1NDMgYm95cwoKYGBge3IgfQphIDwtIDQzNwpiIDwtIDU0MwpgYGAKCkV2YWx1YXRlIGRlbnNpdGllcyBhdCBldmVubHkgc3BhY2VkIHBvaW50cyBiZXR3ZWVuIDAuMzc1IGFuZCAwLjUyNQoKYGBge3IgfQpkZjEgPC0gZGF0YS5mcmFtZSh0aGV0YSA9IHNlcSgwLjM3NSwgMC41MjUsIDAuMDAxKSkKYGBgCgpQb3N0ZXJpb3Igd2l0aCBCZXRhKDEsMSksIGllLiB1bmlmb3JtIHByaW9yCgpgYGB7ciB9CmRmMSRwdSA8LSBkYmV0YShkZjEkdGhldGEsIGErMSwgYisxKQpgYGAKCjMgZGlmZmVyZW50IGNob2ljZXMgZm9yIHByaW9ycwoKLSBCZXRhKDAuNDg4XCoyLCgxLTAuNDg4KVwqMikKLSBCZXRhKDAuNDg4XCoyMCwoMS0wLjQ4OClcKjIwKQotIEJldGEoMC40ODhcKjIwMCwoMS0wLjQ4OClcKjIwMCkKCmBgYHtyIH0KbiA8LSBjKDIsIDIwLCAyMDApICMgcHJpb3IgY291bnRzCmFwciA8LSAwLjQ4OCAjIHByaW9yIHJhdGlvIG9mIHN1Y2Nlc3MKCiMgaGVscGVyZiByZXR1cm5zIGZvciBnaXZlbiBudW1iZXIgb2YgcHJpb3Igb2JzZXJ2YXRpb25zLCBwcmlvciByYXRpbwojIG9mIHN1Y2Nlc3NlcywgbnVtYmVyIG9mIG9ic2VydmVkIHN1Y2Nlc3NlcyBhbmQgZmFpbHVyZXMgYW5kIGEgZGF0YQojIGZyYW1lIHdpdGggdmFsdWVzIG9mIHRoZXRhLCBhIG5ldyBkYXRhIGZyYW1lIHdpdGggcHJpb3IgYW5kIHBvc3RlcmlvcgojIHZhbHVlcyBldmFsdWF0ZWQgYXQgcG9pbnRzIHRoZXRhLgpoZWxwZXJmIDwtIGZ1bmN0aW9uKG4sIGFwciwgYSwgYiwgZGYpCiAgY2JpbmQoZGYsIHByID0gZGJldGEoZGYkdGhldGEsIG4qYXByLCBuKigxLWFwcikpLCBwbyA9IGRiZXRhKGRmJHRoZXRhLCBuKmFwciArIGEsIG4qKDEtYXByKSArIGIpLCBuID0gbikKIyBsYXBwbHkgZnVuY3Rpb24gb3ZlciBwcmlvciBjb3VudHMgbiBhbmQgcGl2b3QgcmVzdWx0cyBpbnRvIGtleS12YWx1ZSBwYWlycy4KZGYyIDwtIGxhcHBseShuLCBoZWxwZXJmLCBhcHIsIGEsIGIsIGRmMSkgJT4lCiAgZG8uY2FsbChyYmluZCwgYXJncyA9IC4pICU+JQogIHBpdm90X2xvbmdlcighYyh0aGV0YSwgbiksIG5hbWVzX3RvID0gImdycCIsIHZhbHVlc190byA9ICJwIikgJT4lCiAgbXV0YXRlKGdycCA9IGZhY3RvcihncnAsIGxhYmVscz1jKCdQb3N0ZXJpb3InLCdQcmlvcicsJ1Bvc3RlcmlvciB3aXRoIHVuaWYgcHJpb3InKSkpCiMgYWRkIGNvcnJlY3QgbGFiZWxzIGZvciBwbG90dGluZwpkZjIkdGl0bGUgPC0gZmFjdG9yKHBhc3RlMCgnYWxwaGEvKGFscGhhK2JldGEpPTAuNDg4LCBhbHBoYStiZXRhPScsZGYyJG4pKQpgYGAKClBsb3QgZGlzdHJpYnV0aW9ucwoKYGBge3IgfQpnZ3Bsb3QoZGF0YSA9IGRmMikgKwogIGdlb21fbGluZShhZXModGhldGEsIHAsIGNvbG9yID0gZ3JwKSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAuNDg4LCBsaW5ldHlwZSA9ICdkb3R0ZWQnKSArCiAgZmFjZXRfd3JhcCh+dGl0bGUsIG5jb2wgPSAxKSArCiAgbGFicyh4ID0gJycsIHkgPSAnJykgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSBOVUxMKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgo=