# How to plot estimate + raw data of a Bayesian zero inflated poisson?

GENERAL QUESTION:

How to back-transform estimates from a zero-inflated poisson to obtain the original scale in R? (I tried exponential like for poisson but the results are wrong)

DETAILED REPRODUCIBLE EXAMPLE:

To be able to modify the plot at will I want reproduce the plot obtained with: plot(marginal_effects(MODEL_NAME), points=T)

As an example I used the data from https://stats.idre.ucla.edu/r/dae/zip/ https://stats.idre.ucla.edu/stat/data/fish.csv

deleting a couple of extreme points for aesthetically plot purposes:

fish<-read.csv(“https://stats.idre.ucla.edu/stat/data/fish.csv”)
fish$$camper<-as.factor(fish$$camper)
fish<-subset(fish, count!=149 & count!=65)


I ran the following ZERO-INFLATED model

mod_pois_zero = brms::brm(count ~ child + camper +(1| persons),
control = list(adapt_delta = 0.999, max_treedepth = 12),
iter = 6000,
cores = 4,
data = fish,
family = “zero_inflated_poisson”,
file=“mod_pois”)


I wrote the following script to plot estimate and raw data as for marginal_effects:

#extract posterior samples
d<-posterior_samples(mod_pois)

#select the variable(s) you want to plot
d_sub<-d[,c("b_Intercept", "b_child")]
#transform it to a matrix
d_sub1<-as.matrix(d_sub)

# to plot my simulation I have to extract the effect and add them to my original dataset
newdat<-data.frame(b_child = seq(-0,3, by = 0.025))
#if you choose the same length of your database you can add it easily
#put the entire range of your data

#create the matrix you are going to fill
Xmat <-model.matrix(~ b_child, data=newdat)
fitmat<-matrix(ncol=nrow(d), nrow=nrow(newdat))
#Fill the matrix
for(i in 1:nrow(d)) fitmat[,i] <- exp(Xmat%*%d_sub1[i,]) #poisson model
#for(i in 1:nrow(d)) fitmat[,i] <- Xmat%*%d_sub1[i,]      #gaussian model

#add the credible interval and fit
newdat$$lower<- apply(fitmat,1,quantile, prob=0.025) newdat$$upper<-apply(fitmat, 1, quantile, prob=0.975)
newdat$fit<-apply(fitmat, 1, quantile, prob=0.5) #newdat$fit<-exp(Xmat%*% fixef(mod_pois)) #fixef does not work in brms

#now you can either merge newdat with your original database or plot them directly

plot(fish$$count ~fish$$child,
xlab =("b"), ylab = ("a"))
lines(newdat$$b_child, newdat$$fit, lwd=2)
lines(newdat$$b_child, newdat$$lower)
lines(newdat$$b_child, newdat$$upper)


the script works fine for norm and poisson models but the estimates with zero_inflated are different from the ones calculated from the automatic plotting of the package (see figures).

plot_brms

my_plot

How can I back-transform zero inflated models? clearly exp() is not enough. Where can I find this information for the rest of the link functions ( e.g. skew normal, etc.)?

On a different note, the default plot plot(marginal_effects(MODEL), points=T) always cut the extreme points, any idea why?

Also if someone has a better strategy to plot raw data and estimates from brms I would be glad to know. Online I could not find any example that I was able to reproduce.

Thanks and best regards

• Note that questions about coding are off-topic here. I'd recommend editing your question to focus on the statistical part. – mkt Aug 2 at 15:17
• @mkt thanks for your comment, I will do it. My short statistical question would be: what operation shall I use to back transform estimates from a zero-inflated poisson model? – have fun Aug 3 at 17:55
• You write: "the estimates with zero_inflated are wrong". Could you explain why you think they are wrong, where do you see that, and what you think they should be equal to? – mastropi Aug 4 at 12:58
• @mastropi you are right, I changed the text. Rather than wrong I meant "different from the ones calculated by the package for the plot". I actually do not know which ones are correct, I assumed mine are wrong. – have fun Aug 5 at 13:59