I have a very general question about whether it is possible (and/or sensible) to specify correlations within a model between random effects from distinct clusters. That is, if I have random effects that do not belong to the same “population” (e.g. male- and female-specific random effects with distinct (co)variance), can I specify correlations across those clusters (e.g. do males’ male-specific random effects correlate with their partner’s female-specific random effects). Clearly this can be done by simply correlating posterior estimates between the groups, but I wonder if there is a more principled way to estimate this correlation parameter in the model, despite the fact that there are already separate correlation matrices for the effects within each group. I have not yet provided code because I believe this may be a general conceptual issue. If not, I will gladly share more code for my specific situation, which is described in detail below.

Currently, I have data on male and female subjects observed in multiple dyads over time. Assume that it is appropriate to use a bivariate mixed model to estimate separate models for each sex. (In my case, I have distinct spatial auto-correlation matrices for each sex that are informing their distinct random effect (co)variances.) I’d like to know whether some random effect values correlate between male and female partners, e.g. I want to know if male random intercepts and slopes tend to positively correlate with their female partners’ random intercepts and slopes.

One approach (ignoring the spatial auto-correlation) would be to simply correlate dyad random effects between the models, assuming dyads are always constant. This would effectively correlate the effects for male and female responses given that “dyad” always indexes the same pair. Suppressing non-essential features, this would look something like the following for the random effects

Male \hspace{0.1cm} response_i = (\beta_0 + \mu_{0_{dyad{M[i]}}}) + (\beta_1 + \mu_{1_{dyadM[i]}}) + \varepsilon _i

Female \hspace{0.1cm} response_i = (\beta_0 + \mu_{0_{dyad{F[i]}}}) + (\beta_1 + \mu_{1_{dyadF[i]}}) + \varepsilon _i

\begin{bmatrix}
\mu_{0_{dyad{M}}} \\
\mu_{1_{dyad{M}}} \\
\mu_{0_{dyad{F}}} \\
\mu_{1_{dyad{F}}}
\end{bmatrix} \sim MVN
\begin{bmatrix}
\sigma^2_{\mu_{0_{dyad{M}}}} & ... & ... & \sigma_{\mu_{0_{dyad{M}}}, {\mu_{1_{dyad{F}} }} } \\
\vdots & \sigma^2_{\mu_{1_{dyad{M}}}} & ... & \vdots \\
\vdots & \vdots & \sigma^2_{\mu_{0_{dyad{F}}}} & \vdots \\
\sigma_{\mu_{1_{dyad{F}}}, {\mu_{0_{dyad{M}} }} } & .... & .... & \sigma^2_{\mu_{1_{dyad{F}}}}
\end{bmatrix}

This approach seems to breaks down, however, if dyads change across observations, as it would not account for the fact that the same individuals are observed across multiple dyads, where it instead becomes necessary to include separate random effects for individual and dyad, which then have distinct meanings.

Male \hspace{0.1cm} response_i = (\beta_0 + \mu_{0_{subject{M[i]}}}+ \mu_{0_{dyad{M[i]}}}) + (\beta_1 + \mu_{1_{subject{M[i]}}} + \mu_{1_{dyadM[i]}}) + \varepsilon _i

Female \hspace{0.1cm} response_i = (\beta_0 + \mu_{0_{subject{F[i]}}}+ \mu_{0_{dyad{F[i]}}}) + (\beta_1 + \mu_{1_{subject{F[i]}}} + \mu_{1_{dyadF[i]}}) + \varepsilon _i

The index *i* can no longer be used to easily estimate the correlation(s) of interest because each male and female subject has multiple partners, so the correlation of male *i* and female *i* is meaningless. Also note that the question is whether pairs tend to have correlated individual random effects (not whether the dyad-specific effects correlate after partialling out individual-level effects).

I have tried modelling this by including an additional correlation matrix in my model, so that there are male and female-specific random effect correlations estimated as well as a separate set of correlations between paired males and females. This tends to result in bizarre and clearly incorrect estimates, which seems to be due to the fact that the correlations are being “pulled” in some sense between the separate matrices. If, however, I simply correlate the posterior estimates of random effects between the two models, when the male and female values are treated as effectively independent during estimation, I recover the expected associations from simulation. Is this a case where it would be wiser to just correlate posterior random effect estimates, or is there a smarter way for me to tell the model that I expect these effects to correlate across models in addition to the random effect correlations already specified within each model?

Thanks in advance for your time.