lme4 - предупреждение glmer: «Downdated VtV не является положительно определенным» - это потому, что переменная отклика слишком однородна - PullRequest
Новая
       66

lme4 - предупреждение glmer: «Downdated VtV не является положительно определенным» - это потому, что переменная отклика слишком однородна

0 голосов
/ 09 сентября 2018

У меня есть данные, которые выглядят как моделирование ниже.В данных представлены 3 разные группы (g1, g2 и g3) субъектов, которым потребовались 1, 2 или 3 процедуры (p1, p2, p3) соответственно в зависимости от тяжести заболевания, варьирующейся между субъектами.Следовательно, существует только p1 для группы g1;p1, p2 для группы g2;p1, p2, p3 для группы g3.Каждое измерение было выполнено во время каждой процедуры.То, что я пытаюсь сделать, - это прогон, чтобы соответствовать смешанной модели логистической регрессии (извините, если это не правильная номенклатура: я рад узнать, если это неправильно), чтобы искать разницу в доле пациентов на p1, которые получилиопределенное лечение (y1 дихотомическая переменная) и контролировать идентификатор пациента, устанавливая его в качестве случайной величины.

С этими смоделированными данными вывод

mod_glmer = glmer(y1 ~ procedure + (1|id), data = data_g2, family = binomial)

дает мне следующую ошибку:

Error in length(value <- as.numeric(value)) == 1L : Downdated VtV is not positive definite

Я предполагаю, что это связано с тем фактом, что на p1 есть только 1 при выполнении следующего:

table(data_g2$y1[data_g2$procedure == "p1"])

, тогда как оба 0 и 1 в p2 при выполнении следующего:

table(data_g2$y1[data_g2$procedure == "p2"])

Однако я также выполнил аналогичную "стандартную" логистическую регрессию и не получаюсообщение об ошибке:

mod_glm = glm(y1 ~ procedure, data = data_g2, family = binomial)

Что я должен сделать, чтобы продемонстрировать, что доля y1 отличается в p1 по сравнению с p2 при использовании glmer?

PS: если вы смоделируете данные, запустите несколько раз, так как это включает генерацию случайного числа (я получаю сообщение один раз каждые 2 запуска, если случайно, каждый раз выдает ошибку с set.seed(1). Попробуйте другойчисло n в set.seed(n), если вы не получите сообщение об ошибке)

Моделирование данных 

library(lubridate)
library(dplyr)
library(lme4)

set.seed(1)

# define the number of subjects in each groups
n_g3 = 20
n_g2 = n_g3 * 10
n_g1 = n_g3 * 10

# 3 different groups
id_p1 = paste0("ID",1:c(n_g1 + n_g2 + n_g3))
id_p2 = paste0("ID",c(n_g1+1):c(n_g1 + n_g2 + n_g3))
id_p3 = paste0("ID",c(n_g1 + n_g2+1):c(n_g1 + n_g2 + n_g3))
id = append(append(id_p1, id_p2), id_p3)

# 3 different groups
groups = c(rep("g1", n_g1), rep("g2", n_g2), rep("g3", n_g3), rep("g2", n_g2), rep("g3", n_g3), rep("g3", n_g3))

# 1st, 2nd or 3rd procedure
procedure = c(rep("p1", n_g1), rep("p1", n_g2), rep("p1", n_g3), rep("p2", n_g2), rep("p2", n_g3), rep("p3", n_g3))

# measurement n??2 dichotomous
y1 = c(rep(1, n_g1 + n_g2 + n_g3),
  (rbinom(c(n_g2 + n_g3),1,0.8)),
  (rbinom(n_g3,1,0.3))
)

data = data.frame(id, groups, procedure, y1)

### create subset g2 of data to make comparisons within the group with 2 procedures
data_g2 = data[which(data$groups == "g2"),]

table(data_g2$y1[data_g2$procedure == "p1"])
table(data_g2$y1[data_g2$procedure == "p2"])

mod_glm = glm(y1 ~ procedure, data = data_g2, family = binomial)
summary(mod)

mod_glmer = glmer(y1 ~ procedure + (1|id), data = data_g2, family = binomial)
summary(mod_glmer)
...