-bayesmh-, likelihood evaluator, and random effects

[Jeph Herrin]

I have been trying to replicate some analyses I did in another software package using -bayesmh- and my own likelihood evaluator. My first attempt went well; I wrote a likelihood evaluator for a two part Poisson (hurdle) model:

program myhurdle
    version 14.0
    args lnf xb xc
    tempname sig
    tempvar lnfj pi linpzero linpnozero munozero logpnozero
    qui gen double `linpzero'=`xb'
    qui gen double `pi'=1/(1+exp(-`linpzero'))
    qui gen double `linpnozero'=`xc'
    qui gen double `munozero'=exp(`linpnozero')
    qui gen double `logpnozero'=log(`pi')-log(1-exp(-`munozero'))-`munozero' ///
                         - lngamma($MH_y2+1)+$MH_y2*log(`munozero')
    qui gen double `lnfj'=log(1-`pi') if $MH_touse
    qui replace `lnfj' = `logpnozero' if $MH_y1>0 & $MH_touse
    summarize `lnfj', meanonly
    if r(N) < $MH_n {
        scalar `lnf'=.
        exit
    }
    scalar `lnf'=r(sum)
end

gen days1=days>0

bayesmh (days1 x1) (days x1), llevaluator(myhurdle) prior({days1:} {days:} , flat) saving(sim, replace) dots

This worked fine. However, I would also like to introduce random effects for each part of the model, so I modified the evaluator

program myhurdleRE
    version 14.0
    args lnf xb xc m1 m2
    tempname sig
    tempvar lnfj pi linpzero linpnozero munozero logpnozero
    qui gen double `linpzero'=`xb'+`m1'
    qui gen double `pi'=1/(1+exp(-`linpzero'))
    qui gen double `linpnozero'=`xc'+`m2'
    qui gen double `munozero'=exp(`linpnozero')
    qui gen double `logpnozero'=log(`pi')-log(1-exp(-`munozero'))-`munozero' ///
                         - lngamma($MH_y2+1)+$MH_y2*log(`munozero')
    qui gen double `lnfj'=log(1-`pi') if $MH_touse
    qui replace `lnfj' = `logpnozero' if $MH_y1>0 & $MH_touse
    summarize `lnfj', meanonly
    if r(N) < $MH_n {
        scalar `lnf'=.
        exit
    }
    scalar `lnf'=r(sum)
end

and tried various version of this:

bayesmh (days1 x1 i.id) (days x1 i.id), noconstant llevaluator(myhurdleRE, parameters({m1}, {m2})) ///
    prior({days1:i.id} , normal({m1},{S1})) ///
    prior({m1}, normal(`B0',10000)) prior({S1},igamma(0.001,0.001)) ///
    prior({days:i.id} , normal({m2},{S2})) ///
    prior({m2}, normal(`C0',10000)) prior({S2},igamma(0.001,0.001)) ///
    block({days1:i.id},split) block({days:i.id},split) block({S1},gibbs) block({S2},gibbs) ///
    saving(simRE, replace) dots mcmcsize(1000)

However, each time, after running silently for several hours, I get no output other than an error message:

option m2 for parameter m1 is not allowed
r(198);

So - my questions are:

1. How to avoid the error? I tried removing the comma in the parameters option, restating the priors, etc, but always get the error.
2. Almost as helpful - is there a way to ask -bayesmh- to somehow parse the options and return any errors right away? Waiting several hours for a syntax error is very inefficient.
3. Have I in fact set up the random effects correctly? I cannot find any information on combining REs with -llevaluator-, so I'm just doing what seems intuitive.

thanks,
Jeph

[Nikolay Balov (StataCorp)]

Jeph,
the right way to supply the parameters {m1} and {m1} to your evaluator is to specify parameters({m1} {m2}), without comma between them.

Your myhurdleRE evaluator however does not seem to be quite correct. I don't see why you need to add m1 and m2 to xb and cb. These are already included indirectly through the factor covariates {days:i.id} and {days1:i.id}. Then, if you won't be using {m1} and {m1} in myhurdleRE you may not need to supply them at all.

I also notice that in your second bayesmh call you need to include prior specifications for {days:x1} and {days1:x1}. You may also have some data specific problems which I cannot see only from your code.

In general, I suggest you first run bayesmh with the dryrun option to verify the model specification. Then you may want to make a short run, for example with burnin(0) mcmcs(1), to verify the correctness of your evaluator. You may also need to sub-sample your data if the latter is too big. Only after you are convinced that the model works, you may let it run at default settings.

[Jeph Herrin]

Thanks Nikolay,

The -dryrun- option is what I was looking for, didn't expect to find it under the Model subheading; I eventually took a tiny sample of the data, which at least let me debug things, but that was suboptimal. And yes, I figured out I need priors for {days1: x1} and {days:x1} in the process.

But I'm not so sure about the model specification anymore. I want the random effects of the IDs to be normal at the logit or loglinear level, and your questions made me realize that this is not what I'm doing. That is, if I wrote

linpzero = `xb'

would this translate to (this is pseudo code, not Stata code):

linpzero = coeff1*x1 + RE_i

where RE_i is a normally distributed random error across IDs?

thanks much,
Jeph

[Nikolay Balov (StataCorp)]

Jeph,

that is correct, the `xb' variable in the evaluator represents coeff1*x1 + RE_i in the first likelihood equation days1 and `xc' represents the corresponding linear predictor for the second likelihood equation for days. There is nothing special you have to do in the evaluator to implement random-effects since that is handled in the prior specifications of {days1:i.id} and {days:i.id}.

[Jeph Herrin]

Nilolay,

Thanks for clarifying, that's exactly what I was trying to sort out.

cheers,
Jeph