У меня возникли проблемы с построением графика предельных эффектов моей отрицательной биномиальной регрессии с нулевой инфляцией, особенно для модели с нулевой инфляцией.
Я хотел бы создать график, который визуализирует предельный эффект моей переменной в модели с нулевой инфляцией, как показано ниже.
Моя модель определяется как:
model_1 <- zeroinfl(soaGEN ~
lguerrillad +
lstrikesd +
ldemond +
lriotsd +
ldmidtoward | lfpsimusa,
dist = "negbin",
data_model)
summary(model_1)
Call:
zeroinfl(formula = soaGEN ~ lguerrillad + lstrikesd + ldemond + lriotsd + ldmidtoward |
lfpsimusa, data = data_model, dist = "negbin")
Pearson residuals:
Min 1Q Median 3Q Max
-0.67403 -0.49186 -0.41840 0.06559 9.58756
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.93654 0.06873 57.273 < 2e-16 ***
lguerrillad 0.87534 0.10735 8.154 3.52e-16 ***
lstrikesd -0.21427 0.11339 -1.890 0.0588 .
ldemond -0.01439 0.10429 -0.138 0.8903
lriotsd -0.07146 0.10447 -0.684 0.4939
ldmidtoward 0.16672 0.10886 1.532 0.1256
Log(theta) -0.55302 0.06590 -8.392 < 2e-16 ***
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3032 0.1683 7.745 9.59e-15 ***
lfpsimusa -0.6360 0.0630 -10.096 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 0.5752
Number of iterations in BFGS optimization: 19
Log-likelihood: -5365 on 9 Df
Однако я не могу изолировать часть модели с нулевой инфляцией от остальной части модели при попытке визуализировать.
sjPlot::plot_model(model_1,type = "pred",terms = "lfpsimusa")
Можно ли каким-либо образом изолировать или извлечь модель с нулевой инфляцией из полной модели, а затем построить ее или какой-либо другой способ получить желаемый результат?
Вы можете использовать пакет marginaleffects для вычисления и построения различных количественных показателей для этих моделей. (Отказ от ответственности: я являюсь сопровождающим.)
Ссылка выше включает виньетку «Приступая к работе», виньетку сюжета и более 20 других виньеток, а также подробную документацию. Графическая виньетка, в частности, показывает, как настраивать графики.
library(marginaleffects)
library(pscl)
data("bioChemists", package = "pscl")
mod <- zeroinfl(art ~ phd * ment + factor(fem) + factor(mar) | 1, data = bioChemists)
plot_predictions(mod, condition = "ment")
plot_slopes(mod, variables = "ment", condition = "phd")
'Большое спасибо! с добавленным type = "zero"
он работает отлично!
Вы можете сделать это, используя ggpredict()
из пакета ggeffects
напрямую:
library(pscl)
#> Classes and Methods for R developed in the
#> Political Science Computational Laboratory
#> Department of Political Science
#> Stanford University
#> Simon Jackman
#> hurdle and zeroinfl functions by Achim Zeileis
library(ggplot2)
data("bioChemists", package = "pscl")
fm_zinb <- zeroinfl(art ~ fem + mar + kid5 + phd | ment, data = bioChemists, dist = "negbin")
library(ggeffects)
g <- ggpredict(fm_zinb, terms = "ment [all]", type = "zi_prob")
plot(g) +
geom_rug(data=bioChemists, aes(x=ment), inherit.aes = FALSE)
Created on 2023-04-05 with reprex v2.0.2
Большое спасибо! Это решение отлично сработало для моей проблемы! Мне было интересно, можно ли добавить к этому графику доверительные интервалы?
По какой-то причине доверительный интервал не создается. Вы можете использовать решение @Vincent выше, попробуйте: marginaleffects::plot_predictions(fm_zinb, condition = "ment", type = "zero") + geom_rug(data=bioChemists, aes(x=ment), inherit.aes = FALSE)
для модели в моем ответе.
Спасибо, что предоставили эти подробности @DaveArmstrong! Обратите внимание, что в последней версии пакета в rug
есть аргумент plot_predictions()
.
Использование
type = "zero"
в функцииplot_predictions()
создаст график вероятности нулевой инфляции.