12
$\begingroup$

Let's stick to an ideal situation with random sampling, Gaussian populations, equal variances, no P-hacking, etc.

Step 1. You run an experiment say comparing two sample means, and compute a 95% confidence interval for the difference between the two population means.

Step 2. You run many more experiments (thousands). The difference between means will vary from experiment to experiment due to random sampling.

Question: What fraction of the difference between means from the collection of experiments in step 2 will lie within the confidence interval of step 1?

That can't be answered. It all depends on what happened in step 1. If that step 1 experiment was very atypical, the answer to the question might be very low.

So imagine that both steps are repeated many times (with step 2 repeated many more times). Now it ought to be possible, I would think, to come up with an expectation for what fraction of repeat experiments, on average, have an effect size within the 95% confidence interval of the first experiment.

It seems that the answer to these questions need to be understood to evaluate reproducibility of studies, a very hot area now.

$\endgroup$
4
  • $\begingroup$ For each original (step 1) experiment $i$, define $x_i$ as the fraction of subsequent (step 2) results that produce findings within the original result's confidence interval. You want to compute the empirical distribution of $x$? $\endgroup$ – Matthew Gunn Jun 13 '16 at 8:55
  • $\begingroup$ Yes, you understand what I am asking $\endgroup$ – Harvey Motulsky Jun 13 '16 at 10:41
  • $\begingroup$ @MatthewGunn asked if you wanted the empirical distribution of the "capture fraction" for future observations. Your post asked "...it ought to be possible, I would think, to come up with an expectation for what fraction of repeat experiments, on average, have an effect size within the 95% confidence interval of the first experiment". This is not a distribution but an expected value (average). $\endgroup$ – user75138 Jun 13 '16 at 13:16
  • $\begingroup$ Whuber's analysis is great, but if you need a citation then here is a paper that discusses exactly this question in great detail: Cumming & Maillardet, 2006, Confidence Intervals and Replication: Where Will the Next Mean Fall?. They call it capture percentage of a confidence interval. $\endgroup$ – amoeba Dec 6 '16 at 22:40
12
$\begingroup$

Analysis

Because this is a conceptual question, for simplicity let's consider the situation in which a $1-\alpha$ confidence interval $$\left[\bar x^{(1)} + Z_{\alpha/2} s^{(1)}/\sqrt{n}, \bar x^{(1)} + Z_{1-\alpha/2} s^{(1)}/\sqrt{n}\right]$$ is constructed for a mean $\mu$ using a random sample $x^{(1)}$ of size $n$ and a second random sample $x^{(2)}$ is taken of size $m$, all from the same Normal$(\mu,\sigma^2)$ distribution. (If you like you may replace the $Z$s by values from the Student $t$ distribution of $n-1$ degrees of freedom; the following analysis will not change.)

The chance that the mean of the second sample lies within the CI determined by the first is

$$\Pr\left(\bar x^{(1)} + \frac{Z_{\alpha/2}}{\sqrt{n}} s^{(1)} \le \bar x^{(2)} \le \bar x^{(1)} + \frac{Z_{1-\alpha/2}}{\sqrt{n}} s^{(1)}\right) =\Pr\left(\frac{Z_{\alpha/2}}{\sqrt{n}} s^{(1)} \le \bar x^{(2)}-\bar x^{(1)} \le \frac{Z_{1-\alpha/2}}{\sqrt{n}} s^{(1)}\right).$$

Because the first sample mean $\bar x^{(1)}$ is independent of the first sample standard deviation $s^{(1)}$ (this requires normality) and the second sample is independent of the first, the difference in sample means $U = \bar x^{(2)} - \bar x^{(1)}$ is independent of $s^{(1)}$. Moreover, for this symmetric interval $Z_{\alpha/2}=-Z_{1-\alpha/2}$. Therefore, writing $S$ for the random variable $s^{(1)}$ and squaring both inequalities, the probability in question is the same as

$$\Pr\left(U^2 \le \left(\frac{Z_{1-\alpha/2}}{\sqrt{n}}\right)^2 S^2\right)= \Pr\left(\frac{U^2}{S^2} \le \left(\frac{Z_{1-\alpha/2}}{\sqrt{n}}\right)^2\right).$$

The laws of expectation imply $U$ has a mean of $0$ and a variance of

$$\operatorname{Var}(U) = \operatorname{Var}\left(\bar x^{(2)} - \bar x^{(1)}\right) = \sigma^2\left(\frac{1}{m} + \frac{1}{n}\right).$$

Since $U$ is a linear combination of Normal variables, it also has a Normal distribution. Therefore $U^2$ is $\sigma^2\left(\frac{1}{n} + \frac{1}{m}\right)$ times a $\chi^2(1)$ variable. We already knew that $S^2$ is $\sigma^2/n$ times a $\chi^2(n-1)$ variable. Consequently, $U^2/S^2$ is $1/n + 1/m$ times a variable with an $F(1,n-1)$ distribution. The required probability is given by the F distribution as

$$F_{1,n-1}\left(\frac{Z_{1-\alpha/2}^2}{1 + n/m}\right).\tag{1}$$


Discussion

An interesting case is when the second sample is the same size as the first, so that $n/m=1$ and only $n$ and $\alpha$ determine the probability. Here are the values of $(1)$ plotted against $\alpha$ for $n=2,5,20,50$.

Figure

The graphs rise to a limiting value at each $\alpha$ as $n$ increases. The traditional test size $\alpha=0.05$ is marked by a vertical gray line. For largish values of $n=m$, the limiting chance for $\alpha=0.05$ is around $85\%$.

By understanding this limit, we will peer past the details of small sample sizes and better understand the crux of the matter. As $n=m$ grows large, the $F$ distribution approaches a $\chi^2(1)$ distribution. In terms of the standard Normal distribution $\Phi$, the probability $(1)$ then approximates

$$\Phi\left(\frac{Z_{1-\alpha/2}}{\sqrt{2}}\right) - \Phi\left(\frac{Z_{\alpha/2}}{\sqrt{2}}\right) = 1 - 2\Phi\left(\frac{Z_{\alpha/2}}{\sqrt{2}}\right) .$$

For instance, with $\alpha=0.05$, $Z_{\alpha/2}/\sqrt{2} \approx -1.96/1.41 \approx -1.386$ and $\Phi(-1.386) \approx 0.083$. Consequently the limiting value attained by the curves at $\alpha=0.05$ as $n$ increases will be $1 - 2(0.083) = 1 - 0.166=0.834$. You can see it has almost been reached for $n=50$ (where the chance is $0.8383\ldots$.)

For small $\alpha$, the relationship between $\alpha$ and the complementary probability--the risk that the CI does not cover the second mean--is almost perfectly a power law. Another way to express this is that the log complementary probability is almost a linear function of $\log\alpha$. The limiting relationship is approximately

$$\log\left(2\Phi\left(\frac{Z_{\alpha/2}}{\sqrt{2}}\right)\right) \approx -1.79712 + 0.557203\log(20 \alpha) + 0.00657704 (\log(20 \alpha))^2 + \cdots$$

In other words, for large $n=m$ and $\alpha$ anywhere near the traditional value of $0.05$, $(1)$ will be close to

$$1 - 0.166 (20\alpha)^{0.557}.$$

(This reminds me very much of the analysis of overlapping confidence intervals I posted at https://stats.stackexchange.com/a/18259/919. Indeed, the magic power there, $1.91$, is very nearly the reciprocal of the magic power here, $0.557$. At this point you should be able to re-interpret that analysis in terms of reproducibility of experiments.)


Experimental results

These results are confirmed with a straightforwward simulation. The following R code returns the frequency of coverage, the chance as computed with $(1)$, and a Z-score to assess how much they differ. The Z-scores are typically less than $2$ in size, regardless of $n, m, \mu, \sigma, \alpha$ (or even whether a $Z$ or $t$ CI is computed), indicating the correctness of formula $(1)$.

n <- 3      # First sample size
m <- 2      # Second sample size
sigma <- 2 
mu <- -4
alpha <- 0.05
n.sim <- 1e4
#
# Compute the multiplier.
#
Z <- qnorm(alpha/2)
#Z <- qt(alpha/2, df=n-1) # Use this for a Student t C.I. instead.
#
# Draw the first sample and compute the CI as [l.1, u.1].
#
x.1 <- matrix(rnorm(n*n.sim, mu, sigma), nrow=n)
x.1.bar <- colMeans(x.1)
s.1 <- apply(x.1, 2, sd)
l.1 <- x.1.bar + Z * s.1 / sqrt(n)
u.1 <- x.1.bar - Z * s.1 / sqrt(n)
#
# Draw the second sample and compute the mean as x.2.
#
x.2 <- colMeans(matrix(rnorm(m*n.sim, mu, sigma), nrow=m))
#
# Compare the second sample means to the CIs.
#
covers <- l.1 <= x.2 & x.2 <= u.1
#
# Compute the theoretical chance and compare it to the simulated frequency.
#
f <- pf(Z^2 / ((n * (1/n + 1/m))), 1, n-1)
m.covers <- mean(covers)
(c(Simulated=m.covers, Theoretical=f, Z=(m.covers - f)/sd(covers) * sqrt(length(covers))))
$\endgroup$
6
  • $\begingroup$ You say that using t instead of z won't make much difference. I believe you but haven't checked yet. With small sample size, the two critical values can differ a lot and the t distribution is the correct way to compute the CI. Why do you prefer to use z?? $\endgroup$ – Harvey Motulsky Jun 14 '16 at 3:13
  • $\begingroup$ It's purely illustrative and $Z$ is simpler. When you use $t$ it is interesting that the curves in the figure start high and descend to their limit. In particular, the chance of reproducing a significant result is then much higher for small samples than for large! Note that there's nothing to check, because you are free to interpret $Z_{\alpha}$ as a percentage point of the appropriate Student t distribution (or of any other distribution you might care to name). Nothing changes in the analysis. If you do want to see the particular effects, uncomment the qt line in the code. $\endgroup$ – whuber Jun 14 '16 at 13:13
  • 1
    $\begingroup$ +1. This is a great analysis (and your answer has way too few upvotes for what it is). I just came across a paper that discusses this very question in great detail and I thought you might be interested: Cumming & Maillardet, 2006, Confidence Intervals and Replication: Where Will the Next Mean Fall?. They call it capture percentage of a confidence interval. $\endgroup$ – amoeba Dec 6 '16 at 22:42
  • $\begingroup$ @Amoeba Thank you for the reference. I especially appreciate one general conclusion therein: "Replication is central to the scientific method, and researchers should not turn a blind eye to it just because it makes salient the inherent uncertainty of a single study." $\endgroup$ – whuber Dec 6 '16 at 22:54
  • 1
    $\begingroup$ Update: Thanks to the ongoing discussion in the sister thread, I now believe my reasoning in the above comment was not correct. 95% CIs have 83% "replication-capture", but this is a statement about repeated sampling and cannot be interpreted as giving a probability conditioned on one particular confidence interval, at least not without further assumptions. (Perhaps both this and previous comments should better be deleted in order not to confuse further readers.) $\endgroup$ – amoeba Dec 14 '16 at 20:22
4
$\begingroup$

[Edited to fix the bug WHuber pointed out.]

I altered @Whuber's R code to use the t distribution, and plot coverage as a function of sample size. The results are below. At high sample size, the results match WHuber's of course.

enter image description here

And here is the adapted R code, run twice with alpha set to either 0.01 or 0.05.

sigma <- 2 
mu <- -4
alpha <- 0.01
n.sim <- 1e5
#
# Compute the multiplier.

for (n in c(3,5,7,10,15,20,30,50,100,250,500,1000))
{
   T <- qt(alpha/2, df=n-1)     
# Draw the first sample and compute the CI as [l.1, u.1].
#
x.1 <- matrix(rnorm(n*n.sim, mu, sigma), nrow=n)
x.1.bar <- colMeans(x.1)
s.1 <- apply(x.1, 2, sd)
l.1 <- x.1.bar + T * s.1 / sqrt(n)
u.1 <- x.1.bar - T * s.1 / sqrt(n)
#
# Draw the second sample and compute the mean as x.2.
#
x.2 <- colMeans(matrix(rnorm(n*n.sim, mu, sigma), nrow=n))
#
# Compare the second sample means to the CIs.
#
covers <- l.1 <= x.2 & x.2 <= u.1
#
Coverage=mean(covers)

print (Coverage)

}

And here is the GraphPad Prism file that made the graph.

$\endgroup$
5
  • $\begingroup$ I believe your plots do not use the t distribution, due to a bug: you set the value of T outside the loop! If you would like to see the correct curves, just plot them directly using the theoretical result in my answer, as given at the end of my R code (rather than relying on the simulated results): curve(pf(qt(.975, x-1)^2 / ((x * (1/x + 1/x))), 1, x-1), 2, 1000, log="x", ylim=c(.8,1), col="Blue"); curve(pf(qt(.995, x-1)^2 / ((x * (1/x + 1/x))), 1, x-1), add=TRUE, col="Red") $\endgroup$ – whuber Jun 17 '16 at 0:53
  • 1
    $\begingroup$ @whuber. Yikes! Of course you are right. Embarrassing. I've fixed it. As you pointed out the coverage is higher with tiny sample sizes. (I fixed the simulations, and didn't try your theoretical function.) $\endgroup$ – Harvey Motulsky Jun 17 '16 at 2:25
  • $\begingroup$ I am glad you fixed it, because it is very interesting how high the coverage is for small sample sizes. We could also invert your question and use the formula to determine what value of $Z_{\alpha/2}$ to use if we wished to assure (before doing any experiments), with probability $p=0.95$ (say), that the mean of the second experiment would lie within the two-sided $1-\alpha$ confidence interval determined from the second. Doing so, as a routine practice, could be one intriguing way of addressing some criticism of NHST. $\endgroup$ – whuber Jun 17 '16 at 13:28
  • $\begingroup$ @whuber I think the next step is to look at the distribution of coverage. So far, we have the average coverage (average of many first experiments, with average of many second experiments each). But depending on what the first experiment is, in some cases the average coverage will be poor. It would be interesting to see the distribution. I'm trying to learn R well enough to find out. $\endgroup$ – Harvey Motulsky Jun 17 '16 at 17:34
  • $\begingroup$ Regarding the distributions, see the paper I linked to in the comments above. $\endgroup$ – amoeba Dec 6 '16 at 22:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.