GLM и GEEGLM работают только с меньшим / конкретным набором данных

Question

К сожалению, я не могу предложить полностью воспроизводимый пример, потому что я не могу поделиться данными. Тем не менее, я надеюсь, что кто-то может помочь мне понять следующее.

Данные
В моем наборе данных 134 столбца и 2521 ряд. Для анализа хочу выполнить ГИ (geepack::geeglm) но проблема также возникает в простом glm, Интересующие нас столбцы в модели:

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   2521 obs. of  6 variables:
 $ SUBJID         : chr  "01" "01" "01" "01" ...
 $ util_trans     : num  0 0 0 0 0 0 0 0.431 0.225 0.139 ...
 $ base_utility   : num  0 0 0 0 0 0 0 0.431 0.431 0.431 ...
 $ trt_01         : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 1 1 1 ...
 $ priorreg_factor: Factor w/ 2 levels "1",">1": 1 1 1 1 1 1 1 1 1 1 ...
 $ avisit_group   : Factor w/ 4 levels "baseline","treatment",..: 1 2 2 2 2 3 3 1 2 2 ...

Модель
Я предоставляю код для glm() так как это очень похоже на geepack::geeglm(),

Подгонка модели следующим образом возвращает ошибку:

glm(util_trans ~ I(base_utility) +
             factor(trt_01) +
             factor(priorreg_factor),
     data = na.omit(db),
     subset = avisit_group == "treatment")
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

Однако, если я выберу только необходимые столбцы, модель будет работать отлично;

glm(util_trans ~ I(base_utility) +
             factor(trt_01) +
             factor(priorreg_factor),
     data = na.omit(db %>% dplyr::select(SUBJID, util_trans, base_utility,
                                          trt_01, priorreg_factor, avisit_group)),
     subset = avisit_group == "treatment")

Call:  glm(formula = util_trans ~ I(base_utility) + factor(trt_01) + 
    factor(priorreg_factor), data = na.omit(db.eq5.seq %>% dplyr::select(SUBJID, 
    util_trans, base_utility, trt_01, priorreg_factor, avisit_group)), 
    subset = avisit_group == "treatment")

Coefficients:
              (Intercept)            I(base_utility)            factor(trt_01)1  factor(priorreg_factor)>1  
                0.09                  0.1                  0.02                  0.2  

Degrees of Freedom: 1118 Total (i.e. Null);  1115 Residual
Null Deviance:      32.89 
Residual Deviance: 22.47    AIC: -1187

Имейте в виду, что я вручную изменил значение коэффициентов, чтобы "анонимизировать" их.

Почему разница в результате, хотя данные и вызов функции остаются неизменными?

Answer 1

Вы вызываете na.omit на весь фрейм данных. Это приведет к тому, что строки будут опущены, если они содержат какие-либо значения NA. Кажется, это приводит к тому, что отбрасывается достаточное количество строк, и у вас остается только один уровень, по крайней мере, для одного из ваших факторов. Вот пример этого в уменьшенном масштабе, чтобы проиллюстрировать

> dat <- data.frame(x = factor(c(1,1,1,2)), y = 1:4, unrelated = c(2,5,3,NA))
> dat
  x y unrelated
1 1 1         2
2 1 2         5
3 1 3         3
4 2 4        NA
> na.omit(dat)
  x y unrelated
1 1 1         2
2 1 2         5
3 1 3         3

> na.omit(dat[,c("x", "y")])
  x y
1 1 1
2 1 2
3 1 3
4 2 4

Обратите внимание, что когда мы использовали na.omit, включая несвязанную переменную, он удалял единственную строку, у которой был уровень "2" для x. Если мы явно выберем столбцы, мы заботимся о них, чтобы сохранить эту строку в данных.