--- title: "Comparing Exact and Approximate Confidence Intervals" subtitle: "STAT/MATH 355: Statistical Theory" author: "Your Name" date: today format: html: toc: true embed-resources: true --- # Setup Let $X_1$,...,$X_n$ be from a Bernoulli distribution, so $$p_X(x_i; p) = p^{x_i}(1-p)^{1-x_i} \quad \text{for} \quad x_i = 0,1 \quad \text{and} \quad 0 < p < 1$$ We will use two different approaches to construct a confidence interval for $p$: - establishing an **approximate** pivotal quantity for $p$ - using our estimator's **exact** dsitribution ## Approximate CI See the course website for exercises walking you through the derivation of an **approximate** 95% CI for $p$. Write down the formula for that CI here: $$\hat{p} \pm ???$$ ## Exact CI Since $X_1, \dots, X_n \overset{iid}{\sim} \text{Bernoulli}(p)$, we know $X = \sum_{i=1}^n X_i \sim \text{Binomial}(n, p)$. The **Clopper-Pearson Interval** uses this result to construct an **exact** CI for $p$. Suppose we observe $X = k$ successes. The Clopper-Pearson CI for $p$ is constructed as follows: - Find the _largest_ $p$ such that $P(X \leq k) \geq \frac{\alpha}{2}$. Call this largest value $p_U$. - Find the _smallest_ $p$ such that $P(X \geq k) \geq \frac{\alpha}{2}$. Call this smallest value $p_L$. - Define the $100(1-\alpha)\%$ confidence interval for $p$ to be: $(p_L, p_U)$. \ # Exercises ## Part 1 (@) Read through the CI construction process described above. Turn to a neighbor and try to articulate the logic behind this process. Why does this yield a valid CI? > YOUR ANSWER HERE (@) Suppose we observe 3 successes in 10 trials. Find the Clopper-Pearson 95% CI for $p$. Use R's `pbinom` function to calculate $P(X \le k)$ = `pbinom(k, size = n, prob = p)` and $P(X \ge k) = 1 - P(X < k) = 1 - P(X \le k - 1)$ = `...` You can use a trial-and-error/grid search approach, or define a function that reflects equality in those inequalities and use R's `uniroot` to find the root of your function. Starter code for the latter is provided below. ```{r exact-CI} #| eval: false ## to use uniroot, we need to set up a function f ## where f(p) = 0 when p is a soln to the eqn we want to solve ## use uniroot to solve for upper bound ## find p such that P(X <= k) = alpha/2 or P(X <= k) - alpha/2 = 0 f_upper <- function(p){ pbinom(3, size = 10, prob = p) - 0.05/2 } pU <- uniroot(f_upper, interval = c(0,1))$root ## search over [0,1] since p is probability ## now use uniroot to solve for lower bound ## find p such that P(X >= k) = alpha/2 or P(X >= k) - alpha/2 = 0 f_lower <- function(p){ ## insert code here } pL <- uniroot(f_lower, interval = c(0,1))$root ## print out CI CI <- c(pL, pU) print(CI) ``` (@) Check your answer to the previous question by running the `binom.test` function, which is a pre-existing function in R that can be used to get an exact confidence interval for $p$. ```{r check-exact-CI} ## get exact 95% CI for p using binom.test binom.test(x = 3, n = 10, conf.level = 0.95) ``` (@) Compare the Clopper-Pearson CI to the approximate CI you derived earlier. Comment on any similarities or differences that you notice. ```{r approximate-pivot-CI} ## your code here ``` > YOUR ANSWER HERE (@) How might you go about deciding which of these CIs to report? > YOUR ANSWER HERE \ ## Part 2 With confidence intervals, our goal is to construct an interval that contains the true parameter as an interior point $100(1-\alpha)$% of the time, ie $$P(LB \le \theta \le UB) = 1-\alpha$$ It's worth confirming whether this is actually true: when we construct a 95% CI according to some procedure, do those intervals actually contain the true parameter 95% of the time? To check this, we need to design and implement a **simulation study**: 1. generate some data where we know the true $\theta$ 2. construct a CI based on this data 3. check if $\theta$ is in this interval 4. repeat steps 1--3 many times 5. across all these reps, how often does the CI actually cover the true $\theta$? This simulation approach yields an empirical estimate of the *coverage* of the CI---we'd like coverage to be close to $100(1-\alpha)$%. Let's perform a simulation study to evaluate the coverage of the exact and approximate CIs we constructed above! (@) Insert code into the `do_one_p` function, below, to calculate an exact and approximate $100(1-\alpha)$ percent confidence interval for $p$ using your work from above. **NOTES**: - I've given you most of the code here. Take your time looking through this to make sure you understand what each line is doing. ASK QUESTIONS about any parts that don't make sense. - Feel free to update as you see fit, especially if you see ways to make this code more efficient, easier to read, etc.! ```{r coverage-simulations} #| eval: false ## function to do one simulation replicate do_one_p <- function(n, p, alpha = 0.05){ # generate sample of size n from Binom(1,p) x <- rbinom(n, size = 1, prob = p) phat <- mean(x) # exact CI exactCI <- ## insert code here # approximate CI approxCI <- ## insert code here # check if true p is in CIs inExact <- (p >= exactCI[1]) & (p <= exactCI[2]) inApprox <- (p >= approxCI[1]) & (p <= approxCI[2]) # return results return(c(inExact, inApprox)) } ## run simulation once to check what code is doing ## (using n = 10, p = 0.3) do_one_p(n = 10, p = 0.3) ## repeat simulation 1000 times ## (using n = 10, p = 0.3) set.seed(455) simres <- replicate(1000, do_one_p(n = 10, p = 0.3)) ## calculate coverage for each type of CI ## insert code here ``` (@) What if we increased $n$? Would you expect the coverage of either interval to be impacted? Write down your hypothesis, and then try running your simulation again with different values of $n$ to test your hypothesis. > YOUR HYPOTHESIS ```{r more-simulations} ## repeat the process above, with different values of n ``` ```{r plot-results} ## make a plot summarizing how coverage changes with n ``` > YOUR CONCLUSIONS (@) Repeat the process above, this time keeping $n$ fixed but changing $p$. What do you learn about how the true probability of success impacts the performance of these confidence intervals? ```{r} ## your code here ``` > YOUR ANSWER HERE