Plotting fitted lines for two different groups from a glm model with an interaction

Question

I have the following model:

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

With the following output:

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

I would like to plot the data, using R base plot, with separately colored dots for each group (treatment and control), and separately colored lines fit to each group. Something like the plot on page 20 here would be ideal, but I don't know how to translate what was done in that example (which uses a Poisson model) to what I have here (a quasi-Poisson model).

A reproducible example is included below.

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

Answer 1

I'm going to post a possible solution using base R graphics without much explaining. The core idea is to use predict to generate the predicted values on the response (original) scale and then plotting them. Everything else is basically cosmetics. Personally, I'd use the excellent visreg package that generates these kinds of graphics with ease. The corresponding code using visreg is posted at the bottom of this answer.

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")
)

Using visreg:

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