使用R语言glm函数时遇'variable lengths differ'错误的模型选择问题

问题场景 #

我在进行模型选择时，尝试用cbind(no,yes)构建二项式GLM模型，运行时出现"variable lengths differ"（变量长度不一致）错误。此前仅使用LD变量拟合回归模型时未出现问题，找不到修改代码的方法解决该问题。

我的代码 #

blh <- read.table('Data.dat', header=T)

blh$Level <- factor(blh$Level, labels=c('0','1','2'))
blh$Male <- factor(blh$Male, labels=c('female','male'))
blh$Night <- factor(blh$Night, labels=c('day','night'))
blh$Wet <- factor(blh$Wet, labels=c('dry','wet'))
blh$LD <- factor(blh$LD, labels=c('no','yes'))
ftable(xtabs(Count~Level+Male+Night+Wet+LD, data=blh))

fit <- glm(cbind(no,yes)~ Level+Male+Night+Wet, data=blh,family=binomial)

数据交叉表结构 #

LD no yes
Level Male   Night Wet          
0     female day   dry     3  12
                   wet     2  13
             night dry     2  13
                   wet     3  12
      male   day   dry    12  30
                   wet     9  33
             night dry    11  31
                   wet    15  27
1     female day   dry    12  36
                   wet     4  20
             night dry    11  37
                   wet     1  23
      male   day   dry    27  39
                   wet    10  23
             night dry    30  36
                   wet    10  23
2     female day   dry     9  27
                   wet     5  13
             night dry     8  28
                   wet     2  16
      male   day   dry    24  60
                   wet    15  27
             night dry    30  54
                   wet    13  29

错误原因 #

核心问题是：原始数据框blh中不存在no和yes这两个列。你看到的no/yes是xtabs生成的交叉表中LD的两个水平对应的计数，但这些值并没有被存储到blh里。原始数据的行数和交叉表中no/yes的长度不匹配，导致cbind(no,yes)找不到对应长度的变量，触发错误。

之前用LD拟合模型正常，是因为LD是原始数据中的列，配合Count权重就能正确映射分组观测数。

两种解决方法 #

方法1：将交叉表转换为带no/yes列的数据框 #

先把xtabs的结果转换成宽格式数据框，提取no和yes作为列，再拟合模型：

# 将交叉表转换为长格式数据框
tab_df <- as.data.frame(xtabs(Count~Level+Male+Night+Wet+LD, data=blh))
# 重塑为宽格式，拆分LD的no/yes为单独列
library(reshape2)
wide_df <- dcast(tab_df, Level+Male+Night+Wet ~ LD, value.var = "Freq")
# 拟合模型
fit <- glm(cbind(no, yes) ~ Level+Male+Night+Wet, data=wide_df, family=binomial)

方法2：直接使用原始数据+权重参数 #

不需要生成no/yes列，直接用原始数据中的LD作为响应变量，Count作为观测权重，更高效：

fit <- glm(LD ~ Level+Male+Night+Wet, data=blh, weights=Count, family=binomial)

这种方法和用cbind(no,yes)的二项式模型是等价的，都是拟合分组二项式回归。

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