Графики с планками погрешностей с использованием объекта logistf

Первоначально я создал логистическую модель с помощью пакета glm, но хотел исправить разделение, поэтому я использовал функцию logistf и теперь пытаюсь переделать свои графики. Я не уверен, как сделать график, подобный приведенному ниже, с объектом logistf. Многие пакеты, похоже, не поддерживают его, я пробовал использовать функцию plot_model() пакетов sjPlot, которая рисует точку для прогнозируемой вероятности, но не добавляет планки погрешностей, как это делается автоматически с объектом glm. Как я могу обойти это? Возможно, есть другой пакет, который облегчит эту задачу, или есть способ вручную добавить планки ошибок?

Код для графика, к которому я хочу добавить планки ошибок:

sjPlot::plot_model(lr3, type="int", mdrt.values = "meansd", show.values = TRUE, value.offset = .3)

Результат моей модели LR3:

logistf(formula = foodbank_cv ~ wave + ff_country + relevel(race_grp, 
    ref = "White") + sex_cv + age_r + relevel(numchildren, 
    ref = "None") + wave * ff_hcondhas + relevel(carer, 
    ref = "Not") + sempderived + wave * cd_ff_furlough + 
    log(ff_hours) + qual + num + relevel(keyworksector, ref = "Not keyworker") + 
    ca_clinvuln_dv + freemeals + ca_blbenefits1 + log(hhincome_week), 
    data = data, firth = TRUE, family = binomial(link = "logit"))

Model fitted by Penalized ML
Coefficients:
                                                                         coef   se(coef)  lower 0.95  upper 0.95        Chisq            p method
(Intercept)                                                      -5.237542354 0.46736532 -6.23016284 -4.30807241          Inf 0.000000e+00      2
wave5                                                            -0.377956413 0.32598420 -1.07410577  0.28545651 1.232122e+00 2.669947e-01      2
wave7                                                            -0.929934987 0.40813067 -1.84652632 -0.12926473 5.260388e+00 2.181615e-02      2
ff_country2                                                      -0.118780142 0.33317501 -0.86893024  0.51197342 1.196576e-01 7.294061e-01      2
ff_country3                                                       0.393456771 0.25097814 -0.15010616  0.88210537 2.077828e+00 1.494527e-01      2
ff_country4                                                      -0.219066153 0.43493435 -1.23008781  0.57774984 2.481153e-01 6.184053e-01      2
relevel(race_grp, ref = "White")Asian or Asian British            0.882833792 0.22906054  0.39628625  1.33641305 1.183859e+01 5.801581e-04      2
relevel(race_grp, ref = "White")Black or Black British            1.759374627 0.27942672  1.16321835  2.29702048 2.678592e+01 2.272869e-07      2
relevel(race_grp, ref = "White")Mixed                             1.786978145 0.27773294  1.19285979  2.32350705 2.763841e+01 1.462461e-07      2
relevel(race_grp, ref = "White")Other                            -0.345106379 1.38712570 -5.19048868  1.62733736 6.509258e-02 7.986208e-01      2
ff_hcondhas                                                       0.691244774 0.26776923  0.14697164  1.25269746 6.228205e+00 1.257311e-02      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Код, который я использовал для создания сюжета о голоде и гонках. Я немного отредактировал вручную, чтобы он выглядел лучше, но в идеале я хочу, чтобы мой график выглядел так:

plot_model(model12, type = "pred", terms = c("race_grp"), mdrt.values = "meansd", axis.textsize = .3, wrap.labels = 5)+ theme_sjplot2() + scale_color_sjplot("simply") + ggplot2::labs(title= "Predicted probabilities of Hunger", x= "Race", y="Percentage")

введите здесь описание изображения

введите здесь описание изображения


person Farida El Kafrawy    schedule 24.05.2021    source источник
comment
Привет, можете ли вы предоставить некоторые образцы данных + код, который вы использовали для создания изображений?   -  person slamballais    schedule 24.05.2021
comment
@slamballais Привет, спасибо за ваш комментарий, я обновил сейчас   -  person Farida El Kafrawy    schedule 24.05.2021


Ответы (1)


Однако я нашел способ обойти эту проблему, но не с помощью пакета logistf. Если кто-то в будущем захочет узнать ответ на этот вопрос, я предлагаю вам использовать пакет brglm. Я проверил, и результаты пакета brglm точно такие же, как и пакета logistf. Вот как я воссоздал сюжет Голода, опубликованный выше:

hi2<- brglm(formula= hungry_cv~ wave + ff_country + race_grp + sex_cv + age_r + numchildren + wave*ff_hcondhas + carer + sempderived + wave*cd_ff_furlough + log(ff_hours) + qual + num + keyworksector + ca_clinvuln_dv + freemeals + ca_blbenefits1 + log(hhincome_week), data=data, family=binomial(logit), method = "brglm.fit", pl = TRUE)
racehunger<- plot_model(hi2, type = "pred", terms = c("race_grp"), mdrt.values = "meansd", axis.textsize = .3, wrap.labels = 5, show.values = TRUE)+ theme_sjplot2() + ggplot2::labs(title= "Predicted probabilities of Hunger", x= "Race", y="Percentage")
racehunger
png(file="racehunger.png", units="in", width=11, height=8.5, res=300)
print(racehunger)
dev.off()

Вывод кода:

введите здесь описание изображения

Я лично очень доволен результатом.

person Farida El Kafrawy    schedule 26.05.2021
comment
Кроме того, попробуйте brglm2 (обновление, хотя я сомневаюсь, что результаты будут отличаться) - person Ben Bolker; 27.05.2021