Why is the confidence interval for difference in two proportions inconsistent with the p-value for this difference in R prop.test()?

Question

I calculated the difference in 2 proportions using the prop.test() using some exemplary data I found in internet.

> (ppt <- prop.test(x = c(11, 8), n = c(16, 21),correct = FALSE))
2-sample test for equality of proportions without continuity correction


data:  c(11, 8) out of c(16, 21)
X-squared = 3.4159, df = 1, p-value = 0.06457
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.001220547  0.614315785
sample estimates:
   prop 1    prop 2 
0.6875000 0.3809524 

The p-value is 0.06457. The confidence interval for the difference in proportions is -0.001220547 to 0.614315785.

If I calculate the p-value from the CI manually using the normal distribution:

> z_score <- qnorm((1 + 0.95) / 2)
> se <- diff(ppt$conf.int) / (2 * z_score)
> z <- pw$estimate / se
> 2 * (1 - pnorm(abs(z)))
[1] 0.05091551

that previous p-value does not agree with the obtained p-value using the provided confidence interval.

When I calculate it from the logistic regression using the marginal effects, I get the same confidence interval that prop.test() gives and the p-value I calculated above.

data <- data.frame(Status = c(rep(TRUE, 11), rep(FALSE, 16-11), rep(TRUE, 8), rep(FALSE, 21-8)), Group = c(rep("Gr1", 16), rep("Gr2", 21)))
> m <- glm(Status ~ Group,family = binomial(), data=data)
> margins::margins_summary(m)
factor     AME     SE       z      p   lower  upper
 GroupGr2 -0.3065 0.1570 -1.9522 0.0509 -0.6143 0.0012

But when I use ANOVA on this model with Rao test I get the p-value from prop.test()

> anova(m, test="Rao")
Analysis of Deviance Table
Model: binomial, link: logit
Response: Status
Terms added sequentially (first to last)
  Df Deviance Resid. Df Resid. Dev    Rao Pr(>Chi)  

NULL                     36     51.266

Group  1   3.4809        35     47.785 3.4159  0.06457 .

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So how is that possible, that prop.test() gives me p-value which does NOT come from the confidence interval provided by the same function?

It looks like the confidence interval is Wald's, but the p-value is Rao score test. What is the sense behind it? I was taught that confidence interval should match the p-value. Here it does not, as they are calculated two different methods. It agrees in this scenario in statistical significance, but what if they do not?

---

EDIT: My question is about the possible statistical reasons for using different methods for calculating the p-value and confidence interval in one method. I guess there are some considerations that statisticians agree on here, so I'm curious why this discrepancy is justified?

---

EDIT: The linked answer: P value and confidence interval for two sample test of proportions disagree does not agree with my findings. It's clear, that the confidence interval I got is NOT THE WILSON CI, it's the WALD, so the linked answer is not correct. Again, my question is what is the reason behind reporting WALD confidence interval and Rao score test?

The proposed answer, which led to close my post is wrong for 2 proportions, which is evidently visible in my findings. This can be found also here: https://stats.stackexchange.com/a/570528/382831

Another proof:

> PropCIs::wald2ci(11, 16, 8, 21, conf.level=0.95, adjust="Wald")
data:
95 percent confidence interval:
 -0.001220547  0.614315785
sample estimates:
[1] 0.3065476

---

EDIT 3: It's NOT about using different distributions (it's the z statistic for the proportions, NOT t statistic, because of the s=p(1-p) mean-variance dependency!), as some people suggested in the comments!

Moreover, I didn't do any mistake in the formula:

se <- diff(ppt$conf.int) / (2 * z_score)

where 2 is NOT 1.96 or anything else rounded up, as some suggested that I use some "approximation". It's 1/2 of the CI width, "half-CI", the "precision" divided by the z_score: precision/z_score = (CI_length/2) / z_score = CI_length/(2*z_score).

Answer 1

The difference is in how the standard error is computed. This can be done based on a pooled estimate of the variance and the assumption of equal probabilities (which makes sense if the null hypothesis is true and what you use when you compute the p-value), or based on non-pooled variance estimate.

$$z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{ \left(\frac{\hat{p}_1(1-\hat{p}_1)}{n_1} + \frac{\hat{p}_2(1-\hat{p}_2)}{n_2} \right)}} = \frac{\frac{11}{16} - \frac{8}{21}}{\sqrt{ \left(\frac{(11/16)\cdot(5/16)}{16} + \frac{(8/21)\cdot(13/21)}{21} \right)}} = 1.952191$$

$$z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p}) \left(\frac{1}{n_1} + \frac{1}{n_2} \right)}} = \frac{\frac{11}{16} - \frac{8}{21}}{\sqrt{\frac{19 \cdot 18}{37^2} \left(\frac{1}{16} + \frac{1}{21} \right)}} = 1.848227$$

These relates to respectively a p-value of 0.05092 and 0.06457

> z = (11/16 - 8/21)/sqrt((11/16)*(5/16) * 1/16+ (8/21)*(13/21) * 1/21)
> 2*(1-pnorm(z))
[1] 0.05091551
> 
> z = (11/16 - 8/21)/sqrt((19/37)*(18/37) * (1/16+1/21))
> 2*(1-pnorm(z))
[1] 0.06456946

---

Also note the references in the prop.test manual

Newcombe R.G. (1998). Two-Sided Confidence Intervals for the Single Proportion: Comparison of Seven Methods. Statistics in Medicine, 17, 857–872. doi: 10.1002/(SICI)1097-0258(19980430)17:8<857::AID-SIM777>3.0.CO;2-E.

Newcombe R.G. (1998). Interval Estimation for the Difference Between Independent Proportions: Comparison of Eleven Methods. Statistics in Medicine, 17, 873–890. doi: 10.1002/(SICI)1097-0258(19980430)17:8<873::AID-SIM779>3.0.CO;2-I.

There are a lot of different approaches possible that may give a different result.