Selection and the replication crisis

Class

Author

Joshua Loftus

To start

Two things, on the shared document.

  1. One ethical issue in data science that concerns you. A sentence is enough.
  2. From memory: what does the power of a test depend on? And which way does selecting on significance push a published estimate?

Problems

Everything can be done by hand with \(\Phi(1) = 0.84\), \(\Phi(2) = 0.977\), \(\phi(0) = 0.399\), \(\phi(1) = 0.242\) and \(\phi(2) = 0.054\). Use \(c = 2\) as the threshold throughout; it is close enough to 1.96.

Parts marked (class) are for the class. The rest are for practice afterwards. Answers are in the collapsed boxes. Attempt each part first.

1. Many tests

An analyst tests \(k\) true null hypotheses, independently, each at level \(0.05\), and reports the smallest \(p\)-value.

  1. (class) What is the probability that the reported \(p\)-value is below 0.05, for \(k = 1\), \(k = 10\) and \(k = 20\)?

  2. A second analyst announces in advance that they will test only hypothesis number 7 of the same 20, and reports that. What is the probability their reported \(p\)-value is below 0.05? What has changed, and what has not?

  1. \(1 - 0.95^k\): \(0.05\) at \(k=1\); \(1 - 0.95^{10} = 0.40\); \(1 - 0.95^{20} = 0.64\).

  2. \(0.05\). Same tests, same data. The reported comparison was fixed before the data could influence the choice, so the advertised error rate is now the error rate actually being run. Pre-specification does not change the power of anything. It changes which comparison you are allowed to call “the” result.

2. Power

An estimator \(\hat\mu\) of \(\mu\) has standard error \(\sigma/\sqrt n\), and \(Z = \hat\mu/\mathrm{se} \sim N(\theta, 1)\) with \(\theta = \mu\sqrt n/\sigma\). A one-sided test rejects when \(Z > 2\).

  1. (class) Compute the power at \(\theta = 1\), \(\theta = 2\) and \(\theta = 3\).

  2. The study at \(\theta = 1\) is repeated with nine times the sample size. What is the new \(\theta\), and the new power?

  3. Which of \(\mu\), \(\sigma\) and \(n\) could the analyst have changed, and which does the power describe?

  1. Power \(= P(Z > 2) = 1 - \Phi(2 - \theta)\): at \(\theta=1\), \(1-\Phi(1) = 0.16\); at \(\theta=2\), \(1-\Phi(0) = 0.50\); at \(\theta=3\), \(1-\Phi(-1) = \Phi(1) = 0.84\).

  2. \(\theta\) scales with \(\sqrt n\), so nine times the sample gives \(\theta = 3\) and power \(0.84\).

  3. The analyst controls \(n\), and sometimes \(\sigma\) through better measurement. \(\mu\) is what it is. Power describes how much the estimate varies relative to the effect: at \(\theta = 1\) the standard error is the same size as the effect, so a significant result needs a lucky draw.

3. The filter

Same setting. Only results with \(Z > 2\) are reported. Recall \[E[Z \mid Z > c] = \theta + \frac{\phi(c-\theta)}{1-\Phi(c-\theta)}.\]

  1. (class) Compute \(E[Z \mid Z > 2]\) when \(\theta = 0\). In words: what does this number say about a true null that got reported?

  2. (class) Compute \(E[Z \mid Z > 2]\) and the ratio \(E[Z \mid Z > 2]/\theta\) for \(\theta = 1\) and for \(\theta = 3\). Which of these two studies reports a number closer to the truth, and what quantity from question 2 explains the difference?

  3. Is the problem in (b) a bias, a variance, or both? Say which word applies to which quantity.

  1. \(\phi(2)/(1-\Phi(2)) = 0.054/0.023 \approx 2.37\). A reported true null is, on average, about 2.4 standard errors from zero, and the filter created that from nothing.

  2. \(\theta = 1\): \(1 + \phi(1)/(1-\Phi(1)) = 1 + 0.242/0.16 \approx 2.5\); ratio about 2.5. \(\theta = 3\): \(3 + \phi(-1)/(1-\Phi(-1)) = 3 + 0.242/0.84 \approx 3.29\); ratio about 1.1. The high-powered study exaggerates by a tenth, the low-powered one by a factor of two and a half. The power from question 2 (0.16 versus 0.84) is the denominator of the bias term.

  3. Both, and they are related. The reported estimate is biased upward. The size of that bias is set by the estimate’s variance relative to the effect, which is what the power measures. The filter converts variance into bias. An unfiltered estimate at \(\theta = 1\) would be noisy but unbiased.

4. What a field’s literature looks like (harder)

In a field, a fraction \(p\) of the hypotheses tested are true. Tests are run at \(\alpha = 0.05\) with power \(1 - \beta\), and every rejection is published.

  1. Write down, using Bayes’ rule, the probability that a published rejection is a false positive, in terms of \(p\), \(\alpha\) and \(\beta\).

  2. Evaluate it for \(p = 0.1\) with power \(0.8\), and for \(p = 0.1\) with power \(0.2\).

  3. A critic says: “this proves most published findings are false.” A defender says: “this proves nothing about any real field.” In two sentences, say what the formula shows and what it does not.

  4. Which of the three ingredients is about who is in the data?

  1. \(P(H_0 \mid \text{reject}) = \dfrac{(1-p)\alpha}{(1-p)\alpha + p(1-\beta)}\).

  2. \(p=0.1\), power \(0.8\): \(0.045/(0.045 + 0.08) = 0.36\). Power \(0.2\): \(0.045/(0.045 + 0.02) = 0.69\).

  3. It shows that if true effects are rare and studies are small, most significant findings are false, with no misconduct required, and that low power makes it much worse. It does not show anything about a particular field, because \(p\) is not observed and the numbers depend on it. The arithmetic is not in dispute; the inputs are.

  4. \(p\). The hypotheses a field chooses to test are a sample, and a field that tests wild ideas has a different literature from one that tests safe ones, even at the same \(\alpha\) and power. Who is in the data, and who decided.

Discussion break

In groups of three or four, then to the room:

What are the causes of the replication crisis? Consider both cultural and statistical sources. Try to name at least one of each, and say which of your causes the problem set has anything to say about.

Coding

Open a fresh R script. You will need ggplot2 and dplyr.

A. Significance filter / file-drawer effect

Generate 100,000 standard normal \(z\)-scores. Every one of them is a test of a true null hypothesis.

z <- rnorm(1e5)
  1. Keep only the \(z\)-scores with \(|z| > 1.96\). What fraction survived? What should it be?
  2. Plot a histogram of the survivors. Describe its shape in one sentence.
  3. Compute the mean of \(|z|\) among the survivors. Compare it with \(\phi(1.96)/(1 - \Phi(1.96))\), which in R is dnorm(1.96) / (1 - pnorm(1.96)).
  4. (Practice.) Now generate \(z \sim N(1, 1)\) instead, so every hypothesis is false with a modest effect. Keep the survivors with \(z > 1.96\) and compute their mean. Compare with the formula \(\theta + \phi(c-\theta)/(1-\Phi(c-\theta))\) at \(\theta = 1\), \(c = 1.96\). What is the exaggeration ratio?
survivors <- z[abs(z) > 1.96]
length(survivors) / length(z)
[1] 0.04955
data.frame(z = survivors) |>
  ggplot(aes(z)) + geom_histogram(bins = 60)

mean(abs(survivors))
[1] 2.331199
dnorm(1.96) / (1 - pnorm(1.96))
[1] 2.337835
z1 <- rnorm(1e5, mean = 1)
surv1 <- z1[z1 > 1.96]
mean(surv1)
[1] 2.491784
1 + dnorm(1.96 - 1) / (1 - pnorm(1.96 - 1))
[1] 2.493194
mean(surv1) / 1
[1] 2.491784

About 5% survive. The histogram has two humps and nothing in the middle. Every survivor looks like strong evidence, and none of them is. The mean of \(|z|\) among survivors is about 2.34. At \(\theta = 1\) the survivors average about 2.49, an exaggeration ratio of about 2.5.

B. The best of several

A function that generates a dataset from a linear model, and a function that runs one experiment and returns the \(p\)-value for the coefficient on \(X\).

generate_data <- function(beta0 = 0, betaX = 0, betaA = 0, betaXA = 0,
                          n = 200, proportion = 1/2) {
  X <- rnorm(n)
  A <- rbinom(n, 1, proportion)
  Y <- beta0 + betaX * X + betaA * A + betaXA * X * A + rnorm(n, sd = 1)
  data.frame(Y = Y, X = X, A = factor(A))
}

one_experiment <- function(beta0 = 0, betaX = 0, betaA = 0, betaXA = 0,
                           n = 200, proportion = 1/2) {
  one_sample <- generate_data(beta0, betaX, betaA, betaXA, n, proportion)
  one_fit <- lm(Y ~ X + A + X * A, one_sample)
  summary(one_fit)$coefficients[2, 4]
}
  1. With all coefficients at zero, run 1000 experiments and compute the rejection rate at level 0.05. Is it what it should be?
  2. Write one_searched_experiment(), which fits the same model, looks at the \(p\)-values for all three non-intercept coefficients, also fits the simpler model Y ~ X, and returns the smallest of the four \(p\)-values. Run 1000 of these under the null. What is the rejection rate now, and why is it not \(1 - 0.95^4\)?
  3. Set betaX = 0.2 and find, by trying values of n, the sample size at which the rejection rate is about 0.8. Then explain the number using \(\theta = \mu\sqrt{n_{\text{eff}}}/\sigma\): what is \(n_{\text{eff}}\) here, and why is it not \(n\)?
  4. (Practice.) Run 1000 experiments with betaX = 0.2 at the \(n\) you found in 3, keep the estimated coefficient on \(X\) only when its \(p\)-value is below 0.05, and compute the average kept estimate. Compare it with 0.2. You will need to modify one_experiment() to return the estimate as well as the \(p\)-value.
rejection_level <- 0.05
n_experiments <- 1000
mean(replicate(n_experiments, one_experiment()) < rejection_level)
[1] 0.049
one_searched_experiment <- function(beta0 = 0, betaX = 0, betaA = 0,
                                    betaXA = 0, n = 200, proportion = 1/2) {
  one_sample <- generate_data(beta0, betaX, betaA, betaXA, n, proportion)
  full_fit <- lm(Y ~ X + A + X * A, one_sample)
  simple_fit <- lm(Y ~ X, one_sample)
  min(summary(full_fit)$coefficients[-1, 4],
      summary(simple_fit)$coefficients[2, 4])
}
mean(replicate(n_experiments, one_searched_experiment()) < rejection_level)
[1] 0.162
sapply(c(200, 300, 400, 500), function(n)
  mean(replicate(n_experiments, one_experiment(betaX = 0.2, n = n)) < rejection_level))
[1] 0.509 0.680 0.824 0.877
one_experiment_estimate <- function(betaX = 0.2, n = 400) {
  one_sample <- generate_data(betaX = betaX, n = n)
  one_fit <- lm(Y ~ X + A + X * A, one_sample)
  summary(one_fit)$coefficients[2, c(1, 4)]
}
results <- t(replicate(n_experiments, one_experiment_estimate()))
colnames(results) <- c("estimate", "p_value")
results <- as.data.frame(results)
mean(results$estimate)                              # all experiments
[1] 0.2002631
mean(results$estimate[results$p_value < 0.05])      # the reported ones
[1] 0.2233063
  1. Close to 0.05.
  2. About 0.16 rather than \(1 - 0.95^4 = 0.19\) (0.160 in 20,000 runs). The four \(p\)-values are not independent: the \(X\) coefficient in the two models and the interaction share information, so the search buys less than four independent tries would. Still more than three times the advertised rate, and in the same direction as the significance filter.
  3. About \(n = 400\). The coefficient on \(X\) in Y ~ X + A + X * A is the slope within the group \(A = 0\), which holds roughly half the sample, so \(n_{\text{eff}} \approx n/2\) and \(\theta \approx 0.2\sqrt{200} \approx 2.8\), where the power table gives 0.80.
  4. The kept estimates average noticeably above 0.2, even at 80% power, because only the draws that came out large enough survive. Reduce \(n\) and the gap widens.

Discussion

More prompts than we will get through. Groups of three or four, then to the room; every group’s answer goes on the shared document.

  1. What are the ethical impacts of the replication crisis? Consider both general impacts and impacts specific to a field you know.
  2. What can data scientists and statisticians do about it? Which of your proposals would have changed the number in problem 3(a), which would have changed the number in problem 4(b), and which would change neither?
  3. Each group takes one setting below and writes three sentences: what the filter is; which way the numbers you get to see are biased and what makes the bias larger or smaller; who benefits from the bias being invisible.
    • An investment platform advertises the past returns of the funds it currently offers. Funds that performed badly were closed and are no longer listed.
    • A newspaper publishes a table of the schools whose exam results improved most this year.
    • A hospital’s website shows the complication rates of its surgeons. Only surgeons who have performed at least 50 operations are shown.
    • A recruiting team reports that candidates hired after a new interview process are performing well. It cannot see the candidates it rejected.
    • A social media feed shows you the ten posts from a group of people that got the most reactions today.
  4. Have you seen “spin” in a paper, a press release or a news story: a result reported in the most favorable of the comparisons that could have been made? What was the comparison that was not reported?

To finish

One sentence, individually: what causes p-hacking, and what would fix it?

For next time

Bring one published claim about a difference between two groups of people, from any source. Write down the comparison that was reported, and one comparison that could have been reported instead. One sentence each. Every group’s examples will go on the screen.