Conf. Interval coverage for a Bernoulli parameter p, and its log-odds transform.

bern_ci.jl

function sim_coverage(p, n, nsims)
    # True log-odds/relative-risk
    rr = log(p) - log(1-p)

    # Simulate out data
    dist = Bernoulli(p)
    xs = [rand(dist, n) for _ in 1:nsims]

    # Bernoulli param estimate and its std. err.
    ps = [mean(x) for x in xs]
    se_p = [sqrt(p*(1-p)/n) for p in ps]

    # Log RR  and its std. err. (via delta method approx.)
    logrrs = [log(p) - log(1-p) for p in ps]
    se_logrr = [sqrt(1/sum(x) + 1/(n-sum(x))) for x in xs]

    # Normal-based confidence intervals
    z1 = quantile(Normal(0,1), .025)
    z2 = quantile(Normal(0,1), .975)
    cis_p = [(p + z1*se, p + z2*se) for (p, se) in zip(ps, se_p)]
    cis_rr = [(rr + z1*se, rr + z2*se) for (rr, se) in zip(logrrs, se_logrr)]

    # Coverage rates: how often does CI contain true values?
    coverage_p = mean([lo <= p && p <= hi for (lo, hi) in cis_p])
    coverage_rr = mean([lo <= rr && rr <= hi for (lo, hi) in cis_rr])

    # Return coverage rates
    coverage_p, coverage_rr
end


const nsims = 100_000
@printf("%5s %5s %10s\n", " ", " ", "Coverage")
@printf("%5s %5s %10s\n", " ", " ", "--------")
@printf("%5s %5s %5s %5s\n", "p", "n", "p", "logrr")
@printf("%5s %5s %5s %5s\n", "-----", "-----", "-----", "-----")

for p in [0.005, 0.01, .5]
    println("")
    for n in [20, 50, 150, 500, 1_000]
        cov_p, cov_rr = sim_coverage(p, n, nsims)
        @printf("%5.3f %5d %5.3f %5.3f\n", p, n, cov_p, cov_rr)
    end
end


#               Coverage
#               --------
#     p     n     p logrr
# ----- ----- ----- -----
# 
# 0.005    20 0.094 0.000
# 0.005    50 0.222 0.196
# 0.005   150 0.530 0.489
# 0.005   500 0.914 0.877
# 0.005  1000 0.870 0.963

# 0.010    20 0.182 0.166
# 0.010    50 0.397 0.308
# 0.010   150 0.779 0.715
# 0.010   500 0.873 0.963
# 0.010  1000 0.927 0.963

# 0.500    20 0.959 0.959
# 0.500    50 0.935 0.967
# 0.500   150 0.939 0.959
# 0.500   500 0.946 0.946
# 0.500  1000 0.946 0.954