如何解决R语言中mblogit的变量长度不一致错误

问题重现 #

你的数据框定义如下：

df <- data.frame(File=c("F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8", "F9", "F10"),
                 Type=c("A", "B", "C", "A", "A", "C", "C", "B", "B", "A"),
                 Time=c("T1", "T2", "T1", "T3", "T1", "T2", "T2", "T3", "T1", "T3"),
                 One_count=c(15, 8, 9, 20, 11, 0, 6, 19, 22, 7),
                 Two_count=c(3, 12, 1, 6, 10, 7, 11, 5, 2, 9),
                 Three_count=c(12, 15, 10, 25, 28, 19, 16, 21, 17, 24),
                 Four_count=c(0, 7, 18, 11, 5, 2, 15, 21, 7, 9),
                 Five_count=c(3, 6, 2, 1, 9, 12, 8, 17, 10, 0))

尝试运行以下多项回归代码时出现错误：

value_count <- cbind('Three_count', 'One_count', 'Two_count', 'Four_count', 'Five_count')
mblogit(formula = value_count ~ Time + Type + Time*Type, data=df, random=~1|File, method="MQL", estimator = "ML", dispersion = T)

错误信息：

Error in model.frame.default(formula = value_count ~ Time + Type + Time *  : 
  variable lengths differ (found for 'Time')

错误原因 #

1. 变量类型错误：cbind('Three_count', ...)生成的是长度为5的字符向量，而数据框中Time等变量的长度是10，两者长度不匹配，导致model.frame无法构建模型矩阵。
2. 数据格式不符合要求：多项Logistic回归要求响应变量是表示类别的因子变量，而你的数据是宽格式（每个类别对应一列计数），需要转换为长格式（每行对应一个类别-观测组合，包含计数信息）才能适配模型要求。

解决方案 #

步骤1：转换数据为长格式 #

使用tidyr包将宽格式数据转换为长格式，提取类别名称并指定Three为参考类别：

library(tidyr)

# 转换为长格式
df_long <- df %>%
  pivot_longer(
    cols = c(One_count, Two_count, Three_count, Four_count, Five_count),
    names_to = "Value",
    values_to = "Count"
  ) %>%
  # 去除类别名称中的"_count"后缀
  mutate(
    Value = gsub("_count", "", Value),
    # 将Value转为因子，指定Three为参考类别
    Value = factor(Value, levels = c("Three", "One", "Two", "Four", "Five"))
  )

步骤2：正确构建多项Logistic回归模型 #

根据是否需要随机效应，提供两种常用实现方案：

场景1：带随机效应的混合模型（glmmTMB）

如果需要保留File作为随机分组变量，使用glmmTMB包实现混合效应多项Logistic回归：

library(glmmTMB)

model <- glmmTMB(
  Value ~ Time + Type + Time*Type + (1|File),
  data = df_long,
  weights = Count,  # 用Count作为频数权重
  family = multinomial(),
  contrasts = list(Value = "contr.treatment")
)

summary(model)

场景2：无随机效应的多项Logistic回归（nnet）

如果不需要随机效应，使用nnet包的multinom函数：

library(nnet)

model <- multinom(
  Value ~ Time + Type + Time*Type,
  data = df_long,
  weights = Count
)

summary(model)

注意事项 #

• 若坚持使用mblogit函数（如brglm2包中的实现），需确保数据为长格式后，将公式调整为「因子型响应变量 ~ 预测变量」的形式，并正确指定权重与随机效应参数。
• 当数据存在零计数或过度离散情况时，可根据分布特征调整模型族（如负二项分布）。

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