Confidence intervals for a population proportion: Why do we use a z-distribution (instead of a t-distribution) even who we estimate $\sigma_{\hat{p}}$ When deriving a confidence interval for a population mean $\mu$ where we do not know the value of the population standard deviation $\sigma$, we use the sample standard deviation $s$ to estimate $\sigma$ and so our margin of error ends up being related to a $t$-distrubition:
$E = t_c \frac{s}{\sqrt{n}}$
Now, when deriving a confidence interval for an unknown population proportion $p$, we substitute $\hat{p}$ for $p$ into the formula for the standard deviation of $\hat{p}$:
$\sigma_{\hat{p}}=\sqrt{\frac{p(1-p)}{n}} \rightarrow \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$
However, when we make this substitution, for some reason we still say that the margin of error on our resulting confidence interval depends on a $z$-distribution!!!
$E = z_c \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$

I don't understand why this is. In both cases, we are using some type of sample standard deviation to estimate the population standard deviation. Why in the former case (for a population mean) do we use a $t$ distribution (which I feel is very appropriate), but in the case of a population proportion we still use a $z$-distribution??
Insight appreciated.

Here's me asking (what I view) to be the exact same question but in a different way:
Does $\frac{\hat{p}-p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}}$ ever have a t-distribution? It seems like it should since we are using $\hat{p}$ for $p$ in the formula for $\sigma_{\hat{p}}$; this closely resembles how we get to a t distribution for sample means (by substituting in $s$ for $\sigma$).
 A: To really understand what is going on, let's consider some examples.
Example 1. Suppose I have a sample whose observations are known to arise from a random process that generates normally distributed outcomes.  However, we do not know the mean of this normal distribution, but we do know its variance.  That is to say, $$X_i \sim \operatorname{Normal}(\mu, \sigma^2)$$ and we know $\sigma^2$ but do not know $\mu$.  A natural question that arises is the use of observed data, say the sample $$\boldsymbol X = (X_1, X_2, \ldots, X_n),$$ to construct an inference about $\mu$.  Several types of inference could be made; e.g., a point estimate, say $\bar X$.  Or we might wish to construct an interval estimate, say $$\left( \bar X - z^*_{\alpha/2} \frac{\sigma}{\sqrt{n}}, \bar X + z^*_{\alpha/2} \frac{\sigma}{\sqrt{n}}\right). \tag{1}$$  (Note that this is a valid estimator because $\sigma$ is known.)  Or we can perform a hypothesis test of some kind, such as $$H_0 : \mu = \mu_0 \quad \text{vs.} \quad H_a : \mu \ne \mu_0.$$  Such a test might employ the statistic $$Z \mid H_0 = \frac{\bar X - \mu_0}{\sigma/\sqrt{n}} \sim \operatorname{Normal}(0,1). \tag{2}$$  In the case of the confidence interval and the hypothesis test, we are using information about the variability of $X$.  Because $\sigma^2$ is known, the statistics in $(1)$ and $(2)$ can be computed from from the sample.  The fact that $(2)$ is standard normal is not an approximation:  it is a mathematical consequence of the fact that the sum of normally distributed random variables is itself normally distributed; i.e., $$X_1 + X_2 + \cdots + X_n \sim \operatorname{Normal}(n \mu, n \sigma^2).$$  Because we posited that the $X_i$ are known to be normal, the statistic $(2)$ is standard normal under the assumption that the null is true; and if the null is false, it is still normal with unit variance, but the mean is some other nonzero number.
To be clear, even if the sample size is $n = 1$, you could still construct $(1)$ and perform the test $(2)$.  Sample size in this case is irrelevant.
Example 2.  Suppose we have the same data as in Example 1 above, but now we neither know the mean nor the variance.  In such a situation, the interval $(1)$ cannot be calculated because it depends on an unknown parameter $\sigma$; similarly, the test statistic $(2)$ is invalid.  It is natural to want to use an estimator for $\sigma$ that is calculated from the sample:  $$s^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar X)^2,$$ the Bessel-corrected unbiased sample variance.  That is to say, we are using the data to estimate the variability that we do not know in advance.  However, in doing this, the statistic $$T \mid H_0 = \frac{\bar X - \mu_0}{s/\sqrt{n}} \tag{3}$$ is not normally distributed:  it was shown that in fact, $(3)$ has what is now known as a Student's $t$-distribution with $n-1$ degrees of freedom.  Again, this is not an approximation for the same reason that $Z \mid H_0$ in Example 1 is exactly normal:  we know the $X_i$ are normally distributed.  A mathematical derivation of this fact is available but is lengthy, hence is not provided here.
By inverting the test statistic $(3)$, we can derive a $100(1-\alpha)\%$ confidence interval for $\mu$ as $$\left(\bar X - t_{n-1,\alpha/2}^* \frac{s}{\sqrt{n}}, \bar X + t_{n-1,\alpha/2}^* \frac{s}{\sqrt{n}}\right), \tag{4}$$ where we have replaced $\sigma$ with $s$, which requires us to replace the standard normal quantile $z_{\alpha/2}^*$ with the Student $t$-quantile $t_{n-1,\alpha/2}^*$ which represents the upper $100(\alpha/2)$ percentile of the Student's $t$-distribution with $n-1$ degrees of freedom.
The sample size in this situation only tells us the degrees of freedom.  Exact inference remains possible no matter the sample size.  If $n$ is very small, the degrees of freedom is small and the quantile $t_{n-1,\alpha/2}^*$ will reflect this. If $n$ is very large, $t_{n-1,\alpha/2}^*$ will be very close to $z_{\alpha/2}^*$.
Example 3.  Now suppose that the $X_i$ are not necessarily normally distributed.  Could we still use the estimates and statistics $(1)-(4)$?  The answer to that question depends on the sample size and the extent of deviation of the $X_i$ from normality.  If the sample size is small, the Central Limit Theorem does not apply, so the use of these as an approximation is potentially problematic, with the possible loss of the nominal coverage probability or an increase in the Type I error.
When the sample size is large, the Central Limit Theorem starts takes effect and we can see that $\bar X$ will tend toward a normal distribution (except in certain unusual cases).  Still, if the variance is unknown and the deviation from normality not extreme, the $t$-statistic will be the more appropriate choice rather than a $z$-statistic:  that is to say, the choice between $z$ or $t$ is driven by whether the variance is known.
Example 4.  Now suppose we have a sample drawn from a Bernoulli distribution, i.e. $$X_i \sim \operatorname{Bernoulli}(p), \quad \Pr[X_i = 1] = p, \quad \Pr[X_i = 0] = 1-p.$$  Here, $p$ is the sole unknown parameter, and now the inference is for $p$.  A test statistic we can use for $$H_0 : p = p_0 \quad \text{vs.} \quad H_a : p \ne p_0$$ could be constructed from the binomial random variable $$B = \sum_{i=1}^n X_i \sim \operatorname{Binomial}(n, p).$$ Under the null hypothesis, $$B \mid H_0 \sim \operatorname{Binomial}(n, p_0) \tag{5}$$ is a statistic.  Then the rejection region for the test is determined by computing $$\begin{align}
R_U &= \operatorname{argmin}_x (\Pr[X \ge x \mid H_0] \le \alpha/2), \\
R_L &= \operatorname{argmax}_x (\Pr[X \le x \mid H_0] \le \alpha/2). \end{align}$$
So for instance, we choose the upper rejection limit $R_U$ to be the smallest nonnegative integer $x$ such that $\Pr[X \ge x \mid H_0]$ is at most half the Type I error $\alpha/2$.  This ensures that the probability of incorrectly rejecting $H_0$ when it is true, is at most $\alpha/2$.
As you can imagine, in practice, this is not simple to calculate unless $n$ is relatively small, since it involves computing the CDF of a binomial random variable.  Fortunately, in the case where $n$ is sufficiently large, as we mentioned in Example 3 above, the Central Limit Theorem applies:  $B$ is binomial with mean $np$ and variance $np(1-p)$, therefore the sample mean $\bar X$ can be approximated by a normal distribution with mean $\bar X/n = p$ and variance $p(1-p)/n$.  But since we do not know $p$, to construct a test statistic, we can choose for example $$Z \mid H_0 = \frac{\bar X - p_0}{\sqrt{\bar X (1-\bar X)/n}}, \tag{6a}$$ which uses the sample mean to to estimate the true variance, or $$Z \mid H_0 = \frac{\bar X - p_0}{\sqrt{p_0 (1-p_0)/n}}, \tag{6b}$$ which uses the hypothesized mean to estimate the true variance.  Both are statistics under the assumption that the null is true.  Both are approximations because the exact sampling distribution is discrete and arises from a binomial random variable.  But a $t$-distribution does not apply here because in either case, what we are using is an asymptotic approximation based on the Central Limit Theorem.
The confidence interval arising from inverting $(6a)$ is called the Wald interval and is also an approximation:  $$\left(\bar X - z_{\alpha/2}^* \sqrt{\frac{\bar X (1 - \bar X)}{n}}, \bar X + z_{\alpha/2}^* \sqrt{\frac{\bar X (1 - \bar X)}{n}}\right). \tag{7a}$$  Inverting $(6b)$ yields a different interval called the Wilson score interval which we will not write here since it is a bit long; see the link for more details.  In either case, these intervals, like the hypothesis test, rely on the CLT for validity; as such, when the sample size is small, these may lead to incorrect inference, and the exact test/intervals should be used instead.
In closing, the answer to your question is that the reasons for a $z$- or $t$-distribution arise from different issues with the relevant statistics, as described in the above examples:  the use of Student's $t$ is a mathematical consequence of the way $Z/\sqrt{V/n}$ is distributed, where $Z$ is standard normal, $V$ is chi-squared with $n$ degrees of freedom, and $Z$ and $V$ are independent.  It is not a consequence of an approximation or the use of the Central Limit Theorem.
A: By the Central Limit Theorem, we know that
$$ \sqrt{n}\left(\bar{X} - \mu\right) \to N(0, \sigma^2).$$
This is limit as $n$ goes to infinity. Therefore, even though we can use it to construct confidence intervals for population means (or population proportion, which is just an average of zeroes and ones), it is only an approximation.
Sometimes, when we have more information about the data, we can do better. It turns out that if the data follows a normal distribution, we have the following:
$$ \frac{\bar{X} - \mu}{\widehat{\sigma}/\sqrt{n}} \sim t_{n-1}.$$
In other words, because we know the exact distribution of the ratio $\frac{\bar{X} - \mu}{\widehat{\sigma}/\sqrt{n}}$, we can construct an exact confidence interval for the population mean.

One more thing to point out: combining the Central Limit Theorem and Slutsky's lemma (which I know is a more advanced concept), we also have
$$ \frac{\bar{X} - \mu}{\widehat{\sigma}/\sqrt{n}} \to N(0, 1).$$
In other words, the Central Limit Theorem also holds when we replace $\sigma$ with $\widehat{\sigma}$. It just so happens that when we know more about the data (e.g. data is normally distributed), we can be more precise and write down the exact distribution. When the data is not normally distributed (e.g. when it follows a Bernoulli distribution, such as when you have a proportion), the best we can do is often only rely on an approximation.

Edit: To answer your edit: no, the ratio $\frac{\hat{p}-p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}}$ never has a t-distribution. The t-distribution is defined as the ratio of a standard normal and the square root of a chi-square divided by its degrees of freedom (see for example this Wikipedia article). As you can see, the fact that we get a t-distribution is quite "rare" in a sense: it requires that the data be normally distributed so that $\sqrt{n}(\bar{X} - \mu)/\sigma$ is a standard normal and $(n-1)\widehat{\sigma}^2/\sigma^2$ follows a chi-square with $n-1$ degrees of freedom.
I think the confusion comes from the fact that see the t-distribution has naturally arising from replacing the unknown variance $\sigma^2$ with the sample variance $\widehat{\sigma}^2$. In fact, the t-distribution is a very special case. In most other cases, the best we can do is use the CLT approximation and use the $z_\alpha$ values.
