也许我不清楚，但是您设置的以下模拟似乎适用于您建议的随机效果结构，即
set.seed(2019)
N <- 50 # number of subjects
# id indicators for pairs
ids <- t(combn(N, 2))
id_i <- ids[, 1]
id_j <- ids[, 2]
# random effects
b_i <- rnorm(N, sd = 2)
b_j <- rnorm(N, sd = 4)
# simulate normal outcome data from the mixed model
# that has two random effects for i and j
y <- 10 + b_i[id_i] + b_j[id_j] + rnorm(nrow(ids), sd = 0.5)
DF <- data.frame(y = y, id_i = id_i, id_j = id_j)
library("lme4")
#> Loading required package: Matrix
lmer(y ~ 1 + (1 | id_i) + (1 | id_j), data = DF)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ 1 + (1 | id_i) + (1 | id_j)
#> Data: DF
#> REML criterion at convergence: 2403.071
#> Random effects:
#> Groups Name Std.Dev.
#> id_i (Intercept) 2.1070
#> id_j (Intercept) 3.0956
#> Residual 0.5053
#> Number of obs: 1225, groups: id_i, 49; id_j, 49
#> Fixed Effects:
#> (Intercept)
#> 9.562
但是，我认为您将无法在主题级别包含协变量，因为您只有配对级别的数据。
编辑：根据评论，数据的对称性变得更加清晰。据我所知，当前的实现lmer()不允许这样的数据。下面的代码使用 STAN 模拟和拟合此类数据的模型。
set.seed(2019)
N <- 50 # number of subjects
# id indicators for pairs
ids <- expand.grid(i = seq_len(N), j = seq_len(N))
ids <- ids[ids$i != ids$j, ]
id_i <- ids$i
id_j <- ids$j
# random effects
b <- rnorm(N, sd = 2)
# simulate normal outcome data from the mixed model
# that has one random effect but accounts for the pairs i and j
y <- 10 + b[id_i] + b[id_j] + rnorm(nrow(ids), sd = 0.5)
library("rstan")
Data <- list(N = nrow(DF), n = length(unique(id_i)),
id_i = id_i, id_j = id_j, y = y)
model <- "
data {
int n;
int N;
int id_i[N];
int id_j[N];
vector[N] y;
}
parameters {
vector[n] b;
real beta;
real<lower = 0> sigma_b;
real<lower = 0> sigma;
}
transformed parameters {
vector[N] eta;
for (k in 1:N) {
eta[k] = beta + b[id_i[k]] + b[id_j[k]];
}
}
model {
sigma_b ~ student_t(3, 0, 10);
for (i in 1:n) {
b[i] ~ normal(0.0, sigma_b);
}
beta ~ normal(0.0, 10);
sigma ~ student_t(3, 0, 10);
y ~ normal(eta, sigma);
}
"
fit <- stan(model_code = model, data = Data, pars = c("beta", "sigma_b", "sigma"))
summary(fit)