3
$\begingroup$

EDIT: Added an reproducible example

For one of my models, it seems the coefficients and the graphed out model do not agree.

I'm working with adverse effects data, in which intense reactions are rare and the number of zeros is huge - thus using ZIP model. This is a reproducible code for the problem:

LogFile <- file("ITreatLog.txt", open="wt")
sink(LogFile, append=FALSE, split=TRUE, type = "output")

#Loading libraries
if (!require("pscl")) {
  install.packages("pscl")
  library("pscl")
}
if (!require("lme4")) {
  install.packages("lme4")
  library("lme4")
}
if (!require("ggplot2")) {
  install.packages("ggplot2")
  library("ggplot2")
}
rm(list=ls())

zinb <- read.csv("Reproducible.csv", fileEncoding="UTF-8-BOM", sep = ";")
zinb <- within(zinb, {
  Group.1 <- factor(Group.1)
})

#Replace negative values (Missing) with sysmis
zinb[zinb==-1]<-NA
zinb[zinb==-2]<-NA

var <- "CA_T"
covar <- "DaySleep.1"

f <- formula("CA_T.1 ~ Group.1 + Age + DaySleep.1")
m1 <- zeroinfl(f, data = zinb)

#Code for printing out a graph of the predictions        
scatterwidth=0.3
GraphData <- expand.grid(seq(4,12,1), factor(0:1), 27)
xname = "Hours of Sleep"
colnames(GraphData) <- c(covar, "Group.1", "Age")                         
GraphData$phat <- predict(m1, GraphData)
GraphData <- subset(GraphData, subset=(phat<=100)) #Intensities over 100 make no logical sense
plotline = ggplot(GraphData, aes(group=Group.1, x = !!ensym(covar), y = phat, colour = factor(Group.1)))+
          geom_point() + geom_line() + scale_linetype_discrete() + labs(x=xname, y="Intensity", color="Group") +
          scale_color_hue(labels=c("Sham", "Active")) + 
          theme(legend.position = "bottom",legend.background=element_rect(fill="transparent"), legend.margin = margin(t=-0.2,r=0,b=0,l=0, unit="cm"))
plotline <- plotline + scale_x_continuous(breaks=seq(4,12,1))

ggsave("Reproducible.png",plot=plotline, width=14, height = 14, unit="cm")

sink()

And the model summary for m1:

Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 2.386609   0.282484   8.449  < 2e-16 ***
Group.11    0.210817   0.061288   3.440 0.000582 ***
Age         0.015874   0.004769   3.328 0.000873 ***
DaySleep.1  0.105312   0.029839   3.529 0.000417 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -4.74756    2.61227  -1.817   0.0692 .
Group.11     0.06700    0.56561   0.118   0.9057  
Age          0.05143    0.04667   1.102   0.2704  
DaySleep.1   0.42271    0.26037   1.624   0.1045 

My interpretation of this is that longer sleep (DaySleep.1) is associated with more intense adverse effect (the estimate beign positive). However, then I graph the model out, and output by above code is this:

enter image description here

I'm stumped. For other models, the graphs look like I'd expect them to from the coefficients. With this one, the coefficients suggest that longer sleep increases the adverse effects, while the graph seems to suggest the opposite association.

What am I doing wrong?

Finally, a copy of the dataset the code uses:

CA_T.1;DaySleep.1;Age;Group.1
0;8;21;0
0;8.33;24;0
30;9;22;1
0;7;27;1
10;8.5;21;0
20;6.5;35;0
35;8;21;1
30;6;25;0
0;9.5;23;1
20;8.5;23;1
30;6.5;23;0
50;7.86;24;1
0;7.5;45;0
-1;8;22;1
80;8;26;1
55;8.5;22;1
-1;8.5;22;0
70;8.3;34;0
70;8.33;29;0
0;4.3;28;1
30;8;22;1
0;7.75;25;1
0;7;27;1
0;8.66;22;0
40;9;24;0
40;5;24;1
-2;7;23;0
0;8;25;0
0;9;29;0
0;9;26;1
5;7.25;25;0
70;7.83;20;1
0;10;26;0
50;7.5;43;1
10;6.45;20;0
50;8.25;30;1
40;6;43;0
0;8.3;31;1
0;8.5;31;0
80;6;27;1
10;6;36;0
0;8;34;1
0;8.5;18;0
0;6;21;1
85;7;22;0
0;8;25;1
0;7.5;40;0
20;7.5;28;0
0;7.5;40;1
40;7;22;1
5;6;25;1
0;7;28;1
0;9.33;28;0
0;6.16;25;0
0;6.25;29;0
-1;6.87;22;0
30;7.16;22;1
-1;7.3;25;1
0;8;27;1
70;7;29;0
$\endgroup$
  • 1
    $\begingroup$ I updated the model code - but basically, it's an zero inflated poisson model, which constitutes of two parts - the poisson and the binomial. I hope my edits gave you the necessary information :) $\endgroup$ – Aaron K. Feb 19 at 21:57
  • $\begingroup$ I am not quite sure and I might be wrong (lack of reproducible example) but perhaps it has something to do with what you specified in the predict() function? Does predict(m1, GraphData, type = "prob") change that behaviour? Also have a look at vignette("countreg", package="pscl") page 23, 4th paragraph. $\endgroup$ – Stefan Feb 20 at 4:45
  • $\begingroup$ Also, I think you might get some more traction on this question if you can provide data or a reproducible example. $\endgroup$ – StatsStudent Feb 20 at 15:42
  • $\begingroup$ @Stefan Now there is an reproducible example. Adding type="prob" to the predict makes the ggplot to fail. Looking into it, without the "prob" it outputs list of the predicted means, with the "prob" it outputs a matrix of propabilities. GGPLOT or I have no idea what to do with that. Thanks for your help this far :) $\endgroup$ – Aaron K. Feb 20 at 21:46
0
$\begingroup$

I had more or less the same doubt regarding the coefficient signs in the "inflated" part using stata (the ZI coefficients in R)... The second part give the effect on the likelihood to get a 0.

In the following example (using Stata), the number of fish caught depends positively of the number of fishmen but the displayed coefficient is negative. The given interpretation is the following: "The inflate coefficient for persons suggests that for each unit increase in person the log odds of an inflated zero decrease by .564."

$\endgroup$
  • $\begingroup$ If you have overfit multi-collinearity can cause coefficients to be uninterpretible. $\endgroup$ – Michael R. Chernick Jun 28 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.