Chapter 2
Choosing causal graphs and priors

Lectures 3–4 · where graphs and priors come from, how to learn a graph from data, and how to assume as little as possible

Chapter 1 asked you to begin every analysis with two ingredients: a causal graph, saying what can affect what, and a set of priors, saying what is plausible before the data arrive. It quietly assumed you already had both. Usually you do, at least in part. An epidemiologist draws the graph she has taught for twenty years; a field biologist knows that penguins weigh a few kilograms. Most of the time, both ingredients come from domain knowledge, and the job is to write that knowledge down honestly.

But what if you know little? You may be handed a dataset with a dozen variables and no consensus diagram in any textbook. You may need a prior for a parameter whose plausible size you have never thought about. Making something up is acceptable if the something-up is principled, and this chapter is about the principles. It has two halves.

• Lecture 3 is about graphs. Every graph makes predictions that data can check. When knowledge is thin, some of the graph can be learned from data by causal discovery algorithms, and we will run them on real measurements from thousands of human immune cells. We will also see, sharply, what observational data can never tell you, and what experiments add.
• Lecture 4 is about priors. We will watch the choice of prior decide the answer in small real datasets, learn the principle of maximum entropy for assuming as little as possible, see why every prior lives on a particular scale, and set out what a publishable paper must report about its priors.

Lecture 3 — Where graphs come from

2.1 Where graphs come from

A causal graph is a summary of scientific knowledge, so it comes from wherever scientific knowledge comes from. In practice, four sources do most of the work.

1. Time order. Causes precede effects. A variable measured at birth cannot be caused by one measured at age forty. Time order alone orients a surprising number of arrows, and it is the most reliable knowledge you have.
2. Design. If you assigned a treatment at random, nothing points into it (Section 1.6). If you selected units on some variable, that variable is conditioned on, whether you like it or not.
3. Mechanism and literature. Biology, physics, institutional rules and previous studies tell you which effects are possible. Each arrow in a published graph should be traceable to one of these.
4. Expert elicitation. When the literature is thin, structured interviews with people who know the system (clinicians, dispatchers, teachers) can fill gaps. Ask experts about missing arrows as well as present ones: "could $A$ affect $B$ by any route other than through $C$?"

A fifth source, data, is the subject of most of this lecture. It is the weakest of the five on its own, and the most useful as a check on the others.

2.2 A graph makes predictions you can test

A graph is not just a picture; it is a list of claims about the data. Every missing arrow implies that some pair of variables should be independent once the right other variables are held fixed. These implied conditional independencies follow from the blocking rules of Section 1.6, and they can be checked against data. If a predicted independence clearly fails, the graph is wrong somewhere, and the failure tells you where to look.

Let us check the textbook signalling network for the Sachs data. The sparsebn package distributes both the measurements and the consensus graph assembled from the biological literature. We use the 2,603 cells that received no targeted intervention.

library(tidyverse)
library(dagitty)
library(bnlearn)     # Bayesian-network structure learning
library(brms)

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

# Sachs et al. (2005) flow-cytometry data, as distributed in the sparsebn package.
# Download once from the CRAN mirror on GitHub, then load from disk.
get_sachs <- function(file) {
  path <- file.path("data", file)
  if (!file.exists(path)) {
    dir.create("data", showWarnings = FALSE)
    download.file(paste0("https://raw.githubusercontent.com/cran/sparsebn/master/data/", file),
                  path, mode = "wb", quiet = TRUE)
  }
  load(path, envir = globalenv())
}
get_sachs("cytometryContinuous.rda")
cells   <- as_tibble(cytometryContinuous$data)      # 7,466 cells x 11 molecules (log scale)
targets <- cytometryContinuous$ivn                  # which molecule, if any, was perturbed
table(map_chr(targets, ~ if (length(.x) == 0) "none" else paste(.x, collapse = "+")))

# the consensus network from the literature, as an edge list
molecules <- names(cells)
consensus <- map_dfr(molecules, function(child) {
  parents <- cytometryContinuous$dag[[child]]
  tibble(from = molecules[parents], to = rep(child, length(parents)))
})
nrow(consensus)
observational <- cells[lengths(targets) == 0, ]

 akt  mek none pip2  pka  pkc 
 911  799 2603  810  707 1636 
[1] 20

The consensus graph has 20 arrows among 11 molecules. Now ask dagitty what it predicts, and test each prediction. For continuous data, localTests computes the partial correlation that each implied independence says should be zero, with a 95% interval.

g_consensus <- dagitty(paste0("dag{", paste(consensus$from, "->", consensus$to,
                                            collapse = "; "), "}"))
length(impliedConditionalIndependencies(g_consensus))

tests <- localTests(g_consensus, data = as.data.frame(scale(observational)), type = "cis")
tests |> arrange(desc(abs(estimate))) |> head(5) |> round(3)
mean(abs(tests$estimate) < 0.05)     # share of predictions that hold to within 0.05

[1] 68
                          estimate p.value   2.5%  97.5%
jnk _||_ p38 | pka, pkc      0.280       0  0.245  0.316
jnk _||_ plc | pkc           0.109       0  0.071  0.147
mek _||_ pip3 | pip2, plc   -0.084       0 -0.122 -0.046
erk _||_ plc | pkc           0.083       0  0.044  0.121
akt _||_ plc | pip3, pkc     0.081       0  0.043  0.119
[1] 0.75

The graph makes 68 testable predictions. Most hold well: about three-quarters of the predicted partial correlations are within 0.05 of zero, which for biological data is remarkably good. But one prediction fails badly. The graph says that Jnk and P38 should be independent once PKA and PKC are held fixed; in the data their partial correlation is 0.28, with a tight interval. With 2,603 cells, that is not noise.

The failure points somewhere specific. Jnk and P38 are both "stress-activated" kinases, and in the cell both are activated by enzymes upstream of them that were not measured in this experiment. An unmeasured common cause, $U\to \text{Jnk}$ and $U\to \text{P38}$, would produce exactly this residual correlation. The test has not just said "the graph is wrong"; it has suggested the repair: add a latent common cause. That is the right attitude to these tests. They are diagnostics for the graph, as posterior predictive checks are diagnostics for the model.

2.3 Fingerprints of structure, and what they cannot tell apart

If graphs make testable predictions, can we run the logic backwards and learn the graph from the independencies in the data? Partly. Start with three variables. Chapter 1 met the four elementary structures; here are the three that can connect $X$ and $Z$ through $Y$ (the chain listed in both directions), with the independencies each one implies.

Three different causal stories (two chains and a fork) leave exactly the same fingerprint in the data. No amount of observational data on $X$, $Y$ and $Z$ can tell them apart. They form a Markov equivalence class. The collider, by contrast, has a fingerprint of its own: independence that disappears when you condition. That asymmetry is the only directional information that observational data give for free, and every discovery algorithm is built on it.

Explore every DAG on three variables in the widget. Set the three edges, and it shows the independencies your graph implies and every other DAG that implies the same ones.

Try the chain $X\to Y\to Z$: its class has three members. Now reverse one arrow to make the collider $X\to Y\leftarrow Z$: its class has one member, itself. Then join $X$ and $Z$ as well. With all three pairs joined, no independence is implied at all, and every acyclic orientation is equivalent: six DAGs that observational data cannot separate.

2.4 Learning a graph: the PC algorithm

The PC algorithm, named after its inventors Peter Spirtes and Clark Glymour, turns the fingerprints into a procedure. It assumes that every independence in the data comes from the graph's structure (an assumption called faithfulness, Section 2.6) and that there are no unmeasured common causes. It then works in three phases.

Step through the algorithm on a six-variable graph. To see the logic clearly, the widget uses a perfect "oracle" for each independence test, read off the true graph. Real tests on finite data make mistakes, which we will meet in Section 2.5.

Notice three things. First, the algorithm never asks more than a handful of questions per edge, because it only conditions on current neighbours; that is what makes it fast enough for hundreds of variables. Second, the single collider at $D$ orients not just its own two edges but, through propagation, every edge downstream of it. Third, the edge $A-B$ is left undirected. Both $A\to B$ and $A\leftarrow B$ produce exactly the same independencies, so the algorithm honestly reports that the data cannot decide.

The PC algorithm is constraint-based: it builds the graph from yes/no answers to independence tests. The other main family is score-based. These algorithms define a score for each candidate graph, such as the Bayesian information criterion or the marginal likelihood of the data under the graph, and search for the highest-scoring graph, typically by adding, removing and reversing one arrow at a time (hill climbing). The fully Bayesian version treats the graph itself as unknown, puts a prior on all graphs, and returns a posterior distribution over graphs, from which you can read the probability that each edge exists. That is the most honest output a discovery method can give, but the number of DAGs grows super-exponentially (there are about $4\times {10}^{18}$ DAGs on 10 nodes), so exact posteriors are possible only for small problems, and approximations use MCMC over graphs or, as we will do next, the bootstrap.

2.5 Discovery on real cells

Now the real test. Run PC on the 2,603 observational cells, and compare its output with the consensus graph. A short helper counts how many consensus edges were found, how many learned edges are not in the consensus, and how many were oriented each way.

edge_key <- function(a, b) paste(pmin(a, b), pmax(a, b))

score_graph <- function(g) {
  a <- as_tibble(arcs(g))
  undirected <- filter(a, paste(to, from) %in% paste(a$from, a$to))
  directed   <- anti_join(a, undirected, by = c("from", "to"))
  learned <- unique(edge_key(a$from, a$to))
  truth   <- edge_key(consensus$from, consensus$to)
  tibble(edges     = length(learned),
         found     = sum(truth %in% learned),          # of 20 consensus edges
         extra     = sum(!learned %in% truth),
         right_way = sum(paste(directed$from, directed$to) %in% paste(consensus$from, consensus$to)),
         reversed  = sum(paste(directed$to, directed$from) %in% paste(consensus$from, consensus$to)),
         undirected = nrow(undirected) / 2)
}

pc_obs <- pc.stable(as.data.frame(observational), alpha = 0.01)   # Gaussian tests
hc_obs <- hc(as.data.frame(observational))                          # hill climbing, BIC score
bind_rows(PC = score_graph(pc_obs), `hill climbing` = score_graph(hc_obs), .id = "method")

# A tibble: 2 x 7
  method        edges found extra right_way reversed undirected
  <chr>         <int> <int> <int>     <int>    <int>      <dbl>
1 PC               11     6     5         2        1          3
2 hill climbing    15     7     8         5        2          0

The results are sobering. From 2,603 cells, PC recovers only 6 of the 20 consensus edges, adds 5 that are not in the consensus, and orients just 2 correctly. Hill climbing finds a few more edges, 7 of 20, but also 8 extras. These are not bad algorithms. They are honest about how little observational data can say: signalling relationships are noisy and nonlinear, the Gaussian tests are approximations, and many of the true edges belong to equivalence classes whose directions no observational method could recover.

Sachs and colleagues had something better than more cells. They had experiments. In some batches a chemical blocked or stimulated a particular molecule, which in graph terms means its incoming arrows were cut (Section 2.6). A molecule whose level was set by the experimenter cannot be caused by its usual parents in those cells, and that breaks the symmetries that make equivalence classes. The bnlearn package can use this information with the modified Bayesian Dirichlet equivalent (mbde) score, which ignores, for each molecule, the cells in which it was manipulated. Following Sachs et al., this score works with the data discretised into three levels (low, medium, high).

get_sachs("cytometryDiscrete.rda")        # same cells, discretized to low / medium / high
cells_d <- as.data.frame(lapply(cytometryDiscrete$data, factor, levels = 0:2))

# for each molecule, the rows (cells) in which it was experimentally perturbed
perturbed <- map(set_names(names(cells_d)), function(m)
  which(map_lgl(cytometryDiscrete$ivn, ~ m %in% .x)))
perturbed <- perturbed[lengths(perturbed) > 0]
lengths(perturbed)

hc_ignore <- hc(cells_d, score = "bde", iss = 1)                         # pretend no experiments
hc_int    <- hc(cells_d, score = "mbde", exp = perturbed, iss = 1)      # use the experiments
bind_rows(`ignoring interventions` = score_graph(hc_ignore),
          `using interventions`    = score_graph(hc_int), .id = "method")

 mek pip2  akt  pka  pkc 
 600  600  600  600 1200 
# A tibble: 2 x 7
  method                 edges found extra right_way reversed undirected
  <chr>                  <int> <int> <int>     <int>    <int>      <dbl>
1 ignoring interventions    24    15     9         7        8          0
2 using interventions       25    17     8        14        3          0

The difference is dramatic. Ignoring which cells were perturbed, hill climbing on the same 5,400 cells finds 15 consensus edges but orients only 7 of them correctly and 8 backwards. Using the experiments, it finds 17 of the 20 consensus edges and orients 14 correctly. Same algorithm, same cells, one extra piece of information: which molecules were set by hand. Compare the learned graphs with the consensus in the widget.

Look also at what the methods get wrong, because one mistake is shared by all of them. Every method in the widget, with or without the experiments, adds an edge between Jnk and P38 that the consensus graph does not have. That is exactly the pair whose predicted independence failed in Section 2.2. A learning algorithm that assumes no hidden variables has only one way to explain a residual correlation, which is to draw an arrow. Two different lines of evidence point to the same gap in the textbook graph, and the most plausible repair is the one suggested earlier: an unmeasured upstream activator of both stress kinases.

2.6 The limits of discovery, and how to push them

Three assumptions sit under every result in this lecture, and a careful reader of a discovery paper checks each one.

Faithfulness. PC assumes that every independence in the data comes from the graph's structure, never from effects that happen to cancel out. Suppose exercise lowers blood pressure directly but also raises appetite, which raises blood pressure by exactly the same amount. Exercise and blood pressure would then be independent even though exercise causes blood pressure along two paths. Exact cancellation is unlikely (for parameters drawn at random it has probability zero), but near-cancellation in finite samples is not, and systems that are regulated to stay stable, like body temperature or blood sugar, cancel by design.

No unmeasured common causes. PC and hill climbing assume every common cause of two measured variables is itself measured. Section 2.2 showed this failing for Jnk and P38. The FCI algorithm (Fast Causal Inference) relaxes the assumption: it returns a graph with marks meaning "this may be due to a hidden common cause", at the price of orienting even fewer edges.

Observational data identify only the equivalence class. Three escapes exist, each buying more structure with more assumptions:

• Interventions set a variable from outside, deleting the arrows into it. Section 2.5 showed what they can do.
• Time order forbids arrows from later variables to earlier ones; encode it as a blacklist.
• Functional assumptions. If effects are linear and the noise is non-Gaussian, the direction of every edge is identifiable from observational data (the LiNGAM method of Shimizu and colleagues, 2006). Related "additive noise" methods use nonlinearity. Intuitively, if $Y=f(X)+\epsilon$ with $\epsilon$ independent of $X$, then usually no model $X=g(Y)+\eta$ with $\eta$ independent of $Y$ exists, so the residuals reveal the direction.

Lecture 4 — Where priors come from

2.7 When the prior decides the answer

In the vaccine analysis of Chapter 1, 170 cases overwhelmed every reasonable prior. That is the usual situation for a headline result. It is not the usual situation for everything else in a paper: subgroups, rare outcomes, small pilot studies, parameters that the data barely inform. Here the choice of prior can decide the answer.

The FDA's briefing document on the same trial (FDA, 2020) reported efficacy within subgroups. Some subgroups had very few cases. Among participants aged 75 and over there were 0 cases in the vaccine arm and 5 in the placebo arm. Among Asian participants, 1 and 4. Among Black or African American participants, 0 and 7. The frequentist 95% confidence interval for VE in the oldest group ran from −12.1% to 100%: the data were consistent with a vaccine that was useless in the elderly.

We can analyse each subgroup exactly as in Section 1.2, under five priors on $\theta$ (the share of cases in the vaccine arm):

• flat, Beta(1, 1): every value of $\theta$ equally plausible;
• Jeffreys, Beta(0.5, 0.5): the "reference" prior of Section 2.9;
• protocol, Beta(0.7, 1): the trial's own prior;
• sceptical, Beta(5, 5): centred on no effect and worth 10 cases;
• rest of trial: centred on the efficacy seen in everyone outside the subgroup, but deliberately down-weighted to be worth only 10 cases.

ve <- function(theta) (1 - 2 * theta) / (1 - theta)
subgroups <- tribble(
  ~group,                       ~vax, ~pbo, ~fda_lo, ~fda_hi,
  "Age 75 and over",               0,    5,  -0.121,   1.000,
  "Asian",                         1,    4,  -1.587,   0.995,
  "Black or African American",     0,    7,   0.304,   1.000,
  "Age 65 to 74",                  1,   14,   0.532,   0.998)
all_vax <- 9; all_pbo <- 169          # whole trial, same analysis population

priors_for <- function(v, p) {
  rest <- (all_vax - v) / (all_vax + all_pbo - v - p)   # theta outside the subgroup
  tribble(~prior,                       ~a,               ~b,
          "flat Beta(1, 1)",              1,                1,
          "Jeffreys Beta(0.5, 0.5)",      0.5,              0.5,
          "protocol Beta(0.7, 1)",        0.700102,         1,
          "sceptical Beta(5, 5)",         5,                5,
          "rest of trial (10 cases)",     10 * rest,        10 * (1 - rest))
}

post <- subgroups |>
  mutate(priors = map2(vax, pbo, priors_for)) |>
  unnest(priors) |>
  mutate(median = ve(qbeta(0.5,   a + vax, b + pbo)),
         lo     = ve(qbeta(0.975, a + vax, b + pbo)),     # VE falls as theta rises
         hi     = ve(qbeta(0.025, a + vax, b + pbo)),
         p_30   = pbeta(0.7 / 1.7, a + vax, b + pbo))     # Pr(VE > 30%)

post |> filter(group == "Age 75 and over") |>
  select(prior, median, lo, hi, p_30) |> mutate(across(where(is.numeric), \(x) round(x, 2)))

ggplot(post, aes(y = fct_rev(fct_inorder(prior)))) +
  geom_linerange(aes(xmin = pmax(lo, -1), xmax = hi), colour = harbour[["sky"]], linewidth = 1.2) +
  geom_point(aes(x = median), colour = harbour[["sky"]], size = 2) +
  geom_linerange(data = subgroups, aes(y = "FDA 95% confidence interval",
                 xmin = pmax(fda_lo, -1), xmax = fda_hi), colour = harbour[["muted"]], linewidth = 1.2) +
  geom_vline(xintercept = 0.3, linetype = "dotted", colour = harbour[["rose"]]) +
  facet_wrap(~ fct_inorder(group), ncol = 2) +
  scale_x_continuous(labels = scales::percent, limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  labs(x = "vaccine efficacy (intervals cut at -100%)", y = NULL) +
  theme(panel.spacing.x = unit(1.6, "lines"))

# A tibble: 5 x 5
  prior                    median    lo    hi  p_30
  <chr>                     <dbl> <dbl> <dbl> <dbl>
1 flat Beta(1, 1)            0.88  0.15  1     0.96
2 Jeffreys Beta(0.5, 0.5)    0.96  0.39  1     0.98
3 protocol Beta(0.7, 1)      0.93  0.33  1     0.98
4 sceptical Beta(5, 5)       0.52 -0.39  0.85  0.75
5 rest of trial (10 cases)   0.98  0.8   1     1   

For the over-75s, the posterior median VE ranges from 52% (sceptical prior) to 98% (rest-of-trial prior), and Pr(VE > 30%) from 0.75 to essentially 1. Every one of these analyses uses the same five cases. The prior is not a technicality here; it is the analysis. Which prior is right? None is right in the abstract. Each answers a different question:

• The flat, Jeffreys and protocol priors ask what these few cases say on their own. Their answers are similar to each other and to the confidence interval: probably effective, but very uncertain.
• The sceptical prior asks whether the data are strong enough to convince someone who doubted the vaccine worked at all. For the oldest group, they are not.
• The rest-of-trial prior asks what we should believe if the vaccine works similarly in all subgroups unless the data say otherwise. Biologically that is the most defensible starting point, and it is the logic of the multilevel model of Chapter 6, which learns from the data how much the subgroups should share rather than fixing it at 10 cases.

2.8 Maximum entropy: assuming as little as possible

Suppose you know something about a quantity, such as its possible values and its average, but nothing more. Which distribution expresses exactly that knowledge, and no more? The answer, developed by the physicist Edwin Jaynes from Claude Shannon's theory of information, is: the distribution with the largest entropy among all those that satisfy what you know.

Why maximise it? Here is the intuition, due to Graham Wallis and retold by Jaynes and McElreath. Imagine distributing $N$ pebbles at random among $K$ buckets, each pebble equally likely to land in any bucket. Some arrangements of counts can happen in many more ways than others. The number of ways to obtain counts $(n_1,\dots ,n_K)$ is the multinomial coefficient $W=N!/(n_1!\cdots n_K!)$, and for large $N$,

So the distribution with the highest entropy is the one that can happen in the most ways. Choosing it means choosing the outcome that is most likely if nothing beyond your constraints is pushing the pebbles around. Any other choice smuggles in extra structure that you have not justified.

Jaynes's dice

Jaynes illustrated the principle with a puzzle. A six-sided die is rolled many times and you are told only that the average score was 4.5, not the fair 3.5. What probability should you give each face? Infinitely many distributions have mean 4.5. The maximum-entropy one solves

faces <- 1:6
mean_given_lambda <- function(lambda) {
  p <- exp(lambda * faces); p <- p / sum(p)
  sum(faces * p)
}
lambda <- uniroot(\(l) mean_given_lambda(l) - 4.5, c(-5, 5))$root
p_maxent <- exp(lambda * faces) / sum(exp(lambda * faces))
round(p_maxent, 4)

entropy <- function(p) -sum(ifelse(p > 0, p * log(p), 0))
c(maxent = entropy(p_maxent),
  `only 4s and 5s` = entropy(c(0, 0, 0, 0.5, 0.5, 0)),
  `3s and 6s` = entropy(c(0, 0, 0.5, 0, 0, 0.5)),
  fair_die = entropy(rep(1/6, 6)))

[1] 0.0544 0.0788 0.1142 0.1654 0.2398 0.3475
        maxent only 4s and 5s      3s and 6s       fair_die 
        1.6136         0.6931         0.6931         1.7918 

The maximum-entropy answer rises smoothly from 5.4% for a one to 34.7% for a six: tilted towards high faces just enough to make the mean 4.5, and otherwise as even as possible. The alternatives with the same mean, such as "only fours and fives", have far lower entropy, because they assert things (no sixes ever!) that the information given does not support. Move the mean in the widget.

Almost every distribution you will use in this course is the maximum-entropy distribution for some simple state of knowledge.

This changes how to think about the normal distribution. It is not used because nature is normal (often it isn't). It is used because if all you are prepared to say about a quantity is its centre and spread, the normal is the distribution that says nothing more. A quick numerical check makes the point: among several familiar distributions with the same variance, the normal has the largest entropy.

# differential entropies of distributions with variance 1
c(normal   = 0.5 * log(2 * pi * exp(1)),
  logistic = 2 + log(sqrt(3) / pi),          # scale s = sqrt(3)/pi
  Laplace  = 1 + log(2 / sqrt(2)),           # scale b = 1/sqrt(2)
  uniform  = log(sqrt(12)))                  # width sqrt(12)

  normal logistic  Laplace  uniform 
   1.419    1.405    1.347    1.242 

The same logic applies to likelihoods, which is why Chapter 4 will choose binomial, Poisson and gamma likelihoods by asking what we know about the outcome rather than by habit. And it applies to priors: if all you will commit to is that a standard deviation is positive with some typical size, the exponential distribution is the honest prior, and that is why $\sigma \sim \text{Exponential}(1)$ appears throughout these notes.

2.9 Every prior lives on a scale

Maximum entropy has a catch, and it is the same catch that made the vaccine prior look strange in Section 1.2. A prior that is "flat" or "uninformative" on one scale is not flat on another. A flat prior on $\theta$, the share of cases in the vaccine arm, is anything but flat on VE.

set.seed(5)
theta <- runif(1e5)
tibble(`theta (share of cases in vaccine arm)` = theta,
       `VE = (1 - 2 theta) / (1 - theta)`   = ve(theta)) |>
  pivot_longer(everything()) |>
  mutate(name = fct_inorder(name)) |>
  filter(value > -3) |>
  ggplot(aes(value)) +
  geom_histogram(bins = 60, fill = harbour[["sky"]], colour = NA) +
  facet_wrap(~ name, scales = "free") +
  labs(x = NULL, y = "draws")
mean(ve(theta) > 0.3)          # prior probability that VE exceeds 30%

[1] 0.4171

A "flat" prior on $\theta$ puts 42% of its weight on VE above 30%, a density that rises steadily towards perfect efficacy, and fully half its weight on a vaccine that does harm (VE below 0), spread along a tail that reaches minus infinity. Had we instead put a flat prior on VE between −1 and 1, the implied prior on $\theta$ would have been far from flat. There is no scale-free way to be ignorant.

Two responses to this problem are worth knowing. The first, due to Harold Jeffreys (1946), is to choose priors that give the same answer whichever parameterisation you use. For the binomial, the Jeffreys prior is Beta(0.5, 0.5), the prior used in Section 2.7. It has an attractive mathematical property (invariance under reparameterisation), but it is often improper or odd-looking in more complex models, and it is rarely used for regression coefficients.

The second response, which is the one this course takes, is practical: choose the scale on which you have knowledge, and put the prior there. You probably know more about efficacy than about $\theta$, and more about a slope in kilograms per millimetre than about a slope in grams per micrometre. Then check, by prior predictive simulation, that the implied prior on everything else you care about is reasonable. Maximum entropy tells you what shape to use on a given scale; the scale itself is a modelling decision you must justify.

2.10 Weakly informative priors in practice

Between "flat" and "informative" lies the default of modern applied work: weakly informative priors. They rule out the absurd, gently pull estimates towards plausible values, and let the data decide everything else. Three habits make them easy to choose.

1. Put variables on interpretable scales. Centre predictors at meaningful values (Section 1.7 centred flipper length at 200 mm). Or standardise outcome and predictors to mean 0 and standard deviation 1, after which a slope of 1 means a one-standard-deviation change in the outcome per standard deviation of the predictor, which is an enormous effect in most fields. Then Normal(0, 1), or Normal(0, 0.5), is a sensible weakly informative prior for every slope.
2. Use maximum-entropy shapes on that scale: normal for locations and slopes, exponential or half-normal for standard deviations, and (Chapter 6) the LKJ distribution for correlation matrices.
3. Check by prior predictive simulation, as in Section 1.7, but now let brms do it.

First, see what priors brms would use if you said nothing. It chooses its defaults from the data, which is convenient but means they are not really priors in the strict sense.

library(palmerpenguins)
pen <- penguins |>
  drop_na(body_mass_g, flipper_length_mm, sex) |>
  mutate(mass = body_mass_g / 1000, flipper = flipper_length_mm - 200)

get_prior(mass ~ 1 + flipper + species + sex, data = pen)[, c("prior", "class", "coef", "source")]

                prior     class             coef  source
 student_t(3, 4, 2.5) Intercept                  default
               (flat)         b                  default
               (flat)         b          flipper default
               (flat)         b          sexmale default
               (flat)         b speciesChinstrap default
               (flat)         b    speciesGentoo default
 student_t(3, 0, 2.5)     sigma                  default

The defaults put flat, improper priors on the slopes, marked (flat), and Student-$t$ priors centred and scaled using the data for the intercept and $\sigma$. Flat slopes are harmless with 333 penguins, but they are not what you want in a paper, and they make prior predictive simulation impossible. So state your own priors and simulate from them alone, with sample_prior = "only", which runs the sampler while ignoring the likelihood.

m_prior <- brm(
  mass ~ 1 + flipper + species + sex, data = pen,
  prior = c(prior(normal(4, 1), class = Intercept),
            prior(normal(0, 1), class = b),
            prior(normal(0, 0.05), class = b, coef = flipper),
            prior(exponential(1), class = sigma)),
  sample_prior = "only",                       # ignore the data: simulate from the prior
  seed = 1, refresh = 0, file = "fits/ch2_prior_only"
)
# pp_check(m_prior) draws this plot in one line; here we build it by hand
set.seed(2)
sims <- posterior_predict(m_prior, ndraws = 60)          # 60 simulated datasets
sim_long <- tibble(draw = rep(1:60, times = ncol(sims)), mass = as.vector(sims))
ggplot() +
  stat_density(data = sim_long, aes(mass, group = draw), geom = "line", position = "identity",
               n = 128, colour = harbour[["sky"]], alpha = 0.35) +
  stat_density(data = pen, aes(mass), geom = "line", n = 128,
               colour = harbour[["ink"]], linewidth = 1.2) +
  coord_cartesian(xlim = c(-4, 12)) +
  labs(x = "body mass (kg)", y = "density")
mean(sims < 0)                                             # share of impossible penguins

[1] 0.02533

The thick line is the observed distribution of masses; the thin lines are datasets the priors consider plausible. They are much wider than the data, as weakly informative priors should be, but they produce no elephants, and only about 2.5% of simulated penguins (the number printed below the plot) weigh less than nothing: an acceptable price for priors this weak. These are the priors we used for model m2 in Chapter 1, and this plot is the evidence that they were reasonable.

2.11 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Depaoli, S. and van de Schoot, R. (2017). Improving transparency and replication in Bayesian statistics: the WAMBS-checklist. Psychological Methods, 22(2), 240–261.
2. FDA (2020). Vaccines and Related Biological Products Advisory Committee Briefing Document: Pfizer-BioNTech COVID-19 Vaccine, 10 December 2020.
3. Gelman, A., Simpson, D. and Betancourt, M. (2017). The prior can often only be understood in the context of the likelihood. Entropy, 19(10), 555.
4. Glymour, C., Zhang, K. and Spirtes, P. (2019). Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10, 524.
5. Jaynes, E. T. (1957). Information theory and statistical mechanics. Physical Review, 106(4), 620–630.
6. Jaynes, E. T. (2003). Probability Theory: The Logic of Science. Cambridge University Press.
7. Jeffreys, H. (1946). An invariant form for the prior probability in estimation problems. Proceedings of the Royal Society A, 186, 453–461.
8. Kruschke, J. K. (2021). Bayesian analysis reporting guidelines. Nature Human Behaviour, 5, 1282–1291.
9. Meek, C. (1995). Causal inference and causal explanation with background knowledge. Proceedings of UAI, 403–410.
10. Sachs, K., Perez, O., Pe'er, D., Lauffenburger, D. A. and Nolan, G. P. (2005). Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721), 523–529.
11. Scutari, M. (2010). Learning Bayesian networks with the bnlearn R package. Journal of Statistical Software, 35(3), 1–22.
12. Shannon, C. E. (1948). A mathematical theory of communication. Bell System Technical Journal, 27, 379–423.
13. Shimizu, S., Hoyer, P. O., Hyvärinen, A. and Kerminen, A. (2006). A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7, 2003–2030.
14. Spirtes, P., Glymour, C. and Scheines, R. (2000). Causation, Prediction, and Search (2nd ed.). MIT Press.
15. Verma, T. and Pearl, J. (1990). Equivalence and synthesis of causal models. Proceedings of UAI, 255–270.