Also, I should probably include full simulation and model code.

First, the simulation:

```
## Simulating data for multiple dyadic responses
library(mvtnorm)
library(MASS)
nodes <- data.frame(row.names = 1:300)
nodes$i_ID <- seq(1:300)
nodes$j_ID <- nodes$i_ID
nodes$k_ID <- rep(1:15, each = 20)
d.full <- expand.grid(nodes$i_ID, nodes$j_ID)
colnames(d.full) <- c("i_ID", "j_ID")
d.nodes <- subset (d.full, i_ID != j_ID)
d.nodes$k_ID <- nodes$k_ID [ match (d.nodes$i_ID, nodes$i_ID)]
d.nodes$l_ID <- nodes$k_ID [ match (d.nodes$j_ID, nodes$j_ID)]
d <- subset (d.nodes, k_ID == l_ID)
d <- subset (d, select = -(l_ID))
d$first.node <- ifelse( d$i_ID < d$j_ID, paste (d$i_ID), paste(d$j_ID))
d$second.node <- ifelse( d$i_ID > d$j_ID, paste (d$i_ID), paste(d$j_ID))
d$ij_ID <- paste(d$first.node, "_", d$second.node, sep = "")
d <- subset (d, select = -c(first.node, second.node))
d$ij_dyad_ID <- match (d$ij_ID, unique(d$ij_ID))
d <- subset (d, select = -(ij_ID))
d <- subset (d, i_ID < j_ID)
d <- d[order(d$i_ID),]
row.names(d) <- d$ij_dyad_ID
## Prepping dyad-level random effects -- that is, the residuals
N <- nrow(d)
SIGMA <- matrix (NA, nrow = 4, ncol = 4)
# covariance matrix is defined by:
# [ s1 c12 r1 r12 ]
# [ s2 r12 r2 ]
# [ s1 c12 ]
# [ s2 ]
s1 <- 1.5
s2 <- .5
c12 <- s1 * s2 * -0.5
r1 <- s1 * s1 * 0.5
r2 <- s2 * s2 * 0.5
r12 <- s1 * s2 * -0.5
SIGMA[1,] <- c(s1, c12, r1, r12)
SIGMA[2,] <- c(c12, s2, r12, r2)
SIGMA[3,] <- c(r1, r12, s1, c12)
SIGMA[4,] <- c(r12, r2, c12, s2)
v_e <- rmvnorm(N, mean = c(0,0,0,0), sigma = SIGMA)
## Set intercept
a1 <- 5
a2 <- 6
## Make i to j response, outcome 1
Y1 <- rep (NA, N)
for ( i in 1:N ) {
Y1[i] <- with ( d ,
{ a1 +
v_e [ ij_dyad_ID[i] , 1 ]
}
)
}
## Make i to j response, outcome 2
Y2 <- rep (NA, N)
for ( i in 1:N ) {
Y2[i] <- with ( d ,
{ a2 +
v_e [ ij_dyad_ID[i] , 2 ]
}
)
}
## Make j to i response, outcome 1
Y3 <- rep (NA, N)
for ( i in 1:N ) {
Y3[i] <- with ( d ,
{ a1 +
v_e [ ij_dyad_ID[i] , 3 ]
}
)
}
## Make j to i response, outcome 1
Y4 <- rep (NA, N)
for ( i in 1:N ) {
Y4[i] <- with ( d ,
{ a2 +
v_e [ ij_dyad_ID[i] , 4 ]
}
)
}
d$Y1 <- Y1
d$Y2 <- Y2
d$Y3 <- Y3
d$Y4 <- Y4
dat_list <- list (
N = N,
Y1 = d$Y1,
Y2 = d$Y2,
Y3 = d$Y3,
Y4 = d$Y4,
nresp = 4,
ij_dyad_ID = d$ij_dyad_ID
)
```

And the current model code is:

```
model_code <- "
data {
int<lower=1> N; // total number of observations
real Y1[N]; // response variable
real Y2[N]; // response variable
real Y3[N]; // response variable
real Y4[N]; // response variable
int<lower=1> nresp; // number of responses
}
transformed data {
vector[nresp] Y[N]; // response matrix
for (n in 1:N) {
Y[n] = [Y1[n], Y2[n], Y3[n], Y4[n]]';
}
}
parameters {
real response_1_Intercept;
real response_2_Intercept;
real<lower=0> sigma_response_1; // residual SD for response 1
real<lower=0> sigma_response_2; // residual SD for response 2
// Covariances
real<lower=-1,upper=1> rho_c12;
real<lower=-1,upper=1> rho_r1;
real<lower=-1,upper=1> rho_r2;
real<lower=-1,upper=1> rho_r12;
}
transformed parameters {
// residual-level effects and covariances
matrix [4,4] RHO;
matrix [4,4] Lrescor;
// Transform residuals from residual matrix
RHO[1,1] = sigma_response_1;
RHO[1,2] = rho_c12;
RHO[1,3] = rho_r1;
RHO[1,4] = rho_r12;
RHO[2,1] = rho_c12;
RHO[2,2] = sigma_response_2;
RHO[2,3] = rho_r12;
RHO[2,4] = rho_r2;
RHO[3,1] = rho_r1;
RHO[3,2] = rho_r12;
RHO[3,3] = sigma_response_1;
RHO[3,4] = rho_c12;
RHO[4,1] = rho_r12;
RHO[4,2] = rho_r2;
RHO[4,3] = rho_c12;
RHO[4,4] = sigma_response_2;
Lrescor = cholesky_decompose(RHO);
}
model {
vector[N] mu_Y1 = rep_vector(0, N) + response_1_Intercept;
vector[N] mu_Y2 = rep_vector(0, N) + response_2_Intercept;
vector[N] mu_Y3 = rep_vector(0, N) + response_1_Intercept;
vector[N] mu_Y4 = rep_vector(0, N) + response_2_Intercept;
// Defining Mu and sigma vectors
vector[nresp] Mu[N];
vector[nresp] sigma = [sigma_response_1, sigma_response_2, sigma_response_1, sigma_response_2]';
// cholesky factor of residual covariance matrix
matrix[nresp, nresp] LSigma = diag_pre_multiply(sigma, Lrescor);
for (n in 1:N) {
Mu[n] = [mu_Y1[n], mu_Y2[n], mu_Y3[n], mu_Y4[n]]';
}
// priors including all constants
target += student_t_lpdf(response_1_Intercept | 3, 5, 10);
target += student_t_lpdf(sigma_response_1 | 3, 0, 10)
- 1 * student_t_lccdf(0 | 3, 0, 10);
target += student_t_lpdf(response_2_Intercept | 3, 6, 10);
target += student_t_lpdf(sigma_response_2 | 3, 0, 10)
- 1 * student_t_lccdf(0 | 3, 0, 10);
target += lkj_corr_cholesky_lpdf(Lrescor | 1);
// likelihood including all constants
{
target += multi_normal_cholesky_lpdf(Y | Mu, LSigma);
}
}
"
```