Social relations model for multiple groups

Question

I want to apply the social relations model (SRM) to data from an experiment in which participants were exchanging resources with one another. A pair of actors i and j is then a primary unit of analysis and I am interested in the probability that each actor shares their resources with the other in the pair over a series of rounds (or trials). The specific implementation of SRM that I'd like to apply to my data is an adaptation of R code developed by Richard McElreath in his book "Statistical Rethinking" (see Chap. 14, Sect. 14.4: Social relations as correlated varying effect). In that chapter, McElreath illustrates the use of SRM on data on food sharing in Nicaragua studied by Koster and Leckie (2014).

Some relevant details about my problem are given below:

1. There were multiple experimental sessions. The experimental protocol, instructions, and tasks were identical across all sessions, but some elements of the experimental design differed across the sessions. These differences are not relevant for the model specification.

2. There were always 10 participants per session. Each participant could pick any of the other 9 participants as exchange partner and share resources with them. This situation is repeated 20 times. Thus, within each pair, actors have 20 opportunities to share resources with one another.

3. Given that a pair is a primary unit of analysis in my study, I want to fit and estimate a social relations model in which actors are nested within sessions. I want to be able to decompose the overall variation in resource sharing into the components associated with actor, partner, dyad, and session. Also, I want to be able to estimate two reciprocity parameters: (a) generalised reciprocity, or the actor-partner covariance of the same individual and (b) pair covariance, or the reciprocity in resource sharing within a given pair. The data I use form a multilevel structure in which resource sharing is nested within pairs and pairs are nested within experimental sessions.

I started by selecting only 1 experimental session and applying the following code to the selected data:

library(rethinking)

## R code 14.31
my_data <- list(
    N = nrow(my_dyads),
    N_actors = 10,
    did = my_data$did,
    hidA = my_data$subject,
    hidB = my_data$alter,
    giveAB = my_data$giveAB,
    giveBA = my_data$giveBA
)

my_model <- ulam(
    alist(
        giveAB ~ dbinom(20, pAB),
        giveBA ~ dbinom(20, pBA),
        logit(pAB) <- a + gr[hidA, 1] + gr[hidB, 2] + d[did, 1] ,
        logit(pBA) <- a + gr[hidB, 1] + gr[hidA, 2] + d[did, 2] ,
        a ~ normal(0, 1),

       ## gr matrix of varying effects
        vector[2]:gr[N_actors] ~ multi_normal(0, Rho_gr, sigma_gr),
        Rho_gr ~ lkj_corr(4),
        sigma_gr ~ exponential(1),

       ## dyad effects
        transpars> matrix[N, 2]:d <-
                compose_noncentered(rep_vector(sigma_d, 2), L_Rho_d, z),
        matrix[2,N]:z ~ normal(0, 1),
        cholesky_factor_corr[2]:L_Rho_d ~ lkj_corr_cholesky(8),
        sigma_d ~ exponential(1),

       ## compute correlation matrix for dyads
        gq> matrix[2,2]:Rho_d <<- Chol_to_Corr(L_Rho_d)
    ), 
    data = my_data, chains = 4, cores = 4, iter = 2000)

I am now looking for a way to extend the code, so that I could take the multilevel character of the data into account.