У меня есть сомнения относительно пакета R "поля". Я оцениваю логистическую модель:
modelo1 <- glm(VD ~ VE12 + VE.cont + VE12:VE.cont + VC1 + VC2 + VC3 + VC4, family="binomial", data=data)
Где:
VD2 - дихотомическая переменная (1 болезнь / 0, а не болезнь)
VE12 - это дихотомическая переменная воздействия (со значениями 0 и 1)
VE.cont переменная непрерывного воздействия
VCx (остальные переменные) являются смешанными переменными.
Моя цель - получить прогнозируемую вероятность заболевания (VD2) для вектора значений VE.cont и для каждой группы VE12, но с учетом VCx переменных. Другими словами, я хотел бы получить линию доза-ответ между VD2 и VE.cont группой VE12, но предполагая одинаковое распределение VCx для каждой линии доза-ответ (т.е. без смешения).
Следуя номенклатуре этой статьи (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4052139/) Я думаю, что мне следует выполнить «предельную стандартизацию» (метод 1), которую можно выполнить с помощью stata, но я не уверен, как я могу это сделать с Р.
Я использую этот синтаксис (с R):
cdat0 <- cplot(modelo1, x="VE.cont", what="prediction", data = data[data[["VE12"]] == 0,], draw=T, ylim=c(0,0.3))
cdat1 <- cplot(modelo1, x="VE.cont", what="prediction", data = data[data[["VE12"]] == 1,], draw=marg"add", col="blue")
но я не уверен, правильно ли я это делаю, потому что этот подход дает результаты, аналогичные использованию модели без смешения переменных и функции predict.glm.
modelo0 <- glm(VD2 ~ VE12 + VE.cont + VE12:VE.cont, family="binomial", data=data)
Возможно, мне следует использовать параметр поля, но я не понимаю результатов, потому что значения, полученные в столбце VE.cont, не находятся в шкале вероятности (от 0 до 1).
x <- c(1,2,3,4,5)
margins::margins(modelo1, at=list("VE.cont"=x, "VE12"=c(0,1)), type="response")
Это пример рисунка, который я хотел бы получить: