The z-test to compare two proportions is $\newcommand{\p}{\hat{p}}\newcommand{\v}{\mathrm{Var}} z=\frac{\p_1-\p_2}{\sqrt{\v(\p_1-\p_2)}}$. Usually it is defined that
$$\v(\p_1-\p_2)=\p(1-\hat{p})(1/n_1+1/n_2),$$
where
$$\p=\frac{n_1 \p_1+n_2 \p_2}{n_1+n_2}.$$
Is there any written reference that legitimizes me instead to use the unpooled variance, that is
$$\v(\p_1-\p_2)=\frac{\p_1(1-\p_1)}{n_1}+\frac{\p_2(1-\p_2)}{n_2}?$$
There is quite a bit of discussion about this on the AP site.
You can use whatever statistic you want, provided that you are clear about what you do and look at the appropriate null distribution to calculate p-values or thresholds.
But some statistics are better than others; in this case you'd be looking for (a) null distribution easily calculated and (b) power to detect difference.
But I don't know why you'd favor the unpooled variance over the pooled variance for the test, though it could be preferred in calculating a confidence interval for the difference.
The unpooled variance tends to be too small. This is because under the null hypothesis there will still be chance variation in the two observed proportions, although the underlying probabilities are equal. This chance variation contributes to the pooled variance but not to the unpooled variance.
As a result, $z$ for the unpooled statistic does not even approximately have a standard normal distribution. For instance, when $n_1 = n_2$ and the true probabilities are both $1/2$, the variance of $z$ is only $1/2$ instead of $1$. By using tables of the standard normal distribution, you will get incorrect p-values: they will tend to be artificially small, too often rejecting the null when the evidence is not really there.
Nevertheless, one wonders whether this could be corrected. It can. The question becomes whether a corrected value of $z$, based on unpooled estimates, could have greater power to detect deviations from the null hypothesis. A few quick simulations suggest this is not the case: the pooled test (compared to a properly adjusted unpooled test) has a better chance of rejecting the null whenever the null is false. Therefore I haven't bothered to work out the formula for the unpooled correction; it seems pointless.
In summary, the unpooled test is wrong, but with an appropriate correction, it can be made legitimate. However, it appears to be inferior to the pooled test.
Consider the estimated difference of proportions $\hat{d} = \hat{p}_1 - \hat{p}_2$. There is nothing fundamentally wrong with using the following unpooled estimate of the variance of $\hat{d}$:
$$ \hat{V}_{U} = \frac{\hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_2} $$
In fact, it is frequently used when constructing confidence intervals. Most elementary textbooks suggest the following approximate $100(1 - \alpha)$ confidence interval for $p_1 - p_2$:
$$\hat{d} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_2}} $$
which directly involves $\hat{V}_U$. However, when testing $p_1 = p_2$, the following pooled estimator is often preferable:
$$ \hat{V}_{P} = \hat{p}(1 - \hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)$$
The pooled estimator is only valid when $p_1 = p_2$, which is why it can't be used while constructing confidence intervals.
Let $n = n_1 + n_2$. The true variance of $\hat{d} = \hat{p}_1 - \hat{p}_2$ is equal to
$$\begin{aligned}V(\hat{d}) &= Var(\hat{p}_1) + Var(\hat{p}_2) \\ &= \frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2} \end{aligned}$$
When $p_1 = p_2 = p$, this reduces to
$$V(\hat{d}) = p(1-p)\left(\frac{1}{n_1} + \frac{1}{n_2}\right) $$
The pooled estimator of the variance of $\hat{d}$ has bias equal to
$$E[\hat{V}_P] - V[\hat{p}] = \frac{-p(1 -p)}{n}\left(\frac{1}{n_1} + \frac{1}{n_2} \right)$$
Similarly, the unpooled estimator of the variance has bias equals to
$$E[\hat{V}_U] - V[\hat{p}] = -p(1-p)\left(\frac{1}{n_1^2} + \frac{1}{n_2^2} \right)$$
Both biases go to zero as $n \to \infty$, but because $n = n_1 + n_2$, the bias of $\hat{V}_P$ will always be less than $\hat{V}_U$. Specifically, when $n_1 = n_2$, the pooled estimate of the variance will be half as biased as the unpooled estimator.
It is unclear which test rejects more often
The z-test statistic for the unpooled test will be more extreme than the pooled z-test statistic when
$$\frac{|\hat{d}|}{\sqrt{\hat{V}_P}} = \frac{|\hat{d}|}{\sqrt{\hat{p}(1 - \hat{p})(\frac{1}{n_1} + \frac{1}{n_2})}} < \frac{|\hat{d}|}{ \sqrt{ \frac{ \hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{ \hat{p}_2 (1 - \hat{p}_2)}{n_2} } } = \frac{|\hat{d}|}{\sqrt{\hat{V}_U}}$$
or equivalently, when the estimated variance of the unpooled test is smaller than that of the pooled test,
$$ \hat{p}(1 - \hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right) > \frac{ \hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{ \hat{p}_2 (1 - \hat{p}_2)}{n_2} $$
It is unclear when this holds. However, for samples of equal size, $n_1 = n_2 = 0.5$, this simplifies to
$$\hat{p}(1 - \hat{p}) > \frac{\hat{p}_1 (1 - \hat{p}_1) + \hat{p}_2 (1 - \hat{p}_2)}{2} $$
Because $f(x) = x(1-x)$ is strictly concave, the previous inequality will be always true. Therefore, for equal sample sizes, the unpooled test always reject as or more often than the pooled test.
Simulations
To simulate the two estimators, I ran 100 millions simulations in R with low sample sizes ($n_1 = n_2 = 20$, $p_1 = p_2 = 0.5$). Under these parameters, the true variance of $\hat{d}$ equals $0.025$.
| Parameter | Pooled | Unpooled |
|---|---|---|
| Variance of $\hat{z}$ | 1.0255 | 1.1146 |
| Mean of $\hat{V}$ | 0.0244 | 0.0237 |
| Rejection rate when $\alpha$ = 0.05 | 0.0425 | 0.0807 |
We see that the bias of the unpooled variance is twice as large as the bias of the pooled variance, leading to a $z$ statistic with a higher variance than a standard normal (1.115 vs 1) and that this leads to a higher null rejection rate than specified ($\alpha$ = 0.05).
However, as we increase the sample size, this difference vanishes. If $n_1 = n_2 = 100$, we get a true standard error of 0.005, and the following results:
| Parameter | Pooled | Unpooled |
|---|---|---|
| Variance of $\hat{z}$ | 1.005 | 1.020 |
| Mean of $\hat{V}$ | 0.004975 | 0.004950 |
| Rejection rate when $\alpha$ = 0.05 | 0.056 | 0.056 |
The bias of $ \hat{V}_U$ is still twice as high as $\hat{V}_B$, but the bias of either is so low that it does not matter, and both tests have a similar rejection rate.
This goes against the results of @whuber, who claims without justification that the variance of the z-test statistic $\hat{z}_{U}$ is 1/2, while in fact it is slightly above 1.
I define $c_{nk} = \frac{n_k}{n}$, $k = 1,2$. The true variance of $\hat{d}$ can be written as
$$\begin{aligned} Var(\hat{d}) &= \frac{1}{n} p(1 - p)\left( \frac{1}{c_{n1}} + \frac{1}{c_{n2}}\right) \end{aligned}$$
If we assume that the two samples grow roughly at the same speed, $c_{nk} \to c_k$, with $0 < c_k < 1$, then from the Central Limit Theorem
$$\frac{\hat{d}}{Var(\hat{d})} = n^{-1/2}\frac{\hat{d}}{\sqrt{ p(1 - p)\left( \frac{1}{c_{n1}} + \frac{1}{c_{n2}} \right)}} \to_d n^{-1/2}\frac{\hat{d}}{\sqrt{ p(1 - p)\left( \frac{1}{c_{1}} + \frac{1}{c_{2}} \right)}} \to_d N(0,1) $$
which justifies the use of the normal distribution.
However since $p$ is unknown, we can't use this to construct the test statistic. Instead, we need to use our sample estimates to get an asymptotically valid estimator.
According to the Strong Law of Large Numbers, we have
$$\begin{aligned} \hat{p}_1 \to_{a.s.} p \\ \hat{p}_2 \to_{a.s.} p \\ \hat{p} = \frac{x_1 + x_2}{n} \to_{a.s.} p \end{aligned}$$
From the Continuous Mapping Theorem, we have
$$\begin{aligned} \hat{p} (1 - \hat{p})\left(\frac{1}{c_{n1}} + \frac{1}{c_{n2}} \right)\to_{a.s.} p(1-p)\left(\frac{1}{c_1} + \frac{1}{c_2} \right) \\ \frac{\hat{p}_1(1 - \hat{p}_1)}{c_{n1}} + \frac{\hat{p}_2(1 - \hat{p}_2)}{c_{n2}} \to_{a.s.} p(1-p)\left(\frac{1}{c_1} + \frac{1}{c_2} \right) \end{aligned}$$
Therefore, from Slutsky's theorem, where can
$$ \begin{aligned} n^{-1/2}\frac{\hat{d}}{\sqrt{ \hat{p}(1 - \hat{p})\left( \frac{1}{c_{n1}} + \frac{1}{c_{n2}} \right)}} &\to_{d} n^{-1/2}\frac{\hat{d}}{\sqrt{ \frac{\hat{p}_1(1 - \hat{p}_1)}{c_{n1}} + \frac{\hat{p}_2(1 - \hat{p}_2)}{c_{n2}}}} \\ &\to_{d} n^{-1/2}\frac{\hat{d}}{\sqrt{ p(1 - p)\left( \frac{1}{c_{n1}} + \frac{1}{c_{n2}} \right)}} \\ &\to_{d} N(0,1) \end{aligned}$$
In order words, asymptotically and when the null hypothesis is true, there is no difference between using the exact variance, the unpooled variance, or the pooled variance.
Under the alternative hypothesis that $p_1 \ne p_2$, we have
$$ \hat{p} \to_{a.s} p_1 c_{1} + p_2 c_{2} = \tilde{p} $$
Therefore, we will asymptotically have, for the pooled estimator
$$ \hat{p} (1 - \hat{p}) \left(\frac{1}{c_{n1}} + \frac{1}{c_{n2}}\right) \to_{a.s.} \tilde{p}(1 - \tilde{p})(\frac{1}{c_1} + \frac{1}{c_2}) $$
Meanwhile, the unpooled estimator will keep converging to the actual variance, because
$$ \frac{ \hat{p}_1 (1 - \hat{p}_1)}{c_{n1}} + \frac{ \hat{p}_2 (1 - \hat{p}_2)}{c_{n2}} \to_{a.s.} = \frac{ p_1 (1 - p_1)}{c_{1}} + \frac{ p_2 (1 - p_2)}{c_{2}} $$
The unpooled estimator will be asymptotically more powerful when
$$\frac{1}{\sqrt{\tilde{p}(1 - \tilde{p})(\frac{1}{c_1} + \frac{1}{c_2})}} < \frac{1}{ \sqrt{ \frac{ p_1 (1 - p_1)}{c_{1}} + \frac{ p_2 (1 - p_2)}{c_{2}} } }$$
or equivalently
$$ \tilde{p}(1 - \tilde{p})(\frac{1}{c_1} + \frac{1}{c_2}) > \frac{ p_1 (1 - p_1)}{c_{1}} + \frac{ p_2 (1 - p_2)}{c_{2}} $$
Either estimator could be more powerful, depending on $p_1$, $p_2$, $c_1$ and $c_2$.
However, when the two sample sizes are equal ($c_1 = c_2 = 0.5$), the condition simplifies to
$$ \tilde{p}(1 - \tilde{p}) > \frac{p_1 (1 - p_1) + p_2 (1 - p_2)}{2} $$
Because $f(x) = x(1-x)$ is strictly concave, the previous inequality will be always true. Therefore, for equal sample sizes, the unpooled test is always asymptotically more powerful.
