Explanation for the Wilson Score Interval? I'm looking at this blog to try to understand the Wilson Score interval. I understand it somewhat, but I'm confused by the part under the title "Excerpt". In particular, I don't understand what he's calling the "Interval equality principal" and how he arrived at the below graph:

Could someone elaborate on it, or really just explain how/why the Wilson Score Interval is arrived at from the basic Wald Interval (normal approximation)? Basically, what I'm trying to understand is why the Wilson Score Interval is more accurate than the Wald test / normal approximation interval?
 A: The explanation of "interval equality principle" was impossible for me to readily understand. However, it is not needed to know why the Wilson score interval works. The Wilson interval is derived from the Wilson Score Test, which belongs to a class of tests called Rao Score Tests. It relies on the asymptotic normality of your estimator, just as the Wald interval does, but it is more robust to deviations from normality. Case in point: Wald intervals are always symmetric (which may lead to binomial probabilties less than 0 or greater than 1), while Wilson score intervals are assymetric. 
Wilson intervals get their assymetry from the underlying likelihood function for the binomial, which is used to compute the "expected standard error" and "score" (i.e., first derivative of the likelihood function) under the null hypotheisis. Since these values will change as you very your null hypothesis, the interval where the normalized score (score/expected standard error) exceeds your pre-specified Z-cutoff for significance will not be symmetric, in general. 
In basic terms, the Wilson interval uses the data more efficiently, as it does not simply aggregate them into a a single mean and standard error, but uses the data to develop a likelihood function that is then used to develop an interval.
A: The difference between the Wald and Wilson interval is that each is the inverse of the other. This is how the Wilson interval is derived! 
As a result we have the following type of equality, which I referred to as the interval equality principle to try to get this idea across.

Wald(tail, α, Wilson(¬tail, α, p)) = p,

and, correspondingly,

Wilson(tail, α, Wald(¬tail, α, P)) = P, 

where tail ε {0=lower, 1=upper}, α represents the error level (e.g. 1 in 100 = 0.01), and p is an observed probability ε [0, 1]. The Wald interval is a legitimate approximation to the Binomial interval about an expected population probability P, but (naturally) a wholly inaccurate approximation to its inverse about p (the Clopper-Pearson interval). 
In fitting contexts it is legitimate to employ a Wald interval about P because we model an ideal P and compute the fit from there. But when we plot observed p, we need to employ the Wilson interval.
I would encourage people to read the paper, not just the excerpt! 
Sean Wallis
A: I think the plot in question originally comes from Wallis (2021) so I recommend you have a look at that book for further explanation on the particulars of that graphical representation.  In any case, the main reason why the Wilson score interval is superior to the classical Wald interval is that is is derived by solving a quadratic inequality for the proportion parameter that leads to an interval that respects the true support of the parameter.  Contrarily, the Wald interval can go outside the true support, and it also has worse coverage properties (see Brown, Cai and DasGupta (2001) for further discussion).
It might help here to show you the derivation of the interval in algebraic terms.  Suppose we have $n$ binary data values giving the sample proportion $p_n$ (which we will treat as a random variable) and let $\theta$ be the true proportion parameter.  The classical Wald interval uses the asymptotic pivotal distribution:
$$\sqrt{n} \cdot \frac{p_n-\theta}{\sqrt{\theta(1-\theta)}} \overset{\text{Approx}}{\sim} \text{N}(0,1).$$
For the Wilson score interval we first square the pivotal quantity to get:
$$n \cdot \frac{(p_n-\theta)^2}{\theta(1-\theta)} \overset{\text{Approx}}{\sim} \text{ChiSq}(1).$$
Let $\chi_{1,\alpha}^2$ denote the critical point of the chi-squared distribution with one degree-of-freedom (with upper tail area $\alpha$).  For any confidence level $1-\alpha$ we then have the probability interval:
$$\begin{align}
1-\alpha 
&\approx \mathbb{P} \Big( n (p_n-\theta)^2 \leqslant \chi_{1,\alpha}^2 \theta(1-\theta) \Big) \\[6pt]
&= \mathbb{P} \Big( n (p_n^2 - 2 p_n \theta + \theta^2) \leqslant \chi_{1,\alpha}^2 (\theta-\theta^2) \Big) \\[6pt]
&= \mathbb{P} \Big( (n + \chi_{1,\alpha}^2) \theta^2 - (2 n p_n + \chi_{1,\alpha}^2) \theta + n p_n^2  \leqslant 0 \Big) \\[6pt]
&= \mathbb{P} \Bigg( \theta^2 - 2 \cdot\frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \cdot \theta + \frac{n p_n^2}{n + \chi_{1,\alpha}^2} \leqslant 0 \Bigg) \\[6pt]
&= \mathbb{P} \Bigg( \bigg( \theta - \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \bigg)^2 \leqslant \frac{\chi_{1,\alpha}^2 (n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2)}{(n + \chi_{1,\alpha}^2)^2} \Bigg) \\[6pt]
&= \mathbb{P} \Bigg( \theta \in \Bigg[ \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \pm \frac{\chi_{1,\alpha}}{n + \chi_{1,\alpha}^2} \cdot \sqrt{n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2} \Bigg] \Bigg), \\[6pt]
\end{align}$$
and substitution of the observed sample proportion (for simplicity I will use the same notation for this value) then leads to the Wilson score interval:
$$\text{CI}_\theta(1-\alpha) = \Bigg[ \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \pm \frac{\chi_{1,\alpha}}{n + \chi_{1,\alpha}^2} \cdot \sqrt{n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2} \Bigg].$$
As you can see, solving the quadratic inequality in the probability interval leads to an interval that respects the true space of possible values of the proportion parameter (i.e., it is between zero and one).  This is a major advantage of this method but it also has better coverage properties in general.
