r - 有没有办法绘制障碍模型结果 pscl 包或绘制计数 - r 中障碍模型的零截断 negbin 部分？-6ren

r - 有没有办法绘制障碍模型结果 pscl 包或绘制计数 - r 中障碍模型的零截断 negbin 部分？

转载作者：行者123 更新时间：2023-12-05 04:20:31

有没有办法在 R 中绘制障碍模型结果？我能够绘制障碍模型的零部分(带有 logit 链接的二项式)(如下)，但我无法弄清楚如何绘制模型的计数部分(带有对数链接的截断负二项式)。我正在为跨栏模型使用 pscl 包。

示例数据(df 名称 = 数据):

structure(list(Byths = c(333L, 107L, 0L, 0L, 684L, 0L, 113L, 
0L, 0L, 20L, 251L, 20L, 0L, 0L, 32L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 182L, 30L, 0L, 0L, 183L, 0L, 0L, 0L, 8L, 213L, 108L, 
310L, 960L, 720L, 0L, 0L, 6L, 72L, 78L, 15L, 196L, 256L, 608L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 105L, 0L, 194L, 0L, 440L, 0L, 0L, 
0L, 0L, 0L, 18L, 0L, 239L, 262L, 0L, 102L, 17L, 0L, 0L, 0L, 0L, 
0L, 68L, 93L, 0L, 226L, 118L, 91L, 330L, 104L, 68L, 224L, 0L, 
0L, 18L, 79L, 71L, 8L, 73L, 38L, 39L, 7L), Season = structure(c(3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("Pre", 
"Peak", "Post"), class = "factor"), depthm = c(20.1170446232626, 
9.1441111923921, 18.2882223847842, 18.2882223847842, 17.373811265545, 
13.7161667885881, 13.7161667885881, 13.106559375762, 13.106559375762, 
8.53450377956596, 8.53450377956596, 12.4969519629359, 20.1170446232626, 
20.1170446232626, 8.83930748597903, 19.2026335040234, 19.2026335040234, 
20.1170446232626, 20.1170446232626, 20.1170446232626, 20.1170446232626, 
20.1170446232626, 14.0209704950012, 14.0209704950012, 14.0209704950012, 
14.0209704950012, 14.0209704950012, 8.53450377956596, 12.8017556693489, 
20.7266520360888, 20.7266520360888, 17.373811265545, 14.0209704950012, 
14.0209704950012, 14.0209704950012, 8.22970007315289, 8.22970007315289, 
18.5930260911973, 17.373811265545, 8.53450377956596, 7.92489636673982, 
7.92489636673982, 7.92489636673982, 7.92489636673982, 7.92489636673982, 
7.92489636673982, 9.44891489880517, 9.44891489880517, 9.44891489880517, 
9.44891489880517, 15.8497927334796, 12.1921482565228, 8.83930748597903, 
8.83930748597903, 18.8978297976103, 8.83930748597903, 8.83930748597903, 
17.0690075591319, 20.4218483296757, 20.4218483296757, 20.1170446232626, 
20.1170446232626, 13.106559375762, 13.106559375762, 13.106559375762, 
13.106559375762, 13.106559375762, 13.106559375762, 8.83930748597903, 
10.6681297244574, 9.75371860521824, 9.75371860521824, 9.75371860521824, 
11.5825408436967, 9.1441111923921, 17.0690075591319, 17.0690075591319, 
17.0690075591319, 17.0690075591319, 17.0690075591319, 17.0690075591319, 
20.1170446232626, 10.0585223116313, 17.373811265545, 20.1170446232626, 
20.4218483296757, 20.4218483296757, 7.92489636673982, 7.92489636673982, 
7.92489636673982, 7.92489636673982, 7.92489636673982, 7.92489636673982, 
7.92489636673982, 7.92489636673982)), row.names = c(5L, 7L, 12L, 
13L, 14L, 16L, 18L, 21L, 22L, 23L, 24L, 26L, 28L, 29L, 32L, 33L, 
34L, 35L, 36L, 37L, 38L, 39L, 46L, 48L, 49L, 50L, 51L, 52L, 53L, 
54L, 55L, 56L, 57L, 59L, 61L, 62L, 63L, 65L, 66L, 71L, 73L, 75L, 
77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 92L, 95L, 98L, 
107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 
118L, 119L, 125L, 137L, 139L, 142L, 143L, 144L, 147L, 150L, 151L, 
152L, 154L, 156L, 157L, 159L, 160L, 161L, 162L, 163L, 165L, 177L, 
178L, 179L, 180L, 181L, 183L, 184L, 185L), class = "data.frame")

型号:

mod <- hurdle(Byths ~ depthm * Season, data = data, dist = "negbin")

summary(mod)

给予:

Call:
hurdle(formula = Byths ~ depthm * Season, data = data, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.8326 -0.5676 -0.2256  0.1534  4.7222 

Count model coefficients (truncated negbin with log link):
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        5.88538    0.63004   9.341  < 2e-16 ***
depthm            -0.05014    0.05577  -0.899 0.368658    
SeasonPeak        -3.03180    0.86315  -3.512 0.000444 ***
SeasonPost        -2.74270    1.95696  -1.402 0.161061    
depthm:SeasonPeak  0.17058    0.07106   2.400 0.016380 *  
depthm:SeasonPost  0.21029    0.13113   1.604 0.108783    
Log(theta)         0.13512    0.19772   0.683 0.494367    
Zero hurdle model coefficients (binomial with logit link):
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        4.02258    1.10922   3.627 0.000287 ***
depthm            -0.32024    0.08317  -3.851 0.000118 ***
SeasonPeak        -2.55477    1.65557  -1.543 0.122798    
SeasonPost        -3.94742    3.35163  -1.178 0.238892    
depthm:SeasonPeak  0.27367    0.12088   2.264 0.023569 *  
depthm:SeasonPost  0.30061    0.21990   1.367 0.171617    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.1447
Number of iterations in BFGS optimization: 15 
Log-likelihood: -344.6 on 13 Df

我对系数取幂来解释它们，但是我想为我的论文做一些数字。

我弄清楚如何绘制障碍模型的零部分的唯一方法是执行以下操作:

#turned counts into presence/absence

structure(list(Byths = c(1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 
1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 
1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1), 
    Season = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("Pre", 
    "Peak", "Post"), class = "factor"), depthm = c(20.1170446232626, 
    9.1441111923921, 18.2882223847842, 18.2882223847842, 17.373811265545, 
    13.7161667885881, 13.7161667885881, 13.106559375762, 13.106559375762, 
    8.53450377956596, 8.53450377956596, 12.4969519629359, 20.1170446232626, 
    20.1170446232626, 8.83930748597903, 19.2026335040234, 19.2026335040234, 
    20.1170446232626, 20.1170446232626, 20.1170446232626, 20.1170446232626, 
    20.1170446232626, 14.0209704950012, 14.0209704950012, 14.0209704950012, 
    14.0209704950012, 14.0209704950012, 8.53450377956596, 12.8017556693489, 
    20.7266520360888, 20.7266520360888, 17.373811265545, 14.0209704950012, 
    14.0209704950012, 14.0209704950012, 8.22970007315289, 8.22970007315289, 
    18.5930260911973, 17.373811265545, 8.53450377956596, 7.92489636673982, 
    7.92489636673982, 7.92489636673982, 7.92489636673982, 7.92489636673982, 
    7.92489636673982, 9.44891489880517, 9.44891489880517, 9.44891489880517, 
    9.44891489880517, 15.8497927334796, 12.1921482565228, 8.83930748597903, 
    8.83930748597903, 18.8978297976103, 8.83930748597903, 8.83930748597903, 
    17.0690075591319, 20.4218483296757, 20.4218483296757, 20.1170446232626, 
    20.1170446232626, 13.106559375762, 13.106559375762, 13.106559375762, 
    13.106559375762, 13.106559375762, 13.106559375762, 8.83930748597903, 
    10.6681297244574, 9.75371860521824, 9.75371860521824, 9.75371860521824, 
    11.5825408436967, 9.1441111923921, 17.0690075591319, 17.0690075591319, 
    17.0690075591319, 17.0690075591319, 17.0690075591319, 17.0690075591319, 
    20.1170446232626, 10.0585223116313, 17.373811265545, 20.1170446232626, 
    20.4218483296757, 20.4218483296757, 7.92489636673982, 7.92489636673982, 
    7.92489636673982, 7.92489636673982, 7.92489636673982, 7.92489636673982, 
    7.92489636673982, 7.92489636673982)), row.names = c(5L, 7L, 
12L, 13L, 14L, 16L, 18L, 21L, 22L, 23L, 24L, 26L, 28L, 29L, 32L, 
33L, 34L, 35L, 36L, 37L, 38L, 39L, 46L, 48L, 49L, 50L, 51L, 52L, 
53L, 54L, 55L, 56L, 57L, 59L, 61L, 62L, 63L, 65L, 66L, 71L, 73L, 
75L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 92L, 95L, 
98L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 
117L, 118L, 119L, 125L, 137L, 139L, 142L, 143L, 144L, 147L, 150L, 
151L, 152L, 154L, 156L, 157L, 159L, 160L, 161L, 162L, 163L, 165L, 
177L, 178L, 179L, 180L, 181L, 183L, 184L, 185L), class = "data.frame")

然后绘制模型的障碍部分...

ggplot(dataPA, aes(x = depthm, y = Byths)) + 
  geom_point() +
  stat_smooth(method = "glm", method.args = list(family="binomial"), se = TRUE) + facet_wrap(~ Season)

有没有办法绘制障碍模型(截断的 negbin)的计数模型部分？还是总体上绘制障碍模型结果的更好方法？

我应该为此使用 predict.hurdle 并使用预测值吗？如果是这样，我该怎么做，我发现的示例不够清晰，我无法理解，而且我已经尝试了几周但没有奏效。

谢谢!

最佳答案

绘制模型的计数部分

您可以轻松地完成与障碍部分相同的操作，但现在使用模型的计数部分。您只需要使用负二项式模型，但只需要使用 Byths > 0 的数据。另一个问题是 glm() 没有负二项式模型的功能，因此您需要使用 MASS::glm.nb()

给你:

ggplot(data[data$Byths != 0, ], aes(x = depthm, y = Byths)) + 
  geom_point() +
  stat_smooth(method = MASS::glm.nb, se = TRUE) + facet_wrap(~ Season)

现在我想确保这会得到与使用函数作者的 predict 功能相同的结果。好消息是我们这样做了，但在上面我们很容易得到我们预测的置信区间。

newdat <- expand.grid(
  depthm = seq(8, 20, .25),
  Season = c("Pre", "Peak", "Post")
)
newdat$count_mean <- predict(mod, newdat, type = "count")
ggplot(newdat, aes(x = depthm, y = count_mean)) + 
  geom_point() + 
  facet_wrap(~ Season) +
  ylim(c(0, 1250))

绘制整体模型

如果这是我的研究，我可能希望整体查看我的模型结果。所以在这里我绘制了每个值的预期计数——结合模型的计数和障碍部分，我们的整体预测。我还想看看整个预测分布，而不仅仅是均值。所以在这里，我绘制了整个分布(在选定值处)，以便模型预测更加清晰。为此，我使用了 {ggdist} 包(我试图直接在那里输入模型预测，但不知道该怎么做，所以我从预测分布中采样了 10,000 次).

给你:

library(ggdist)

mean_dat <- expand_grid(
  depthm = seq(8, 20, 0.25),
  Season = c("Pre", "Peak", "Post")
) %>% 
  mutate(
    mean_count = predict(mod, ., type = "response")
  ) 
  
simdat <- expand_grid(
  depthm = c(8, 12, 16, 20),
  Season = c("Pre", "Peak", "Post")
  ) %>% 
  mutate(
    counts = predict(mod, ., type = "prob", at = 0:1000)
  ) %>% 
  nest(counts) %>% 
  mutate(
    sim_count = map(data, ~ sample(0:1000, size = 10000, replace = TRUE, prob = t(.x))),
    .keep = "unused"
  ) %>% 
  unnest()


ggplot(mean_dat, aes(x = depthm, y = mean_count)) + 
  geom_line() + 
  geom_line() + 
  stat_slab(aes(y = sim_count), fill = "skyblue", alpha = .5, data = simdat) +
  facet_wrap(~ Season)

黑线显示预测变量值的预期计数，蓝色分布显示某些选定值的预测分布(样本来自)。

现在，在我看来，这种可视化在显示模型的实际预测方面做得更好。大多数值仍然集中在 0 附近，但存在很大差异，尤其是当 depthm 变大时。

关于r - 有没有办法绘制障碍模型结果 pscl 包或绘制计数 - r 中障碍模型的零截断 negbin 部分？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/74440752/

文章推荐： Flutter 在气泡聊天消息上对齐时间戳文本

文章推荐： python - 有没有办法欺骗类的 isinstance 结果？

文章推荐： python - FutureWarning in using iteritems() in use .iloc() pandas

c++ - 使用 SDL_Renderer 绘制 2D 内容，使用 SDL_GLContext 绘制 OpenGL 内容
我学习 SDL 二维编程已有一段时间了，现在我想创建一个结合使用 SDL 和 OpenGL 的程序。我是这样设置的: SDL_Init(SDL_INIT_VIDEO); window = SDL_Cr
绘制 map 投影类型
尝试查找可在地块中使用的不同类型项目的列表来自不同样本的投影类型: projection = list(type = "equirectangular") projection = list(typ
Java 绘制 GIF
我正在尝试使用 Java Graphics API 绘制 GIF，但无法使用下面的代码成功绘制 GIF。仅绘制 GIF 的第一张图像或缩略图，但不播放。 public void paintCompon
Java JFrame 绘制
我目前正在使用 JFrame 并尝试绘制一个矩形，但我不知道如何执行代码 paint(Graphics g)，如何获取 Graphics 对象？ package com.raggaer.frame;
java - 绘制 ImageView
这个领域的新手，希望得到一些帮助。我有一个"Missile.java" 类，我在那里画东西。我想绘制一个 ImageView，我正在使用以下代码: ImageView v = (ImageView)
HTML5 Canvas - 绘制
下面列出了圆形的例子这是我的 JavaScript 代码。最佳答案假设您的 randomColor 是正确的，您只需要: 从 canvas.onclick 中移除 context.clearR
Android在ImageView上缩放、拖动、绘制
我在绘制和缩放 ImageView 时遇到问题。请帮帮我.. 当我画一些东西然后拖动或缩放图像时 - 绘图保留在原处，如您在屏幕截图中所见。而且我只需要简单地在图片上绘图，并且可以缩放和拖动这张图片。
c# - 绘制/绘制外部形式
我们可以在形式之外绘制图像和文本...我的意思是在字面上... 我知道问这个问题很愚蠢但是我们能不能... 最佳答案您可以通过创建表单并将其 TransparentColor 属性设置为背景色来“作
java - 绘制/布局期间的对象分配？
我在绘制/布局期间收到 3 个对象分配警告 super.onDraw(canvas); canvas.drawColor(Color.WHITE); Paint textPaint = new Pai
python - 绘制 Pandas 时间序列数据框的线性回归线的置信区间
我有一个示例时间序列数据框: df = pd.DataFrame({'year':'1990','1991','1992','1993','1994','1995','1996',
r - 绘制 R 数据框中所有列的分布
我试图想出一种简洁的方法来绘制 R 数据框中所有列的 GridView 。问题是我的数据框中既有离散值又有数值。为简单起见，我们可以使用 R 提供的名为 iris 的示例数据集。我会使用 par(mf
r - 绘制 "list"的密度
我有一个由 10 列和 50 行组成的 data.frame。我使用 apply 函数逐列计算密度函数。现在我想绘制我一次计算的密度。换句话说，而不是绘图... plot(den[[1]]) plo
r - 绘制 PCA 的所有组件
我想知道我们如何才能在第一个和第二个组件之外绘制个人，如下所示: 最佳答案这可能有效: pc.cr <- princomp(USArrests, cor = TRUE) pairs(pc.cr$lo
pandas - 绘制 Pandas DataFrame时缺少xticklabels的第一个值
我是Pandas和matplotlib的新手，想绘制此DataFrame season won team matches pct_won 0 20
python - 绘制 Distplot 子图
我正在尝试为 distplot 子图编写一个 for 循环。我有一个包含许多不同长度列的数据框。 (不包括 NaN 值) fig = make_subplots( rows=len(asse
r - 绘制 3d 密度
我想创建一个具有密度的 3d 图。我使用函数 density 首先为特定的 x 值创建一个二维图，然后该函数创建密度并将它们放入 y 变量中。现在我有第二组 x 值并将其再次放入密度函数中，然后我得
python - 绘制 OpenStreetMap 关系不会生成连续线
全部，我一直在研究全局所有 MTB 步道的索引。我是 Python 人，所以对于所有涉及的步骤，我都尝试使用 Python 模块。我能够像这样从 OSM 立交桥 API 中获取关系: from O
r - 绘制 SVM 分类图时出错
我正在使用 e1071 包中的支持向量机对我的数据进行分类，并希望可视化机器实际如何进行分类。但是，在使用 plot.svm 函数时，出现无法解决的错误。脚本: library("e1071") d
r - 绘制 XTS 对象时的变化
我制作了以下图表，它是使用 xts 对象创建的。我使用的代码很简单 plot(graphTS1$CCLL, type = "l", las = 2, ylab = "(c)\nCC for I
uml - 绘制 UML 状态图
在绘制状态图时，您如何知道哪些状态放在框中，哪些状态用于转换箭头？我注意到转换也是状态。我正在查看 this page 上的图 1 : 最佳答案转换不是状态。转换是将对象从一种状态移动到下一种状态

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城