random effects probit with missing data

  • Manny Gomez

    Manny Gomez - 2014-05-07

    I'm trying to estimate a random effects probit with missing response (y) using a latent Normal approach (Albert and Chib 1993).
    The model is below:


    for (i in 1:N){ # N - No. of individuals

      z[i] ~ dnorm(mu[i],1)T(low[y[i]+1],high[y[i]+1])
      mu[i] <- beta[1] + u[cluster[i]] +
        beta[2]*treat[i] + beta[3]*x3[i] + beta[4]*x4[i]        

    low[1]<--20 ; low[2]<-0 ; high[1]<-0 ;high[2]<-20

    Level-2 random effects

    for (j in 1:J){ # J - No. of studies
    u[j] ~ dnorm(zero, tau.u)

    Prior on level-2 variance

    tau.u ~ gamma(0.001, 0.001)
    sigma2.u <-inverse(tau.u)

    Prior on betas

    beta[1:4] ~ dmnorm(b0[1:4] , B0[1:4,1:4])

    When I have fully-observed outcome, the model runs fine, but with missing outcome, it fails to compile (even when I provide starting values for the missing 'y'):

    Compilation error on line 4.
    Unable to resolve node y[102]
    This may be due to an undefined ancestor node or a directed cycle in the graph

    Any help would be much appreciated.

  • Martyn Plummer

    Martyn Plummer - 2014-05-08

    You have mixed censoring and truncation. What you really want is this:

    y[i] ~ dinterval(z[i], 0)
    z[i] ~ dmnorm(mu[i], 1)

    But in any case you don't need to set up the Albert and Chib latent variable model yourself. Just load the glm module, define your model as a binary glm:

    y[i] ~ dbern(p[i])
    probit(p[i]) <- mu[i]

    and the Albert-Chib sampler will do block updating of beta[] and u[] together.

  • Manny Gomez

    Manny Gomez - 2014-05-08

    Thanks Martyn.
    I think I need the 1st approach because I would like to calculate the residual of the latent variable:
    Can I calculate this with the glm model you suggest?

    Anyway, I was not able to run your proposed latent variable model. See error below:

    << Error in node y2[1]
    Observed node inconsistent with unobserved parents at initialization.
    Try setting appropriate initial values. >>

    I tried different initial values for tau and beta, but did not work.
    Any ideas?

    • Krzysztof Sakrejda

      JAGS is very nicely telling you that the initial values for y2[1] (or its parents) are not consistent with whatever data is given as 'z'. Please check them. If you are still getting this error and you are sure the values are right, post the relevant values (z, the value of the node with the error, the values of all parents for the node with the error, and the model).

      I'm almost certain the different values you tried did not constitute a set that would produce a non-zero posterior probabilty for your data. Very (very) rarely there is an actual problem with JAGS... but for all the many times I've had to fix this problem in my models, it was never a bug.

      Perry deValpine did have a package out a little while ago which can be used to visualize BUGS models and it might help you out if you're not able to find the issue.

      • Manny Gomez

        Manny Gomez - 2014-05-08

        Thanks Krzysztof.

        'z' is a latent variable, and by definition not observed in the data. Only 'y' is observed
        if z[i]>0 then y=1, and vice-versa.

        Martyn suggested the code above, but perhaps something is missing. Imagine the simplest linear predictor with only the intercept and random effect. If I set initial values such as z>0 then when y[i]=0 it will give me that error message. When initial values so that z<0, then when y[i]=1 there will also be an inconsistency.

        Does that make sense?

        • Krzysztof Sakrejda

          "If I set initial values such as z>0 then when y[i]=0 it will give me that error message. "

          Yes, you must set initial values for all z[i]'s such that they are consistent with the observed y's. Settting them to something like
          z <- y - 0.5 in R should make it go (even if it's a bad starting point).

          • Manny Gomez

            Manny Gomez - 2014-05-09

            Thanks so much for your patience and help, it worked.
            I should have been smarter than that.


            • Krzysztof Sakrejda

              You're welcome! Glad it worked out.


Log in to post a comment.