4

我创建了一个零膨胀负二项式模型，并想调查有多少零被划分为采样或结构零。我如何在 R 中实现这一点。我不清楚 zeroinfl 页面上的示例代码。

data("bioChemists", package = "pscl")

fm_zinb2 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin")

table(round(predict(fm_zinb2, type="zero"))) 
>   0   1 
> 891  24 

table(round(bioChemists$art))
    >   0   1   2   3   4   5   6   7   8   9  10  11  12  16  19 
    > 275 246 178  84  67  27  17  12   1   2   1   1   2   1   1

这是在告诉我什么？

当我对我的数据执行相同操作时，我得到的读数只是在 1? 谢谢

4

1 回答 1

11

详细信息在 Zeileis (2008) 的论文中，网址为https://www.jstatsoft.org/article/view/v027i08/v27i08.pdf

predict收集关于函数对库中每个模型的作用的所有解释需要做一些工作（几年后，您的问题仍未得到解答）pscl，并将其隐藏在数学表达式中（第 19,23 页）似然函数（方程 7、8）。我将您的问题解释为您想要/需要知道如何使用不同type的预测：

  • 预期计数是多少？( type="response")
  • 超零的（条件）预期概率是多少？( type="zero")
  • 任何计数的（边际）预期概率是多少？( type="prob")
  • 最后，有多少预测零是多余的（例如抽样）而不是基于回归的（即结构性的）？

要读入 pscl 包附带的数据：

data("bioChemists", package = "pscl")

然后拟合一个零膨胀负二项式模型：

fm_zinb2 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin")

如果您希望预测预期值，则使用

predict(fm_zinb2, type="response")[29:31]
       29        30        31 
0.5213736 1.7774268 0.5136430

因此，在这个模型下，29 岁和 31 岁的生物化学家在博士的最后 3 年发表的文章的预期数量是一半，而 30 岁的生物化学家则接近 2 篇。

但是我相信您追求的是过零的概率（在零点质量中）。此命令执行此操作并打印第 29 到 31 行中项目的值（是的，我去钓鱼了！）：

predict(fm_zinb2, type="zero")[29:31]

它产生这个输出：

        29         30         31 
0.58120120 0.01182628 0.58761308

因此，第 29 项是多余零的概率（您将其称为抽样零，即非结构性零，因此没有由协变量解释）为 58%，第 30 项为 1.1%，第 31 项是 59%。因此，这两位生物化学家预计发表的文章为零，这超出了可以通过各种协变量的负二项式回归来解释的那些。

您已经将整个数据集的这些预测概率制成表格

table(round(predict(fm_zinb2, type="zero"))) 
  0   1 
891  24

因此，您的输出告诉您，只有 24 位生物化学家可能是超零，即预测的超零概率超过二分之一（由于四舍五入）。

如果您将百分比刻度上的 10 个点制成表格，可能会更容易解释

table(cut(predict(fm_zinb2, type="zero"), breaks=seq(from=0,to=1,by=0.1)))

给予

 (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6] 
     751        73        34        23        10        22 
(0.6,0.7] (0.7,0.8] (0.8,0.9]   (0.9,1] 
        2         0         0         0

因此，您可以看到 751 位生化学家不太可能成为超零，但 22 位生化学家有 50-60% 的机会成为超零，只有 2 位有更高的机会（60-70%）。没有人极有可能成为一个过分的零。以图形方式，这可以显示在直方图中

hist(predict(fm_zinb2, type="zero"), col="slateblue", breaks=seq(0,0.7,by=.02))

您将每个生化学家的实际计数制成表格（无需四舍五入，因为这些是计数）：

table(bioChemists$art)
  0   1   2   3   4   5   6   7   8   9  10  11  12  16  19 
275 246 178  84  67  27  17  12   1   2   1   1   2   1   1

发表19篇论文的特约生化学家是谁？

most_pubs <- max(bioChemists$art)
most_pubs
extreme_biochemist <- bioChemists$art==most_pubs
which(extreme_biochemist)

您可以获得每个生化学家拥有任意数量的酒吧的估计概率，恰好为 0 并且达到最大值，这里是难以置信的 19！

preds <- predict(fm_zinb2, type="prob")
preds[extreme_biochemist,]

你可以为我们的一位特殊的生物化学家看这个，他有 19 篇出版物（这里使用 base R 绘图，但 ggplot 更漂亮）

expected <- predict(fm_zinb2, type="response")[extreme_biochemist]
# barplot returns the midpoints for counts 0 up to 19
midpoints<-barplot(preds[extreme_biochemist,], 
  xlab="Predicted #pubs", ylab="Relative chance among biochemists")
# add 1 because the first count is 0
abline(v=midpoints[19+1],col="red",lwd=3)
abline(v=midpoints[round(expected)+1],col="yellow",lwd=3)

这表明，尽管我们预计 biochemist 915 有 4.73 篇出版物，但在此模型下，2-3 家酒吧的可能性更大，远不及实际的 19 家酒吧（红线）。

#pubs 获得 biochemist profile 29 的机会

回到问题，对于 biochemist 29，超零的概率是

pzero <- predict(fm_zinb2, type="zero")
pzero[29]
       29 
0.5812012

总体（边际）为零的概率是

preds[29,1]
[1] 0.7320871

所以预测的零概率与结构（即由回归解释）的比例是：

pzero[29]/preds[29,1]
       29 
0.7938962

或者，在超出零的可能性之外，零的附加概率是：

preds[29,1] - pzero[29]

       29 
0.1508859

biochemist 29 的实际出版物数量是

bioChemists$art[29]
[1] 0

因此，预测 biochemist 发表文章为零的原因很少能通过回归 (20%) 得到解释，而且大多数情况下没有解释（即过量，80%）。

总的来说，我们看到对于大多数生物化学家来说，情况并非如此。我们的生物化学家 29 是不寻常的，因为他们零酒吧的机会大多是多余的，即回归无法解释。我们可以通过以下方式看到：

hist(pzero/preds[,1], col="blue", xlab="Proportion of predicted probability of zero that is excess")

这给了你：

超过生化学家的零预测概率的比例

于 2017-10-18T20:58:48.817 回答