R-非贝叶斯方法中的混合效应多项式 Logistic 回归模型？

Question

我正在尝试在 R 中运行至少有 1 个（但最好是 2 个）随机效应的多项逻辑回归，但非常不成功。似乎很少有软件包可以处理多项式模型，和具有随机效应。

我已经（匆忙）匿名化了我的数据，所以假设我正在调查 80 年代不同比赛中运动员的网球搭档的分布情况。每行 (n=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

我想运行这个模型：

对手 ~ 种族 + 排名 + Opp_rank + Rel_rank + Opp_rel_rank + 胜利 + opp_wins +(1|运动员) + (1|日期)

关键数据集特征/障碍： 1.我的分类响应变量是一个有 14 个级别的因子 2.运动员有时会在同一天打多场比赛，因此运动员是伪复制的来源 3. 日期也是伪复制的一个来源，因为同一天的比赛都具有相同的排名和胜利值。 4. 大多数软件包似乎只允许 1 个随机效果

任何帮助或见解将不胜感激。甚至告诉我我想做什么也是不可能的！ 如果需要的话，我也愿意将我的响应变量更改为数字变量，只要它能提供有关对手的信息。例如，我可以使用对手的等级或年龄。使用明确的回应似乎是我许多问题的核心

我已经尝试了所有我能找到的相关建模包。 multinom、mlogit 和 mclogit 软件包都无法处理随机效应。包 glmmTMB 和 glmmADMB 都无法处理多项模型。 glmer、lme4或任何其他类似的众所周知的软件包也不能。 *编辑; mclogit 对我来说不起作用，因为数据重组对于我的数据集来说太复杂了。实际上有 100 列，因此扩展它对我的计算机来说太费力了

似乎满足我要求的唯一传统方法是来自 mixcat 包的 npmlt ，但我每次都会遇到相同的错误，这似乎与我的随机效应有关。我尝试删除所有 NA 值，仅使用一种随机效应，将随机效应合并为一个，但没有任何效果。在下面的代码中，“comb”是我在测试过程中组合随机效应的列。我读到在运行这个包之前需要附加数据。

> 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()

我害怕的另一个似乎确实合适的选择是使用 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)

我尝试的最后一个贝叶斯选项是MCMCglmm，但这需要复杂的数据重组，而我没有能力完成。

Answer 1

数据库下载

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

将连续变量集中在平均值和标准差上

data$Rank2 <- (data$Rank - mean(data$Rank))/sd(data$Rank)
data$Rel_rank2 <- (data$Rel_rank - mean(data$Rel_rank))/sd(data$Rel_rank)
data$Opp_rank2 <- (data$Opp_rank - mean(data$Opp_rank))/sd(data$Opp_rank)
data$Opp_rel_rank2 <- (data$Opp_rel_rank - mean(data$Opp_rel_rank))/sd(data$Opp_rel_rank)
data$wins2 <- (data$wins - mean(data$wins))/sd(data$wins)
data$opp_wins2 <- (data$opp_wins - mean(data$opp_wins))/sd(data$opp_wins)

将分类变量转换为因素

data$Date2 <- as.factor(data$Date)
data$Athlete2 <- as.factor(data$Athlete)
data$Opponent2 <- as.factor(data$Opponent)
data$Ethnicity2 <- as.factor(data$Ethnicity)

混合模型方程！！一些变量彼此共线，因此提供相同的信息。我只包含了一些变量，因为样本太小，无法包含模型中的所有变量。

library(mclogit)

modele <- mblogit(formula = Opponent2 ~ Ethnicity2 + Rank2 + Opp_rank2, random = list(~1|Athlete2, ~ 1|Date2), data = data, method=c("PQL"), estimator=c("ML"))

获取混合模型参数

coefs <- summary(modele)$coefficients
LLs <- coefs[,1] + qnorm(.025)*coefs[,2]
ULs <- coefs[,1] + qnorm(.975)*coefs[,2]
OR <- exp(coefs[,1])
ORLL <- exp(LLs)
ORUL <- exp(ULs)
HHES <- coefs[,1]/1.81 # Hasselblad and Hedges Effect Size

获取混合模型的系数

round(cbind(coefs, LLs, ULs), 3)

获取混合模型的优势比

round(cbind(OR, ORLL, ORUL, HHES), 3)