Skip to content

Zero Inflated Models - Plotting Results

2 messages · lieslpe, Achim Zeileis

lieslpe
# 
Dear SamiC,
I am also attempting to plot my zero inflated model on my data. Did you find
a solution? Does anyone else on this list have a solution?

Thanks,
Liesl 

Message from SamiC Jun 30, 2011:
I am fitting a zero inflated negative binomial model to my data.  I have
pretty much got my selected model and now i am wanting to plot the model to
the origional data.

Example.

f8<-formula(Birds~Tide+Fish+Zooplankton+Chlorophylla|StratificationIndex)
Nb8<-zeroinfl(f8, dist="negbin", link= "logit", data=ocean)

Tide is a factor, so i assume i have to plot a different graph for each
level of the factor.

what i essentially want to do is plot a graph for each variables against
birds with the fitted line of the model.

I have been using the predict function, but i get the same trend for every
graph and variable.  I was reading that the predict function gives a mean
across all the values (or something to this effect).

Does anyone know how to code for this in R.  from the above model i want to
plot birds~fish (the original data) and then fit the line of the
relationship from the model.




--
View this message in context: http://r.789695.n4.nabble.com/Zero-Inflated-Models-Plotting-Results-tp3636373p4644800.html
Sent from the R help mailing list archive at Nabble.com.
Achim Zeileis
#
On Tue, 2 Oct 2012, lieslpe wrote:

            

      
      
      
    

        
If you want to compare observed and fitted frequencies for the counts 0, 
1, 2, ..., then a common approach is use a rootogram or some variant of 
it. Below is a code snippet that works for "zeroinfl" and "hurdle" objects 
from pscl. But it's fairly straightforward to extend this yourself.

If you want to visualize dependence on covariates, you can either use 
scatterplots for continuity-corrected logs (as the countreg vignette of 
pscl does) or separately visualize factor(y > 0) ~ x and log(y) ~ x, 
subset = y > 0 or so.

hth,
Z

rootogram <- function(obj, max = NULL, ...) {
   y <- model.response(model.frame(obj))
   tab <- table(factor(y, levels = 0:max(y)))
   tab2 <- colSums(predprob(obj))
   if(is.null(max)) {
     max <- if(all(tab2 >= 1)) max(y) else max(ceiling(mean(y)), min(which(tab2 < 1)) - 1)
   }
   max <- min(max, length(tab) - 1) + 1
   obsrvd <- sqrt(tab[1:max])
   expctd <- sqrt(tab2[1:max])
   res <- obsrvd - expctd
   x <- barplot(obsrvd, offset = -res, xlab = "Count", ylab = "sqrt(Frequency)")
   lines(x, expctd, col = 2, type = ifelse(max > 25, "l", "b"), pch = 19)
   abline(h = 0)
   invisible(cbind(observed = tab, expected = tab2))
}