I am trying to implement a generalized mixed-effects model specified as:
Dependent variable $y = \log(\frac{p}{1 - p})$ where $p$ is a quantity measured for a pair of individuals ($i$ and $j$).
$E(y_{ij} \mid X_{ij}, a_i, b_j) = \beta X_{ij} + a_i + b_j,$
where $X_{ij}$ represents my covariates (predictor or independent variables) for individuals $i$ and $j$, $\beta$ represents my fixed effects regression coefficients, and $a_i$ and $b_j$ represent the random effects (crossed effects) for the two individuals involved (every data point represents some feature values or measurements for a pair of individuals and we want to control for the idiosyncratic contributions of the individuals involved in a pair by conditioning on these individuals via the random crossed effects).
$ (a_i, b_i) \sim MVN(0, \Sigma), \quad \Sigma = \ \left[ {\begin{array}{cc} \sigma_a^2 & \rho\sigma_a\sigma_b \\ \rho\sigma_a\sigma_b & \sigma_b^2 \\ \end{array} } \right]. $
This model is fit using the Laplace approximation, and the authors of the paper whose model I am following use the lme4 package in R (I am open to both R and a way to do this in Python).
I am unable to understand how to implement these crossed effects for the individuals with the generalized mixed-effects model. I have my $y$ and $X$ in place and have the IDs for each pair and each individual making up a particular pair. I feel like I need to use the individual IDs for my random effects here, but I am not sure how. Most default examples seem to add a random intercept or random slopes corresponding to the covariates specified in $X$ but this does not seem to be the case here.
If anyone has any pointers about how to implement the above-specified model using lmer, it would help a lot. I hope I specified the model clearly enough and would be happy to provide an application scenario or point to the paper which describes and uses this model if that helps.
lme4. $\endgroup$