[R-sig-ME] Predictive growth curve for a ZIP model in MCMCglmm

[Christopher David Desjardins]

Hi,
I have the following model:

q2.schint<- MCMCglmm(sus~trait-1 + at.level(trait,1):grade + 
at.level(trait,1):I(grade^2) + at.level(trait,1):male.f + 
at.level(trait,1):ethnic.f + at.level(trait,1):sped.f 
+at.level(trait,1):risk.f + at.level(trait,1):male.f*grade + 
at.level(trait,1):ethnic.f*grade + at.level(trait,1):sped.f*grade 
+at.level(trait,1):risk.f*grade + at.level(trait,1):male.f*I(grade^2) + 
at.level(trait,1):ethnic.f*I(grade^2) + 
at.level(trait,1):sped.f*I(grade^2) 
+at.level(trait,1):risk.f*I(grade^2), random=~us(at.level(trait,1)):id.f 
+ us(at.level(trait,1)):schn.f + us(at.level(trait,1)):schn.f:id.f, 
data=suslm, rcov=~idh(trait):units, family="zipoisson", 
prior=prior,nitt=30000, thin=50, burnin=10000)

I would like to generate group level growth curves by "risk.f". risk.f 
is a categorical variable with 4 levels. I curious how I might do this? 
Additionally, how can I find out the proportion of zeros predicted by 
the Poisson component and the proportion of zeros predicted by the 
Poisson+ZIP component?

Jarrod Hadfield has provided the following code in the past but I am not 
sure how to generalize to my example.

Thanks!
Chris


beta<-colMeans(m3$Sol)
c2<-16*sqrt(3)/(15*pi)             # see Diggle 2004
zi<-inv.logit(beta[2]/sqrt(1+c2))  # scale zero-inflation as if residual 
variance was set to zero
pred<-exp(beta[1]+beta[3]*I(1:12)+beta[4]*I((1:12)^2)+0.5*mean(rowSums(m3$VCV[,1:2]))) 
# poisson predictions
pz<-ppois(0, pred)  # proportion of zeros predicted from the poisson
tz<-zi+(1-zi)*pz  # total proportion of zeros predicted
Ezip0<-function(mu){mu/(1-exp(-mu))} # the expected value of a zero 
truncated Poisson with mu = mean prior to trunaction.

# plot for non-zero data points
plot(suslm$sus[which(suslm$sus>0)]~as.factor(suslm$grade)[which(suslm$sus>0)],xlab="Grade",ylab="Number 
of Suspensions")
points(tapply(suslm$sus[which(suslm$sus>0)], 
as.factor(suslm$grade)[which(suslm$sus>0)], mean)~I(1:12), pch=16)  # 
overlay observed means
lines(Ezip0(pred)~I(1:12)) # predictions goodish - see below