在R中实现含随机效应的多项Logistic回归（非贝叶斯方法）

我正尝试在R中运行至少包含1个（理想为2个）随机效应的多项Logistic回归，但始终未能成功——目前支持同时处理多项模型和随机效应的R包极少。

我将数据匿名化，模拟80年代网球运动员各场比赛的对手分布：数据集共1050行，每行代表一场比赛，包含比赛日期、目标运动员姓名、对手姓名、运动员种族、联赛排名、相对排名、对手排名、对手相对排名、当月总胜场数、对手当月总胜场数。

数据样例如下：

> head(df,5)
        Date Athlete Opponent      Ethnicity Rank Rel_rank Opp_rank Opp_rel_rank wins opp_wins
1 1987-09-10    Emma     Nora          Asian    2     0.10        8         0.83   12        7
2 1982-09-14  Olivia     Zoey      Caucasian    5     0.50        3         0.30    5        6
3 1988-11-21     Ava     Mila      Caucasian    9     0.93        9         0.93   11        9
4 1988-09-10  Sophia     Mila South American    8     0.83        7         0.83    7       10
5 1989-01-30  Amelia   Gianna      Caucasian    1     0.30        6         0.66   11       12

目标模型公式：

Opponent ~ Ethnicity + Rank + Opp_rank + Rel_rank + Opp_rel_rank + wins + opp_wins +(1|Athlete) + (1|Date)

数据集关键特征/难点 #

• 因变量Opponent为分类变量，是包含14个水平的因子
• 部分运动员单日参赛多场，Athlete存在伪重复问题
• Date存在伪重复，同日比赛的排名和胜场数值一致
• 多数R包仅支持单个随机效应

我已尝试所有能找到的相关建模方法，均遇阻碍：

1. 常规统计包 #

• multinom、mlogit、mclogit：不支持随机效应；其中mclogit需复杂的数据重组，实际数据集有100列，扩展后计算机负担过重，无法使用
• glmmTMB、glmmADMB、glmer、lme4：不支持多项模型

2. mixcat包（传统方法） #

唯一看似符合要求的是mixcat包中的npmlt函数，但运行始终出现随机效应相关错误。已尝试删除NA值、仅用单个随机效应、合并随机效应等操作，均无效。

运行代码及报错：

> all <- na.omit(all)
> all$Athlete <- as.factor(all$Athlete)
> all$Ethnicity <- as.factor(all$Ethnicity)

> attach(all)
The following objects are masked from all (pos = 3):

    Rel_rank, Opp_rank, Date, Wins, Ethnicity, Opponent, Opp_wins, comb, Rel_rank,
    Opp_rel_rank, Athlete

> model.po <- npmlt(formula = Opponent ~ Ethnicity + Rel_rank*Rank + Opp_rel_rank*Opp_rank + Wins + Opp_wins, 
+       formula.npo = ~ 1,       random = ~ 1 | Athlete,    k = 15)
Error in model.matrix.default(random, data = a) : 
  model frame and formula mismatch in model.matrix()

3. brms贝叶斯方法 #

可行但我不太想用：默认参数运行耗时一整晚，有效样本量（ESS）过低；已调整树深度和迭代次数，但尚未重新运行。此外我对贝叶斯模型验证信心不足，也不清楚如何为此类数据设置先验。

运行代码：

model <- brm(Opponent ~ Ethnicity + Opp_rank + Rel_rank 
             + Opp_rel_rank + wins + opp_wins + (1 | Athlete) + (1|Date),
             data = all, family = categorical(), 
             control = list(adapt_delta = 0.9, max_treedepth = 15),
             iter = 4000, chains = 4)

4. MCMCglmm贝叶斯方法 #

需复杂的数据重组，暂无能力完成。

我希望得到任何帮助或见解，哪怕告知此需求无法实现也可以；我也愿意将因变量替换为能有效反映对手特征的数值型变量（比如对手排名、年龄等），分类因变量似乎是问题核心。

内容的提问来源于stack exchange，提问作者MrCitywide