I have a logistic mixed-effects model below.
$\mathbf{y}=(y_1, y_2, ..., y_n)^T$ is a n dimensions vector.
$\mathbf{p}=(p_1,p_2, ...,p_n)^T$ is a n dimensions vector.
$logit(\mathbf{p}) = (logit(p_1),logit(p_2), ...,logit(p_n))^T \sim N(0, \sigma_1^2\Sigma_1+\sigma_2^2\Sigma_2)$.
$\Sigma_1$ and $\Sigma_2$ are known n by n matrix. $\sigma_1^2$ and $\sigma_2^2$ are unknown variance components.
$y_i \sim Bernoulli(p_i)$
Now I have an observed $\mathbf{y}$. I would like to estimate the n by n covariance matrix $var(\mathbf{y})$.
From wikipedia https://en.wikipedia.org/wiki/Estimation_of_covariance_matrices, I can easily estimate the n by n covariance matrix $var(\mathbf{y})$ if I have several $\mathbf{y}$. But in my case, I only have an observed $\mathbf{y}$. I am just wondering if there is an easy way to estimate $var(\mathbf{y})$? I just need the mathematical formula.