Logistic Regression - Covariance Matrix of “Response” Residual Vector

Suppose that I have a logistic regression model

$logit(\pi_i) = \beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}, \quad i=1,\ldots,N$

Where:

There are $N$ "explanatory variable patterns" (possible combinations of the covariates),

$\mathbf{X}$ is the $N$ x $p$ design matrix,

At each "explanatory variable pattern", $i$, $y_i \sim Bin(n_i,\pi_i) \quad i=1,\ldots, N$

$\mathbf{y} = (y_1,\ldots,y_N)$ is the vector of binomial observed values,

$\mathbf{\pi} = (\pi_1,\ldots,\pi_N)$ is the vector of true probabilities,

$\mathbf{\hat{p}} = (\hat{p}_1,\ldots,\hat{p}_n)$ is the vector of fitted proportions, based on the logistic regression model, and

$\mathbf{N} = diag(n_1,\ldots,n_N)$ (the number of observations at each explanatory variable pattern/combination). Then, the "fitted values" can be referred to as $\mathbf{N \hat{p}}$.

Then, what I am interested in is $$Var(\mathbf{y} - N \mathbf{\hat{p}}) = Var(\mathbf{y}) + N \cdot Var(\mathbf{\hat{p}}) \cdot N' - Cov(\mathbf{y},\mathbf{\hat{p}})\cdot N'$$ That is, I am interested in the covariance matrix of the vector of "response" residuals from a logistic regression model. The name "response" comes from the similar name that is given when extracting residuals from a GLM in R.

For example, it is pretty easy to get $$\widehat{Var}(\mathbf{y}) = diag(n_1\hat{p}_1(1-\hat{p}_1),\ldots,n_N\hat{p}_N(1-\hat{p}_N))$$ Since the observations are binomial and independent.

Also, $Var(\mathbf{\hat{p}})$ is pretty straightforward, since $\widehat{Var}(\mathbf{\hat{\beta}}) = (X'VX)^{-1}$, (from the Fisher Information), where $V = \widehat{Var}(\mathbf{y})$, and then I can use the Delta Method to obtain an approximation to $Var(\mathbf{\hat{p}})$.

Really, what I can't get at is $Cov(\mathbf{y},\mathbf{\hat{p}})$, to finish the problem.

Maybe there is another way around all of this with a clever trick in another direction?