-1

I am attempting to plot fitted negative binomial results from the count portion of a hurdle regression model. Data (reproducable subset):

        Age     gender familysupport    bullying     Suicide. SuicideBinary NegBinSuicide
1   -0.771845      0       at risk    -0.34840000        1             1             1
2    0.228155      0       at risk     0.05160000        0             0            NA
3    0.228155      0     resilient     0.45160000        1             1             1
4    4.228155      0     resilient     0.25160000        0             0            NA
5   -2.771845      1     resilient    -0.44840000        0             0            NA
6   -2.771845      0       at risk    -0.64840000        0             0            NA
7   -0.771845      0     resilient    -0.04840000        0             0            NA
8   -0.771845      0     resilient    -0.14840000        0             0            NA
9   -0.771845      1       at risk    -0.64840000        0             0            NA
10   0.228155      0       at risk     0.05160000        0             0            NA
11   0.228155      0       at risk    -0.24840000        0             0            NA
12  -2.771845      0       at risk     0.15160000        0             0            NA
13  -0.771845      0     resilient    -0.14840000        0             0            NA
14  -1.771845      0       at risk    -0.44840000        0             0            NA
15   4.228155      0       at risk    -0.24840000        2             1             2
16   0.228155      0     resilient     0.05160000        1             1             1
17  -2.771845      0     resilient     0.05160000        0             0            NA
18   4.228155      0       at risk    -0.44840000        0             0            NA
19  -2.771845      1       at risk     0.25160000        0             0            NA
20  -1.771845      1       at risk    -0.54840000        0             0            NA
21  -0.771845      0     resilient    -0.14840000        0             0            NA
22  -0.771845      0       at risk    -0.34840000        0             0            NA
23  -2.771845      0     resilient    -0.14840000        0             0            NA
24   0.228155      1     resilient    -0.44840000        0             0            NA
25   0.228155      0       at risk    -0.14840000        2             1             2
26  -0.771845      0       at risk     1.95160000        0             0            NA
27  -2.771845      1       at risk    -0.44840000        1             1             1
28  -1.771845      0       at risk    -0.04840000        0             0            NA
29   2.228155      0     resilient    -0.44840000        0             0            NA
30  -0.771845      0       at risk    -0.34840000        0             0            NA
31   4.228155      0       at risk     0.15160000        4             1             4
32  -0.771845      1     resilient     0.15160000        1             1             1
33   0.228155      0     resilient     0.45160000        0             0            NA
34  -0.771845      0       at risk     0.15160000        0             0            NA
35  -0.771845      0       at risk    -0.04840000        0             0            NA
36   4.228155      0     resilient    -0.54840000        0             0            NA
37   0.228155      0     resilient     0.05160000        0             0            NA
38   1.228155      0       at risk    -0.34840000        1             1             1
39   2.228155      0       at risk     0.25160000        0             0            NA
40  -2.771845      0       at risk    -0.34840000        1             1             1
41   0.228155      0       at risk     1.75160000        2             1             2
42   4.228155      0       at risk     0.65160000        0             0            NA
43   0.228155      0     resilient     0.25160000       NA            NA            NA
44  -1.771845      0     resilient    -0.24840000        0             0            NA
45  -2.771845      0       at risk    -0.04840000        3             1             3
46  -0.771845      0     resilient     0.25160000        0             0            NA
47   3.228155      0     resilient     0.45160000        0             0            NA
48  -0.771845      0     resilient     0.85160000        0             0            NA
49  -2.771845      1       at risk     0.25160000        2             1             2
50  -0.771845      0       at risk     0.15160000        0             0            NA
51   1.228155      0     resilient    -0.44840000        0             0            NA
52   0.228155      0       at risk    -0.34840000        0             0            NA
53  -2.771845      0       at risk    -0.64840000       NA            NA            NA
54  -1.771845      0       at risk    -0.14840000        6             1             6
55   1.228155      0       at risk    -0.64840000        3             1             3
56   0.228155      0     resilient    -0.64840000        0             0            NA
57   2.228155      0     resilient    -0.64840000        0             0            NA
58   1.228155      0     resilient     1.05160000        0             0            NA
59   0.228155      0       at risk     0.25160000        0             0            NA
60  -0.771845      0       at risk     0.15160000        0             0            NA
61  -1.771845      0     resilient    -0.64840000        0             0            NA
62   1.228155      0       at risk    -0.44840000        0             0            NA
63   1.228155      0       at risk    -0.64840000        0             0            NA
64  -0.771845      0     resilient    -0.04840000        0             0            NA
65  -2.771845      0       at risk    -0.64840000        0             0            NA
66   1.228155      0       at risk     0.15160000        2             1             2
67   2.228155     NA     resilient    -0.64840000        0             0            NA
68   0.228155      0       at risk    -0.04840000       NA            NA            NA
69  -0.771845      0       at risk     0.05160000        0             0            NA
70  -0.771845      0       at risk    -0.64840000        0             0            NA
71  -0.771845      0     resilient    -0.64840000        0             0            NA
72  -0.771845      1       at risk     2.05160000       50             1            50
73   0.228155      0     resilient    -0.44840000        0             0            NA
74   1.228155      0     resilient     2.95160000        0             0            NA
75   0.228155      0     resilient     1.25160000        3             1             3
76   1.228155      0       at risk     0.45160000        2             1             2
77   0.228155      0     resilient             NA        0             0            NA
78   2.228155      0       at risk    -0.04840000        0             0            NA
79   2.228155      0       at risk    -0.64840000        2             1             2
80   0.228155      1     resilient     0.35160000        0             0            NA
81  -0.771845      0     resilient     0.25160000        2             1             2
82  -1.771845      1     resilient    -0.44840000        0             0            NA
83   0.228155      0       at risk    -0.64840000        0             0            NA
84   2.228155      0     resilient     0.01826667        0             0            NA
85   4.228155      0     resilient    -0.14840000        0             0            NA
86  -2.771845      0       at risk     0.25160000        0             0            NA
87   0.228155      0       at risk    -0.42617778        1             1             1
88   1.228155      0     resilient    -0.64840000        0             0            NA
89   0.228155      0     resilient    -0.04840000        0             0            NA
90   0.228155      0     resilient     0.15160000        0             0            NA
91   0.228155      0       at risk    -0.64840000        1             1             1
92   0.228155      0       at risk    -0.64840000        0             0            NA
93   4.228155      0     resilient    -0.34840000        0             0            NA
94   4.228155      0     resilient    -0.54840000        0             0            NA
95   1.228155      0     resilient     0.75160000        0             0            NA
96   3.228155      0       at risk    -0.24840000        0             0            NA
97  -2.771845      0     resilient    -0.64840000        0             0            NA
98  -1.771845      0     resilient     0.25160000        0             0            NA
99  -2.771845      0     resilient     0.35160000        2             1             2
100 -1.771845      0       at risk    -0.64840000        0             0            NA

However, when I use the following code:

library(VGAM)
##model
mod <- vglm(NegBinSuicide ~ Age + gender + bullying*familysupport, family=posnegbinomial())

library(visreg)
##plot with INTERACTION TERM
visreg(mod, "bullying", by="familysupport", xlab = "bullying", ylab = "Count model (number of suicide attempts)")

I get the following error message:

Error: $ operator not defined for this S4 class

I'm not sure what it means. Can anyone provide insight into how I might remedy this?

Ultimately, I'm trying to plot hurdle regression output for both components, as my interaction term was significant in each. (perhaps not given the sample data above)

library(pscl)
##Model
FullModel <- hurdle(Suicide. ~ Age + gender + bullying*familysupport | Age + gender + bullying*familysupport, dist = "negbin", link = "logit")

I want to create separate plots for each hurdle component. I've been able to do so for the logit portion (estimated separately) in visreg using data fitted in glm in MASS (which produced logistic results identical to the logit portion from the hurdle model), but using glm.nb for the NB portion produced different estimates. Hence my decision to switch to vglm from VGAM--estimates were same as hurdle, but plotting error emerged.

Any insight into how to work around the error, or how to plot these data otherwise would be greatly appreciated.

Nick
  • 11
  • 2
  • 1
    Please provide a [minimal reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) and indicate the non-base packages used. – digEmAll Aug 25 '14 at 14:52
  • You can find possible solution to this answer here: http://stackoverflow.com/questions/10961842/how-to-define-the-subset-operators-for-a-s4-class. You might also try searching SO for related questions by just pasting in your error message to the search bar. – Tommie C. Aug 25 '14 at 15:00
  • If you're going to muck with S4 classes, it really will be worth your time to read thru the R-intro and other manuals provided at CRAN – Carl Witthoft Aug 25 '14 at 15:34

1 Answers1

1

I suspect that the visreg package doesn't support VGAM::vglm fits. The help page for ?visreg says:

fit: The fitted model object you wish to visualize. Any object with 'predict' and 'model.frame' methods are supported, including lm, glm, gam, rlm, coxph, and many more.

Knowing that the pscl package has well-supported negative binomial hurdle models, I tried this:

library(pscl)
## using the built-in "bioChemists" data set
hh <- hurdle(art~fem+mar,dist="negbin",data=bioChemists)
library(visreg)
visreg(hh)

which seems to work fine and might work for your case.

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453