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.