I am fitting a negative binomial-2 regression model where there is a multivariate normal random effects term. I would like to find an equation for the covariance of two outcomes. In "the multivariate Poisson-log normal distribution", Atchison and Ho found a solution for the Poisson distribution, but I would like to extend this to the negative binomial and I cannot find any previous work on this.
$$ Y_{i} \sim \text{negative-binomial}(\mu_{i} = \exp(X_i^T \beta + \phi_i), \theta) \text{ for } i = 1..n\\ % \log \left[\mathbb{E}(Y_{i}) \right] = X_{i}^T \beta^{\text{c}}_i + \phi_{i}\\ \phi \sim \mathcal{N}(0, \Sigma) \\ \Sigma \in \mathbb{R}^{n x n} % \phi_{.,t}|\phi_{.,t-1} \sim \mathcal{N}(\alpha^{\text{c}} \phi_{.,t-1}, \tau^2 \tilde{Q}^{-1}) \\ $$
What is $Cov(Y_i, Y_j)$?
