Predictions from zero-inflated or hurdle models

Hi,

6/ you are correct - my mistake.


5b/ I'm still not sure what you mean by the uncertainty, but ignoring
the random effects will result in under prediction. For a Poisson
log-normal the mean  is  exp(a_1+v_1/2)  rather than  exp(a_1) for the
standard Poisson.

Cheers,

Jarrod



Quoting Ruben Arslan <rubenarslan at gmail.com> on Tue, 17 Mar 2015
19:21:06 +0000:

Hi,

two more questions, sorry.

6.
could it be that where you wrote
normal.evd<-function(x, mu, v){
exp(-exp(x))*dnorm(x, mu, sqrt(v))
}

you're not actually passing the variance but the standard deviation? so

the

functions should be

normal.evd<-function(x, mu, sd){
exp(-exp(x))*dnorm(x, mu, sd)
}
normal.zt<-function(x, mu, sd){
exp(x)/(1-exp(-exp(x)))*dnorm(x, mu, sd)
}

Because at least in the example, you were taking the square root here
already: sd_1<-mean(sqrt(m1$VCV[,1]))
Maybe this was a copy-paste from the source where you do it differently?

Or

I misunderstood...

An addition to question 5:
5. When I use the correct inverse link function in my old code, but don't
do double integration, I do get a "reasonable" amount of uncertainty

around

the reference factor. The predicted means are too low, which I put down

to

this approach ignoring the VCV. I was just wondering how this could

happen

(actually seeing the "reasonable" amount of uncertainty was what made me
hopeful that my approach might not be entirely wrong).

Best,

Ruben


On Tue, Mar 17, 2015 at 7:53 PM Ruben Arslan <rubenarslan at gmail.com>

wrote:

Hi Jarrod,

thanks for the extensive reply! This helps a lot, though it sounds like

I

was hubristic to attempt this myself.
I tried using the approach you mapped out in the function gist
<https://gist.github.com/rubenarslan/aeacdd306b3d061819a6> I posted. I
simply put the pred function in a loop, so that I wouldn't make any
mistakes while vectorising and since I don't care about performance at

this

point.

Of course, I have some follow up questions though.. I'm sorry if I'm

being

a pain, I really appreciate the free advice, but understand of course if
other things take precedence.

1. You're not "down with" developing publicly are you? Because I sure
would like to test-drive the newdata prediction and simulation

functions..

2. Could you make sure that I got this right: "When predictions are to

be

taken after marginalising the random effects (including the `residual'
over-dispersion) it is not possible to obtain closed form expressions."
That is basically my scenario, right? In the example I included, I also
had a group-level random effect (family). Ore are you talking about the
"trait" as the random effect (as in your example) and my scenario is
different and I cannot apply the numerical double integration procedure

you

posted?
To be clear about my prediction goal without using language that I might
be using incorrectly: I want to show what the average effect in the
response unit, number of children, is in my population(s). I have data

on

whole populations and am using all of it (except individuals that don't
have completed fertility yet, because I have yet to find a way to model
both zero-inflation and right censoring).

3. "Numerical integration could be extended to double integration in
which case covariance between the Poisson part and the binary part could
be handled." That is what you posted an example of and it applies to my
scenario, because I specified a prior R=list(V=diag(2), nu=1.002, fix=2)
and rcov=~idh(trait):units, random=~idh(trait):idParents?
But this double integration approach is something you just wrote
off-the-cuff and I probably shouldn't use it in a publication? Or is

this

in the forthcoming MCMCglmm release and I might actually be able to

refer

to it once I get to submitting?

4. Could I change my model specification to forbid covariance between

the

two parts and not shoot myself in the foot? Would this allow for a more
valid/tested approach to prediction?

5. When I use your method on my real data, I get less variation around

the

prediction "for the reference level" than for all other factor levels.
My reference level actually has fewer cases than the others, so this
isn't "right" in a way.
Is this because I'm not doing newdata prediction? I get the "right"
looking uncertainty if I bootstrap newdata predictions in lme4,
Sorry if this is children's logic :-)
Here's an image of the prediction
<http://rpubs.com/rubenarslan/mcmcglmm_pred> and the raw data
<http://rpubs.com/rubenarslan/raw>.

Many thanks for any answers that you feel inclined to give.

Best wishes,

Ruben

On Tue, Mar 17, 2015 at 5:31 PM Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:

Hi,

Sorry - I should have replied to this post earlier. I've been working
on predict/simulate methods for all MCMcglmm models (including
zero-inflated/altered/hurdle/truncated models) and these are now
largely complete (together with newdata options).

When predictions are to be taken after marginalising the random
effects (including the `residual' over-dispersion) it is not possible
to obtain closed form expressions. The options that will be available
in MCMCglmm are:

1) algebraic approximation
2) numerical integration
3) simulation

1) and 2) are currently only accurate when the random/residual effect
structure implies no covariance between the Poisson part and the
binary part.

1) is reasonably accurate for zero-inflated distributions, but can be
pretty poor for the remainder because they all involve zero-truncated
Poisson log-normal distributions and my taylor approximation for the
mean is less than ideal (any suggestions would be helpful).

2) could be extended to double integration in which case covariance
between the Poisson part and the binary part could be handled.

In your code, part of the problem is that you have fitted a zapoisson,
but the prediction is based on a zipoisson (with complementary log-log
link rather than logt link).

In all zero-inflated/altered/hurdle/truncated models

E[y] = E[(1-prob)*meanc]

where prob is the probabilty of a zero in the binary part and meanc is
the mean of a Poisson distribution (zipoisson) or a zero-truncated
poisson (zapoisson and hupoisson). If we have eta_1 as the linear
predictor for the poisson part and eta_2 as the linear predictor for
the binary part:

In zipoisson: prob = plogis(eta_2)     and meanc = exp(eta_1)
In zapoisson: prob = exp(-exp(eta_2))  and meanc =
exp(eta_1)/(1-exp(-exp(eta_1)))
In hupoisson: prob = plogis(eta_2)     and meanc =
exp(eta_1)/(1-exp(-exp(eta_1)))
In ztpoisson: prob = 0                 and meanc =
exp(eta_1)/(1-exp(-exp(eta_1)))

In each case the linear predictor has a `fixed' part and a `random'
part which I'll denote as `a' and `u' respectively. Ideally we would
want

E[(1-prob)*meanc] taken over u_1 & u_2

if prob and meanc are independent this is a bit easier as

E[y] = E[1-prob]E[meanc]

and the two expectations ony need to taken with repsect to their
respective random effects. If we have sd_1 and sd_2 as the standard
deviations of the two sets of random effects then for the zapoisson:

   normal.evd<-function(x, mu, v){
      exp(-exp(x))*dnorm(x, mu, sqrt(v))
   }
   normal.zt<-function(x, mu, v){
     exp(x)/(1-exp(-exp(x)))*dnorm(x, mu, sqrt(v))
   }

   pred<-function(a_1, a_2, sd_1, sd_2){
     prob<-1-integrate(normal.evd, qnorm(0.0001, a_2,sd_2),
qnorm(0.9999, a_2,sd_2), a_2,sd_2)[[1]]
     meanc<-integrate(normal.zt, qnorm(0.0001, a_1,sd_1),
qnorm(0.9999, a_1,sd_1), a_1,sd_1)[[1]]
     prob*meanc
   }

#  gives the expected value with reasonable accuracy.  As an example:

   x<-rnorm(300)
   l1<-rnorm(300, 1/2+x, sqrt(1))
   l2<-rnorm(300, 1-x, sqrt(1))

   y<-rbinom(300, 1, 1-exp(-exp(l2)))
   y[which(y==1)]<-qpois(runif(sum(y==1), dpois(0,
exp(l1[which(y==1)])), 1), exp(l1[which(y==1)]))
   # cunning sampler from Peter Dalgaard (R-sig-mixed)

   data=data.frame(y=y, x=x)
   prior=list(R=list(V=diag(2), nu=0.002, fix=2))

   m1<-MCMCglmm(y~trait+trait:x-1, rcov=~idh(trait):units, data=data,
family="zapoisson", prior=prior)

   b_1<-colMeans(m1$Sol)[c(1,3)]
   b_2<-colMeans(m1$Sol)[c(2,4)]
   sd_1<-mean(sqrt(m1$VCV[,1]))
   sd_2<-mean(sqrt(m1$VCV[,2]))

   # note it is more accurate to take the posterior mean prediction
rather than the prediction from the posterior means as I've done here,
but for illustration:

   x.pred<-seq(-3,3,length=100)
   p<-1:100
   for(i in 1:100){
     p[i]<-pred(a_1 = b_1[1]+x.pred[i]*b_1[2], a_2 =
b_2[1]+x.pred[i]*b_2[2], sd_1=sd_1, sd_2=sd_2)
   }

   plot(y~x)
   lines(p~x.pred)

Cheers,

Jarrod




Quoting Ruben Arslan <rubenarslan at gmail.com> on Tue, 17 Mar 2015
13:33:25 +0000:

Dear list,

I've made a reproducible example of the zero-altered prediction,
in the hope that someone will have a look and reassure me that I'm

going

about this the right way.
I was a bit confused by the point about p_i and n_i being correlated

(they

are in my case), but I think this was a red herring for me
since I'm not deriving the variance analytically.
The script is here: https://gist.github.com/rubenarslan/

aeacdd306b3d061819a6

and if you don't want to run the simulation fit yourself, I've put an

RDS

file of the fit here:
https://dl.dropboxusercontent.com/u/1078620/m1.rds

Best regards,

Ruben Arslan

On Tue, Mar 10, 2015 at 1:22 PM Ruben Arslan <rubenarslan at gmail.com>

wrote:

Dear Dr Bolker,

I'd thought about something like this, one point of asking was to
see whether
a) it's implemented already, because I'll probably make dumb

mistakes

while

trying b) it's not implemented because it's a bad idea.
Your response and the MCMCglmm course notes make me hope that it's

c)

not

implemented because nobody did yet or d) it's so simple that

everybody

does

it on-the-fly.

So I tried my hand and would appreciate corrections. I am sure

there is

some screw-up or an inelegant approach in there.
I included code for dealing with mcmc.lists because that's what I

have

and

I'm not entirely sure how I deal with them is correct either.

I started with a zero-altered model, because those fit fastest and
according to the course notes have the least complex likelihood.
Because I know not what I do, I'm not dealing with my random

effects at

all.

I pasted a model summary below to show what I've applied the below
function to. The function gives the following quantiles when applied

to 19

chains of that model.
         5%         50%         95%
5.431684178 5.561211207 5.690655200
         5%         50%         95%
5.003974382 5.178192327 5.348246558

Warm regards,

Ruben Arslan

HPDpredict_za = function(object, predictor) {

if(class(object) != "MCMCglmm") {
if(length( object[[1]]$Residual$nrt )>1) {
object = lapply(object,FUN=function(x) { x$Residual$nrt<-2;x })
}
Sol = mcmc.list(lapply(object,FUN=function(x) { x$Sol}))
vars = colnames(Sol[[1]])
} else {
Sol = as.data.frame(object$Sol)
vars = names(Sol)
}
za_predictor = vars[ vars %ends_with% predictor & vars %begins_with%
"traitza_"]
za_intercept_name = vars[ ! vars %contains% ":" & vars %begins_with%
"traitza_"]
 intercept = Sol[,"(Intercept)"]
za_intercept = Sol[, za_intercept_name]
l1 = Sol[, predictor ]
l2 = Sol[, za_predictor ]
if(is.list(object)) {
intercept = unlist(intercept)
za_intercept = unlist(za_intercept)
l1 = unlist(l1)
l2 = unlist(l2)
}
 py_0 = dpois(0, exp(intercept + za_intercept))
y_ygt0 = exp(intercept)
at_intercept = (1-py_0) * y_ygt0

py_0 = dpois(0, exp(intercept + za_intercept + l2))
y_ygt0 = exp(intercept +  l1)
at_predictor_1 = (1-py_0) * y_ygt0
print(qplot(at_intercept))
print(qplot(at_predictor_1))
df = data.frame("intercept" = at_intercept)
df[, predictor] = at_predictor_1
print(qplot(x=variable,
y=value,data=suppressMessages(melt(df)),fill=variable,alpha=

I(0.40),

geom =

'violin'))
print(quantile(at_intercept, probs = c(0.05,0.5,0.95)))
print(quantile(at_predictor_1, probs = c(0.05,0.5,0.95)))
invisible(df)
}

summary(object[[1]])

 Iterations = 100001:299901
 Thinning interval  = 100
 Sample size  = 2000

 DIC: 349094

 G-structure:  ~idh(trait):idParents

                      post.mean l-95% CI u-95% CI eff.samp
children.idParents       0.0189   0.0164   0.0214     1729
za_children.idParents    0.2392   0.2171   0.2622     1647

 R-structure:  ~idh(trait):units

                  post.mean l-95% CI u-95% CI eff.samp
children.units        0.144    0.139    0.148     1715
za_children.units     1.000    1.000    1.000        0

 Location effects: children ~ trait * (maternalage.factor +

paternalloss +

maternalloss + center(nr.siblings) + birth.cohort + urban + male +
paternalage.mean + paternalage.diff)

                                           post.mean  l-95% CI

u-95%

CI

eff.samp   pMCMC
(Intercept)                                 2.088717  2.073009

2.103357

  2000 <0.0005 ***
traitza_children                           -1.933491 -1.981945

-1.887863

  2000 <0.0005 ***
maternalage.factor(14,20]                   0.007709 -0.014238

0.027883

  1500   0.460
maternalage.factor(35,50]                   0.006350 -0.009634

0.024107

  2000   0.462
paternallossTRUE                            0.000797 -0.022716

0.025015

  2000   0.925
maternallossTRUE                           -0.015542 -0.040240

0.009549

  2000   0.226
center(nr.siblings)                         0.005869  0.004302

0.007510

  2000 <0.0005 ***
birth.cohort(1703,1722]                    -0.045487 -0.062240

-0.028965

  2000 <0.0005 ***
birth.cohort(1722,1734]                    -0.055872 -0.072856

-0.036452

  2000 <0.0005 ***
birth.cohort(1734,1743]                    -0.039770 -0.056580

-0.020907

  2000 <0.0005 ***
birth.cohort(1743,1750]                    -0.030713 -0.048301

-0.012214

  2000   0.002 **
urban                                      -0.076748 -0.093240

-0.063002

  2567 <0.0005 ***
male                                        0.106074  0.095705

0.115742

  2000 <0.0005 ***
paternalage.mean                           -0.024119 -0.033133

-0.014444

  2000 <0.0005 ***
paternalage.diff                           -0.018367 -0.032083

-0.005721

  2000   0.007 **
traitza_children:maternalage.factor(14,20] -0.116510 -0.182432

-0.051978

  1876   0.001 ***
traitza_children:maternalage.factor(35,50] -0.045196 -0.094485

0.002640

  2000   0.075 .
traitza_children:paternallossTRUE          -0.171957 -0.238218

-0.104820

  2000 <0.0005 ***
traitza_children:maternallossTRUE          -0.499539 -0.566825

-0.430637

  2000 <0.0005 ***
traitza_children:center(nr.siblings)       -0.023723 -0.028676

-0.018746

  1848 <0.0005 ***
traitza_children:birth.cohort(1703,1722]   -0.026012 -0.074250

0.026024

  2000   0.319
traitza_children:birth.cohort(1722,1734]   -0.279418 -0.329462

-0.227187

  2000 <0.0005 ***
traitza_children:birth.cohort(1734,1743]   -0.260165 -0.312659

-0.204462

  2130 <0.0005 ***
traitza_children:birth.cohort(1743,1750]   -0.481457 -0.534568

-0.426648

  2000 <0.0005 ***
traitza_children:urban                     -0.604108 -0.645169

-0.562554

  1702 <0.0005 ***
traitza_children:male                      -0.414988 -0.444589

-0.387005

  2000 <0.0005 ***
traitza_children:paternalage.mean           0.006545 -0.018570

0.036227

  2000   0.651
traitza_children:paternalage.diff          -0.097982 -0.136302

-0.060677

  2000 <0.0005 ***


On Mon, Mar 9, 2015 at 10:12 PM Ben Bolker <bbolker at gmail.com>

wrote:

Ruben Arslan <rubenarslan at ...> writes:

Dear list,

I wanted to ask: Is there any (maybe just back of the envelope)

way

to

obtain a response prediction for zero-inflated or hurdle type

models?

I've fit such models in MCMCglmm, but I don't work in ecology

and my

previous experience with explaining such models to "my audience"

did not

bode well. When it comes to humans, the researchers I presented

to

are

not

used to offspring count being zero-inflated (or acquainted with

that

concept), but in my historical data with high infant mortality,

it

is

(in

modern data it's actually slightly underdispersed).

Currently I'm using lme4 and simply splitting my models into two

stages

(finding a mate and having offspring).
That's okay too, but in one population the effect of interest is

not

clearly visible in either stage, only when both are taken

together

(but

then the outcome is zero-inflated).
I expect to be given a hard time for this and hence thought I'd

use

a

binomial model with the outcome offspring>0 as my main model, but

that

turns out to be hard to explain too and doesn't
really do the data justice.

Basically I don't want to be forced to discuss my smallest

population

as a

non-replication of the effect because I was insufficiently able

to

explain

the statistics behind my reasoning that the effect shows.

  I think the back-of-the envelope answer would be that for a

two-stage

model with a prediction of p_i for the probability of having a

non-zero

response (or in the case of zero-inflated models, the probability

of

_not_ having a structural zero) and a prediction of n_i for the
conditional
part of the model, the mean predicted value is p_i*n_i and the
variance is _approximately_ (p_i*n_i)^2*(var(p_i)/p_i^2 +

var(n_i)/n_i^2)

(this is assuming
that you haven't built in any correlation between p_i and n_i,

which

would be hard in lme4 but _might_ be possible under certain

circumstances

via a multitype model in MCMCglmm).

  Does that help?

      [[alternative HTML version deleted]]

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.