Updating population variance assumption after sampling

Question

I'm trying to get my head around how to think about a situation where an assumption on the population variance in the sample sizing stage of an experiment goes wrong.

Let's say I'm estimating a binomial proportion and I want a 95% confidence interval to have an error of ±5 percentage points. I can estimate the required sample size under the most conservative assumption on population variance (the case where what I'll denote π below, the true population proportion of successes, is 50% or 0.5) using the standard formula:

julia> using Distributions

julia> function sample_size(;π = 0.5, α = 0.95, ε = 0.05)
           Χ² = quantile(Chisq(1), α)

           n = round(Int, Χ² * π * (1 - π) / ε^2, RoundUp)
       end
sample_size (generic function with 1 method)

julia> sample_size()
385

But now let's assume I was overly conservative and the true population proportion is actually only 10%, so the population variance is materially lower. In that case, my sample size formula tells me that I could have gotten away with a much smaller sample:

julia> sample_size(π = 0.1)
139

but it's too late, I've already sampled. Now I could invert the sample size formula to calculate the new (narrower) width of a 95% confidence interval under the new assumptions, or Monte Carlo it to find:

julia> quantile([mean(rand(Bernoulli(0.1), 385)) for _ ∈ 1:1_000_000], [0.025, 0.975])
2-element Vector{Float64}:
 0.07012987012987013
 0.12987012987012986

So at 95% confidence I would expect a ±3 percentage point error around the estimated population success rate when the true population rate is 0.1. I can look at my sample:

julia> mean(rand(Bernoulli(0.1), 385))
0.0987012987012987

and now conclude what? Based on the design of my experiment I could faithfully report "the estimated rate of successes in the population is 9.8%, with a 95% confidence interval from 4.8% to 14.8%". But of course this is based on assuming that the population rate is 50%, and I've drawn 9.8%, and:

julia> quantile([mean(rand(Bernoulli(0.5), 385)) for _ ∈ 1:1_000_000], [0.0001, 0.9999])
2-element Vector{Float64}:
 0.4051948051948052
 0.5948051948051948

it's basically impossible to get only 37 successes in 385 draws if there really are 50% successes in my population, so I know (or am at least very very confident in assuming) that the actual rate of successes is not 50%.

What do I do with this knowledge?

I suppose if I was Bayesian I could just say I start from $Beta(1,1)$ and update this by adding 1 to $\alpha$ or $\beta$ depending on whether I got a success, in which case I'd get:

julia> quantile(Beta(1 + 385/2, 1 + 385/2), [0.025, 0.975])
2-element Vector{Float64}:
 0.4502761549342859
 0.549723845065714

if my sample contained 50% successes and

julia> quantile(Beta(1 + 37, 1 + 385-37), [0.025, 0.975])
2-element Vector{Float64}:
 0.07061116109463793
 0.12970001401043718

in the example where I only get 37 successes in 385 draws - happy days, but how does that square with the frequentist assumption driven approach?

Answer 1

1. As a general rule, it does not matter at all what your assumptions were before you collected your data. The data is the data, and you should let it tell its story. E.g. you may have assumed normality of a dataset before sampling; but if, after collecting the data you realize the data is severely skewed, you need to adjust accordingly. Yes, the test you wanted to use (e.g. 1-sample t) may no longer be the right one, the power of your test may no longer be what you wanted, but... You may mention this in the discussion section (e.g. that your test was under powered because initial assumptions were not met; but then collect more data). But, in your specific instance, it is the other way around; your test is even better powered than you thought. So nothing to worry about, nor to disclose. Just deal with the data as it is (not as you thought it could be, which in your case was a worse case scenario).
2. You seem to want to use normal approximations to estimate the CI of a binomial proportion. But at 10% success, that binomial PMF does not look normal at all. So you should instead use an "exact" binomial CI; Clopper-Pearson being the "treaditional" one, but you could also use Blaker, which is also "exact" (i.e. based on the exact binomial PMF), but a bit less conservative than Clopper-Pearson in the double sided case.
3. Typically, when dealing with CI's of proportions, one is concerned only about one tail. Is the success rate too low? (or is the failure rate too high?), but rarely both. But you seem to care about both? That will give you bounds which are larger than maybe you wanted? (e.g. the low bound of the success rate will really be a 97.5% bound, not 95%). Just food for thought...
4. For example, using Clopper-Pearson, one gets a 95% double sided CI of [6.86%, 13.00%], and a single sided low bound of 7.25% (since they are "successes, we care about the low bound). Blaker 95% double sided gives [6.91%, 12.90%]. Now if you use Agresti-Coull (which seems to be what you did? But I am not familiar with Julia...) you get [7.06%, 12.99%], or 7.40% for a single-sided low bound. Depending on your context, these (small) differences may be inconsequential, or they may be meaningful.