使用glm模型绘制两个不同组别的拟合线并考虑交互作用

3

我有以下的模型:

mod <- glm(data=data, events ~ treatment * size, family = quasipoisson)

以下是输出结果:

Call:
glm(formula = events ~ treatment * size, 
family = quasipoisson, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4842  -1.4939  -0.4199   0.5921   4.0068  

Coefficients:
                                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         0.070077   0.202376   0.346   0.7299    
treatmenttreatment                 -0.042710   0.315499  -0.135   0.8926    
size                                0.009464   0.002061   4.591   1.36e-05 *** 
treatmenttreatment:size             0.010270   0.005145   1.996   0.0488   *  
    ---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 2.039369)

    Null deviance: 252.78  on 97  degrees of freedom
Residual deviance: 191.43  on 94  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

我想使用R基础绘图功能绘制数据图,对于每个组(治疗和对照),使用不同颜色的点,并使用不同颜色的线拟合每个组。类似第20页这里的图表是理想的,但我不知道如何将该示例中使用的Poisson模型转换为我这里的准Poisson模型。下面包含可重复的示例。
individual  treatment  size     events
1           control    10.97    0
2           treatment  0        0
3           control    10.31    1
4           treatment  17.77    1
5           control    11.37    0
6           control    3.857    1
7           treatment  3.8      0
8           treatment  2.029    0
9           treatment  0.8571   0
10          control    0        0
11          control    0        0
12          control    7.943    0
13          control    0        0
14          treatment  8.514    0
15          control    0        0
16          treatment  28.69    1
17          treatment  39.03    4
18          treatment  33.49    0
19          control    2.514    0
20          control    2.771    1
21          treatment  3.257    0
22          control    24.6     1
23          control    1.714    1
24          treatment  9.343    1
25          treatment  10.86    2
26          treatment  28.77    3
27          control    89.97    6
28          control    17.17    0
29          control    4.057    0
30          control    20.4     2
31          treatment  28.49    3
32          treatment  28.66    1
33          treatment  30.66    1
34          control    8.114    0
35          treatment  29.03    2
36          treatment  0        0
37          control    6.543    0
38          treatment  18.86    1
39          control    42.37    3
40          treatment  9.257    3
41          treatment  29       3
42          control    13.46    0
43          control    8.143    0
44          control    0.08571  0
45          treatment  5.2      0
46          control    17.23    0
47          control    17.23    0
48          control    18.97    0
49          treatment  18.4     6
50          treatment  104.6    3
51          control    23.29    3
52          control    3.486    3
53          control    28.2     2
54          control    23       0
55          treatment  37.4     2
56          treatment  16.2     0
57          control    16.03    3
58          treatment  0        0
59          treatment  57.8     6
60          treatment  68.37    5
61          control    4.229    0
62          treatment  45.14    9
63          treatment  33.54    1
64          treatment  55.71    0
65          treatment  12.86    1
66          control    2.429    0
67          treatment  0        0
68          treatment  23.31    4
69          treatment  6.229    2
70          control    21.57    3
71          control    46.11    3
72          treatment  60.29    3
73          control    42.63    2
74          control    61.37    2
75          control    26.8     0
76          treatment  37.57    3
77          treatment  57.83    9
78          control    2.229    0
79          treatment  18.14    1
80          control    19.89    0
81          treatment  35.74    2
82          control    243.6    6
83          control    8.314    0
84          treatment  31.97    1
85          control    84.2     5
86          control    15.91    4
87          treatment  94.66    4
88          treatment  6.429    0
89          treatment  36.2     3
90          control    32.23    6
91          treatment  36.09    3
92          control    43.94    9
93          control    20.86    1
94          control    59.86    4
95          control    7.086    2
96          treatment  3.257    1
97          treatment  18.85    0
98          treatment  25.43    2

@Glen_b 在这里的回答 https://stats.stackexchange.com/questions/177921/plotting-the-results-of-glm-in-r/177926#177926 看起来很有帮助,但我还没有成功地让它工作。 - JKO
几天前,我展示了如何在https://stats.stackexchange.com/a/373259/919上进行三方交互。 - whuber
1个回答

2
我将使用基本的R图形而不进行太多解释来发布一种可能的解决方案。核心思想是使用predict在响应(原始)比例上生成预测值,然后绘制它们。其他所有内容基本上都是外观美化。个人而言,我会使用出色的visreg包轻松生成这种类型的图形。使用visreg的相应代码发布在本答案底部。

Quasipoisson_picture

mod <- glm(data=dat, events ~ treatment * size, family = quasipoisson)

plot(
  events~size
  , xlim = range(dat$size)
  , ylim = range(dat$events)
  , pch = 1
  , col = "#008FD0"
  , las = 1
  , data = subset(dat, treatment %in% "control")
)

points(
  events~size
  , pch = 1
  , col = "#F07E00"
  , data = subset(dat, treatment %in% "treatment")
)

legend(
  "topleft"
  , legend = c("treatment", "control")
  , pch = c(1, 1)
  , lwd = c(2, 2)
  , col = c("#F07E00", "#008FD0")
  , bty = "n"  
)

xx <- seq(min(dat$size), max(dat$size), length.out = 1000)

pred_frame <- expand.grid(
  size = xx
  , treatment = c("control", "treatment")
)

pred_frame$preds <- predict(mod, newdata = pred_frame, type = "response")

lines(
  preds~size
  , col = "#F07E00"
  , lwd = 2
  , data = subset(pred_frame, treatment %in% "treatment")
)

lines(
  preds~size
  , col = "#008FD0"
  , lwd = 2
  , data = subset(pred_frame, treatment %in% "control")
)

使用 visreg
visreg(
  mod
  , "size"
  , by = "treatment"
  , overlay = TRUE
  , scale = "response"
  , band = FALSE
  , ylim = range(dat$events)
)

@COOSerdash 这个脚本一般来说看起来表现不错,但是它似乎把点绘制在错误的位置。例如,对于另一个数据集,同样的函数会把几个点放在 y > 10 的位置,尽管所有的 y 都小于 10。这是为什么?我该怎么解决? - JKO
@JKO 你是在说哪个脚本?是 visreg 还是手动的那个? - COOLSerdash
visreg 脚本。 - JKO
@JKO 我也是这么想的。visreg 绘制的是偏残差而不是原始数据。你可以在这里找到有关它的教程和更多解释。如果你需要原始数据,我建议使用基本图形。 - COOLSerdash
1
一个更新 - 您可以简单地设置 partial = F,然后像在 R 基本绘图中一样添加点。 - JKO
显示剩余2条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接