OpenBUGS error messages. Something went right.

Although my relationship with OpenBUGS is still in the early stages, one glaringly obvious and quite infuriating shortcoming has been the lack of detail in its error messages. The initiated amongst you are no doubt familiar with the phrase “Sorry, something went wrong in procedure [blah]”. These error messages had me chasing my tail for hours in vain attempts to find offending segments of code. Surely, I thought, there must be a way to display more detailed information about the cause of an error.

Default OpenBUGS error messages

A default OpenBUGS error message

As it happens, there is.

The solution — enabling “Oberonmicrosystems’ extremely powerful trap handler” — is in fact mentioned briefly in the OpenBUGS Developer Manual. It’s straightforward enough… simply locate the file Bugs/Code/Traphandler.ocf, and rename it (e.g. to Traphandler_simple.ocf). Once the file is renamed, quit and restart OpenBUGS, which will now use the Oberonmicrosystems trap handler instead of the frustrating default handler. The result: OpenBUGS now spews forth copious amounts of information relating to an error, including the procedure that ran into the error (and its source code), and the values of all local variables in that procedure. If at any time you feel the need to revert to the original trap handler, just restore the original file name.

The Oberonmicrosystems trap handler

The Oberonmicrosystems trap handler

In the above screenshot, the trap handler indicates that OpenBUGS tripped up when it attempted to sample F[26, 5]. In this case it was because the precision parameter for the normal distribution was set to infinity (clearly the standard deviation of the distribution can’t be calculated if 1/σ2 is equal to infinity). For more elusive errors, the source code (below), which has the problem line highlighted, can help to reveal the problematic procedures, parameters or data.

Oberonmicrosystems traphandler - procedure source code

Procedure source code

More information about OpenBUGS trap windows can be found in the OpenBUGS manual.

Advertisements

30 thoughts on “OpenBUGS error messages. Something went right.

  1. Hi, This is great news. OBUGS sometimes gets on my nerves. The last simulation I ran it just said “Sorry, something went wrong.” Very frustrating, after waiting for several hours for the simulation to finish.

    Do you know how can I do this in Ubuntu? I search for the OpenBUGS files or folder, but I could only find one library (ELF-file) and the compiled C program in /usr/loca/bin/. There seems to be nothing else.

    Thanks again!

    • Luke,

      Sorry for the slow reply. When using Linux, I run the Windows release of OpenBUGS (through Wine) rather that the Linux release. I do this because it provides the GUI, which I find useful for initial model development/monitoring (live trace plots, etc.). (Typically I then move to running models directly from R, where more often than not I use JAGS and rjags.) If you run through Wine, then the file I refer to above can be found at, e.g. /home/john/.wine/drive_c/Program Files/OpenBUGS/OpenBUGS322/Bugs/Code/Traphandler.ocf.

      OpenBUGS for Linux provides only the CLI, in which case I suspect this post does not apply. Apologies for the confusion!

  2. thank you very much for this post actually, i have same problem and the problem line highlighted in my code is “i < LEN(m.parents)" what is that means?

    • Hi thanoon. It’s been a long time since I’ve used OpenBUGS (I use JAGS from R these days), so I can’t tell what’s wrong with your model. You could try to run it in JAGS to see if you get a more informative (or at least different) error message. The code is usually equivalent. Sorry I can’t be of more help.

  3. i am using also R2JAGS and also have same problem “Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
    Error in node y1[1,1]
    Unobserved node inconsistent with unobserved parents at initialization.
    Try setting appropriate initial values”
    i used inits=null to generate the initial values so where is the problem??

    • From the JAGS manual, if you use inits=NULL, JAGS will use a

      ‘typical value’ from the prior distribution. The exact meaning of ‘typical value’ depends on the distribution of the stochastic node, but is usually the mean, median, or mode.

      The error suggests that one of your priors might be inappropriate, allowing values that are incompatible with the prior of a downstream node. Consider this crude example:

      M <- function() {
        for (i in 1:1000) {
          y[i] ~ dbin(p[i], 100)
          p[i] ~ dnorm(0, 0.0001)
        }
      }
      
      library(R2jags)
      jags(list(y=rbinom(1000, 100, 0.8)), NULL, 'p', M)
      

      By assigning a normal prior to p, it allows p to take values outside of the [0, 1] interval. This is a problem, since p is the probability parameter for the binomial distribution of Y. The above model will throw the error Error in node y[1] ... Observed node inconsistent with unobserved parents at initialization.

      So make sure your priors are sensible, and specify initial values that are consistent with the model.

  4. Thank you very much for your explanation, this is the priors and the precision, I think the problem in the priors but how can i discover this problem and where is it please?

    #priors on loadings and coefficients
    for(i in 1:18) {
        mu.y1[i]~dnorm(0.0,4.0)
    }
    for(i in 1:15) {
         lam1[i]~dnorm(0.8,4.0)
    } 
    for(i in 1:5) {
         gam1[i]~dnorm(0.6,4.0)
    }
    for(i in 1:18) {
         mu.y2[i]~dnorm(0.0,4.0)
    } 
    for(i in 1:15) {
         lam2[i]~dnorm(0.8,4.0)
    } 
    for(i in 1:5) {
         gam2[i]~dnorm(0.6,4.0)
    }
    
    #priors on precisions
    for(j in 1:P) {
         psi1[j]~dgamma(10,8)
         sgm1[j]<-1/psi1[j]
    }
    
    psd1~dgamma(10,8)
    sgd1<-1/psd1
    phi1[1:2,1:2]~dwish(R[1:2,1:2], 30)
    phx1[1:2,1:2]<-inverse(phi1[1:2,1:2])
    
    for(j in 1:P) {
        psi2[j]~dgamma(10,8)
        sgm2[j]<-1/psi2[j]
    }
    
    psd2 ~ dgamma(10,8)
    sgd2 <- 1/psd2
    phi2[1:2,1:2] ~ dwish(R[1:2,1:2], 30)
    phx2[1:2,1:2] <- inverse(phi2[1:2,1:2])
    
    }', #end of model
    

    Many thanks in advance

  5. thank you very much for your response
    this is the dist. of y1 and y2:

    for (j in 1:P) {
        y1[i,j]~dnorm(mu1[i,j],psi1[j]) T(thd[j,z1[i,j]], thd[j,z1[i,j]+1])
    }
    for (j in 1:P) {
        y2[i,j]~dnorm(mu2[i,j],psi2[j]) T(thd[j,z2[i,j]], thd[j,z2[i,j]+1])
    }
    
    • When JAGS selects initial values for unobserved nodes (when you set inits=NULL), it selects “typical values” from the distribution, but doesn’t seem to take into account any truncation. If JAGS is using mu1[i, j] as the prior for y1, there’s a chance that this is outside the limits imposed by the truncation.

      The solution is to specify initial values for any unobserved nodes of y1 and y2.

  6. thank you for your help
    Actually, this is my problem now i have the following initial values as following”

    Init 1

    testJAGSinits1 = list(
    mu.y1=rep(0.2, 18), lam1=rep(0.4, 15), psi1=rep(1.0, 18),
    psd1=1.0, gam1=rep(1.0,5),
    phi1=structure(
    .Data=c(1.0, 0.0,0.0, 1.0),
    .Dim=c(2,2)),

                  mu.y2=rep(0.6, 18), lam2=rep(0.5, 15), psi2=rep(1.0, 18), psd2=1.0,
                  gam2=rep(1.0,5), 
                  phi2=structure(
                       .Data=c(1.0, 0.0,0.0, 1.0),
                       .Dim=c(2,2)),
    
                  xi1=structure( 
                       .Data=rep(0, 400),
                       .Dim=c(200,2)),
    
                  xi2=structure( 
                       .Data=rep(0, 400), 
                       .Dim=c(200,2)))
    

    Init 2

    testJAGSinits2 = list(
    mu.y1=rep(0.5,18), lam1=rep(0.5, 15), psi1=rep(0.5, 18), psd1=0.6,
    gam1=c(0.0, 0.0, 0.0,0.0,0.0),

                  phi1=structure( .Data=c(0.5, 0.0, 0.0, 0.5), 
                                  .Dim=c(2,2)), 
    
                  mu.y2=rep(0.5, 18), lam2=rep(0.5, 15), psi2=rep(0.5, 18),  psd2=0.6, 
                  gam2=c(0.0, 0.0, 0.0,0.0,0.0),
    
                  phi2=structure( .Data=c(0.5, 0.0, 0.0, 0.5), 
                                  .Dim=c(2,2)), 
                  xi1=structure( 
                       .Data=rep(0, 400),
                       .Dim=c(200,2)),
    
                  xi2=structure( 
                       .Data=rep(0, 400),
                       .Dim=c(200,2)))
    

    How can i set initial values for any unobserved nodes of y1 and y2?

    many thanks in advance

    • I assume y1 is a matrix of data with missing values… If this is the case, for inits you should provide a matrix with the same dimensions as y1, containing NA values where y1 has data, and valid initial values where y1 is NA. Same concept for y2.

  7. In my own case y1 defined with another observed data matrix z1 and i am using mean of the matrix y as a dependent var. to write structural equations. I tried to set y1 and y2 in the initial values but the problem still same as a previous. I am sure the problem in the priors dist. of y1 and y2 respectively.
    many thanks

  8. I am running a survival model below with time for the dead following a weibull distribution. On udating, I get this error, sorry something went wrong in procedure updater: Mode in module updater rejection. How can I solve the error. Thx Betty

    model
    {
    for(i in 1 : 7878) {
    #likelihood
    t[i] ~ dweib(r, mu[i])I(tcen[i],)
    log(mu[i]) <- alpha + b.visit_num * visit_num[i] +b.reggregg[i]+b.visit1 *visit1[i]
    + b.place
    place[i] + b.ass_deliveryass_delivery[i]+
    b.respond_slept *respond_slept[i]+
    b.respond_type
    respond_type[i]+
    b.wealth_indwealth_ind[i]+
    b.age_group
    age_group[i]+
    b.age_firstage_first[i]+
    b.educ_level
    educ_level[i]+
    b.readread[i]+
    b.live_births
    live_births[i]+
    b.birth_ordbirth_ord[i]+
    b.smoke
    smoke[i]+
    b.occupat*occupat[i]

    }

    Priors

    r~ dexp(0.001)
    alpha ~ dnorm(0.0, 0.0001)
    b.regg~ dnorm(0.0, 0.0001)
    b.visit_num ~ dnorm(0.0, 0.0001)
    b.visit1 ~ dnorm(0.0, 0.0001)
    b.place~ dnorm(0.0, 0.0001)
    b.ass_delivery~ dnorm(0.0, 0.0001)
    b.respond_slept ~ dnorm(0.0, 0.0001)
    b.respond_type ~ dnorm(0.0, 0.0001)
    b.wealth_ind ~ dnorm(0.0, 0.0001)
    b.age_group ~ dnorm(0.0, 0.0001)
    b.age_first ~ dnorm(0.0, 0.0001)
    b.educ_level ~ dnorm(0.0, 0.0001)
    b.read ~ dnorm(0.0, 0.0001)
    b.live_births ~ dnorm(0.0, 0.0001)
    b.birth_ord ~ dnorm(0.0, 0.0001)
    b.smoke ~ dnorm(0.0, 0.0001)
    b.occupat~ dnorm(0.0, 0.0001)
    }

    • Sorry Betty, but I don’t know what might cause that error. I suggest simplifying your model (remove most of the predictors) and then, if the reduced model runs fine, iteratively add predictors back in and see if there’s a particular predictor that breaks it.

      • Thanks John let me try that and see. Will let you know the progress.

        Thank you Betty

  9. Dear John, I face the following error on loading data, array index is
    greater than array upper bound for w. I am assessing random effects
    due to clusters, w[dhsclust[i]] , where w is parameter in the model
    and dhsclust is the cluster variable. I have N=404 clusters in my data
    and M=7878 observations. I have defined M and N in the data list as
    N=404 and M=7878. Also in the initials I have initialized w as a
    vector of 404 zeros i.e. w=c(0,0,0,0,0,……………….) but in the
    data dhsclust vector has 7878 observations with cluster numbers
    repeating them selves. Where could I have gone wrong. Please advise.
    Below is my model, data and initials.

    Thank you
    Betty

    model {

    for (i in 1:M)
    {

    failure[i] ~ dbern(p[i])
    logit(p[i]) <- b0 + part1[i] + part2[i] +part3[i] +part4[i] +part5[i]
    + part6[i] +part7[i]
    +w[dhsclust[i]]

    part1[i]<-b[1]*regg[i]

    part2[i]<-b[2]*hhnet_own[i]+b[3] *respond_slept[i]

    part3[i]<-b[4]*ass_delivery[i] + b[5] *equals(visit_num[i],2)+ b[6] *
    equals(visit_num[i],3)+b[7] *equals(visit_num[i],4)

    part4[i]<-b[8]equals(wealth_ind[i],2)+b[9]equals(wealth_ind[i],3) +
    b[10]equals(wealth_ind[i],4)+ b[11]equals(wealth_ind[i],5)

    part5[i]<-b[12]equals(age_group[i],2) + b[13]equals(age_group[i],3)
    + b[14]equals(age_group[i],4) + b[15]equals(age_group[i],5) +
    b[16]equals(age_group[i],6) + b[17]equals(age_group[i],7) +
    b[18]equals(age_first[i],2) +
    b[19]
    equals(age_first[i],3)+ b[20]equals(age_first[i],4) +
    b[21]
    equals(educ_level[i],2)+b[22]equals(educ_level[i],3) +
    b[23]
    equals(educ_level[i],4)+b[24]read[i]+b[25]preg_term[i]

    part6[i]<-b[26]mode_delivery[i] + b[27]equals(birth_ord[i],2) +
    b[28]equals(birth_ord[i],3) + b[29]child_age[i] + b[30]*sex_child[i]

    part7[i]<- b[31]smoke[i]+ b[32]equals(occupat[i],2) +
    b[33]equals(occupat[i],3) + b[34]equals(occupat[i],4)

    }

    b0 ~ dnorm(0.0,1.0E-6)
    for(j in 1:34) {
    b[j]~dnorm(0,0.01)
    }

    for(j in 1:34) {
    or[j]<-exp(b[j])
    }

    #Geostatistical random effects

    w[1:N]~spatial.exp(mu[], longnum[], latnum[], tau.sp, phi, 1)
    #N=number of clusters

    for (i in 1:N) {
    mu[i]<-0
    }

    tau.sp~dgamma(2.01,1.01)

    sigma_sp<-1/tau.sp

    phi~dgamma(0.01,0.01) I(0.3948661,7160.152)
    phi.inv<-1/phi
    Range<-3/phi

    #(phi in Winbugs is rho in the notes and tau=1/sigma^2)
    }

    data
    list( M=7878, N=404, dhsclust=c(1,1,1,1,1,2,2, 2, 2, 2, 2, 4,
    4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 6, 6, 6, 6,
    6, 8, 8, 8, 8, 8, 8, 8, 8, 8,
    8, 8, 8, 8, 9, 9, 9, 9, 9, 9,
    9, 9, 9, 9, 9, 9, 9, 9, 11, 11,
    11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
    11, 11, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 14, 14, 14, 14, 14, 14, 14,
    14, 14, 14, 14, 14, 14, 16, 16, 16, 16,
    16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
    16, 16, 16, 17, 17, 17, 17, 17, 17, 17,
    17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
    17, 17, 17, 18, 18, 18, 18, 18, 18, 18,
    18, 18, 18, 18, 18, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 21, 21, 21, 21,
    21, 21, 21, 21, 21, 21, 21, 21, 21, 21,
    22, 22, 22, 22, 22, 22, 22, 22, 22, 22, …………….

    initials
    list(
    b0=1, b=c(0.01,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ,0,0,0,0,0,0,0,0,0,0,0,0,0),
    tau.sp=1, phi=100,
    w=c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…………….

  10. Hello,

    I managed to display the error messages in the way you describe above; however, I still don’t understand what those codes mean. My model is really simple (just learning):

    ### CHEM.EXP ###
    model {
      chem.exp &lt;- h20.in*N.mgl    # chemical exposure is L/d * unit/L
      h20.in ~ dlnorm(1.5,1.563)     # water intake for an adult
      N.mgl ~ dlnorm(41.742,0.0001)  # distribution for nitrate 
    }
    ## Data
    # none
    ## Inits
    list(N.mgl=1, h20.in=1)
    

    I get the error: ‘something went wrong in procedure Monitor. Summary in module SamplesMonitors’

    I have tried reducing and increasing the iterations and changing the inits and allowing it to generate them itself. The model checks out ok but this error shows when I try to look at the stats, density, etc after I update my model.

    Many thanks.

    L

    • Looking at the trap that appears when clicking the stats button, we can see: (1) the error “undefined real result”, which suggests we may have numerically extreme or impossible values (see this link); (2) the SamplesMonitors.Monitor.Summary section indicates that .mean (the mean of the samples for one of the monitored nodes) is inf; (3) we can see which node this refers to by clicking the blue diamond beside “.monitor”, then the diamond beside “.name”, and finally, the diamond beside “.string”. This will either say “N.mgl” or “chem.exp”; (4) returning to the main trap window, we can click on the diamond beside “.sample”, which will take us to an array containing the sampled values that OpenBUGS is trying to summarise. Click the black arrows to expand the array in order to see the values. There are many “inf” values in there.

      Now, for this particular error, a simpler way to debug would be to just monitor the nodes for a few samples, then click “coda” in the Sample Monitor Tool. This will show you the same array of numbers as seen above. The trap messages are really an unnecessary complication (though the “undefined real result” error is a good indication of where the problem might lie).

      The problem derives from the fact that you have a massive variance for N.mgl. You might want to double-check your numbers… the median of N.mgl is exp(41.742), which is 1.34e18. Then, since your (log) sd is 100 (i.e., sqrt(1/0.0001)), you are sampling values that lead to numerical overflow. The mean and precision of your log-normal distribution are the mean and precision of the variable’s log. Do you expect log(N.mgl) to be centred on 41.742? Do you expect log(N.mgl) to have a standard deviation of 100?

      • Thank you so much, John.

        Thank you for noticing the error in the precision. I double checked and it should be N.mgl ~ dlnorm(41.742, 5.207) – mean, precision. The range in the data is quite large (0.1 to 1121.8 mg/L) and it is strongly right skewed. So I expect the mean to be 41.742 on the lognormal distribution (not transformed) and the standard deviation to be 138.57 according to my descriptive statistics. It is having trouble sampling from this distribution but I am unsure what to do about it? The maximum is a reasonable concentration in the data but is extreme when compared to other values in the set.

        I read that the particular error I received may occur when there are no boundaries provided so I bound the distributions and it ran initially; however, once I changed the precision it provides a new error ‘update error for node <N.mgl> algorithm univariate forward updater error can not sample node too many iterations’ so I removed the boundaries on it. I tried changing the init and generating inits but it still says it cannot do it and if I can get it to it cannot provide a density without error. During this time, I freeze the other lines with hashtags to focus on the N.mgl.

        ###CHEM.EXP###
        model {
          chem.exp <- h20.in * chem.c # chemical exposure 
          h20.in <- d.rate * risk.p # drinking rate times probability 
          risk.p ~ dbern(0.92)  # probability if they think it is safe
          d.rate ~ dlnorm(1.5, 1.563)I(0,5.0) #water intake 
          N.mgl ~ dlnorm(41.742, 5.207) # distribution for N 
          chem.c <- N.mgl
        }
        ## Data
        # none
        ## Inits
        list(N.mgl=10, d.rate=0.5,risk.p=0)
        

        I have other questions and it is so helpful after literally days of struggling on my own to be able to converse on the topic. I don’t want to exhaust your helpfulness given the specific nature of this webpage. Is there a question and answer forum that you would suggest that is current and active? Or that address error issues or model development?

        Thanks again. I will continue to review these errors. I am sure they will make more sense the more I do it!

        Sincerely,

        Lorelei

      • You can try asking over at StackOverflow. But I think the issue here is that the parameters of the lognormal distribution are not what you think they are. The lognormal distribution is described by a two parameters (location and scale), but these are, respectively, the mean (mu) and standard deviation (sigma; expressed as precision in BUGS, i.e., 1/sigma^2) of the natural log of the variable, which is normally distributed. For example, LN(10, 2) refers to a lognormally distributed variable that, if logged, has a normal distribution with mean 10 and standard deviation 2. The arithmetic mean of this lognormally-distributed variable is actually exp(mu + 0.5*sigma^2). It’s pretty confusing – see the Wikipedia page about the lognormal distribution for more info.

        The bottom line is that if your data are lognormally distributed, you can approximate the mean and sd of the corresponding normal distribution by taking the mean and sd of the log of your data (note mean and sd of the log, not the log of the mean and sd). Pass those to BUGS (converting sd to precision first). The other option is to estimate the mean and sd from the data themselves with BUGS.

        Is this part of a bigger model where you’re actually estimating parameters, though, or are you just using BUGS to simulate random data from those distributions? If the latter, you’ll save yourself some pain by just doing something like the following in R… no need for BUGS.

        n <- 1000 # number of samples you want
        d.rate <- rlnorm(n, 1.5, 1.563)
        risk.p <- rbinom(n, 1, 0.92)
        N.mgl <- rlnorm(n, 41.742, 5.207)
        h20.in <- d.rate * risk.p 
        chem.exp <- h20.in * N.mgl 
        
  11. In have been following your blog for some time; I have been using Winbugs for many years now; I few days back wanted to run one of my programs and it got all sort of error messages that I don’t understand. Yesterday used open bugs, and got stuck. I reduced the code to a very simplistic one, still the same error pops up: Sorry, something went wrong in procedure Poisson Maths random. I tried to read your solution above, I could not follow it.

  12. Dear John,
    I get the error, Expected left pointing arrow <- or twiddles ~, when I
    check model in winbugs or openbugs. But I put the <- or ~ and still the
    error occurs. Please help.

    I am fitting a negative binomial and below is the code. The error points at
    part9[i] ~ b[9]*zgdp[i] and part10[i] ~ b[10]*equals(year[i], 1).

    I checked on the google for an answer and the solution I found was that “I
    normally get this error when I use = instead of <- or ~”. In my case I have
    tried <- and ~ and the error still happens.

    Please help.
    Thanks Betty

    ####Code 
    model{      
    for(i in 1:N){  #N is Number of observations in the 2 years
    matdeaths[i]~dnegbin(p[i],r) #p[i] is probability of a maternal death in    
    #district i
    p[i]<-r/(r+mu[i]) # r is overdispersion parameter
    
    #Log of mu , mu[i] is no. livebirths in area i
    
    log(mu[i]) <- log(livebirths[i]) + intercept + part1[i] + part2[i] + 
    part3[i] + part4[i] + part5[i] + part6[i] + part7[i] + part8[i] + 
    part9[i] + part10[i] + w[district_id[i]]           
    
    part1[i] <- b[1]*zprop_4ANC_ [i]
    part2[i] <- b[2]*zprop_IPT2[i] 
    part3[i] <- b[3]*znumber_TT2_[i] 
    part4[i] <- b[4]*zprop_pregHIV__ART[i] 
    part5[i] <- b[5]*zprop_skdelivery[i] 
    part6[i] <- b[6]*zprop_PNC[i] 
    part7[i] <- b[7]*zprop_PNC_6days[i] 
    part8[i] <- b[8]*zcpr zHIVprev_preg[i] 
    part9[i] ~ b[9]*zgdp[i]
    part10[i] ~ b[10]*equals(year[i], 1) # compares the effect of 2        
    #years on the outcome
    
    }
    
    #Priors
    intercept~ dnorm(0,0.01)
    intercept_or<-exp(intercept)
    
    for (j in 1:10){
    b[j]~dnorm(0,0.01) 
    or[j]<-exp(b[j])
    }
    
    #Overdispersion parameter
    inv_r~dgamma(2.01,1.01)
    r<-1/inv_r 
    
    ######### Spike and slab prior 
    #continuous variables 
    for (j in 1:9) {
    b[j]~dnorm(0,tau1[j])
    tau1[j]<-Indc[j]*tau[j] + (1-Indc[j])*u0*tau[j]
    tau[j]~dgamma(2.01,1.01)
    Indc[j]~dbern(pind[j])
    pind[j]~dbeta(1,1)
        }
    
    #categorical variables 
    #year
    b[10]~dnorm(0,tau2[10])
    tau2[10]<-Indyr*tau[10] + (1-Indyr)*u0*tau[10]
    tau[10]~dgamma(2.01,1.01)
    pind1~dbeta(1,1)
    Indyr~dbern(pind1)
    
    #Geostatistical random effects 
    w[1:M]~spatial.exp(mu1[], longnum[], latnum[], tau.sp, phi, 1)
    
    for (i in 1:M) {
    mu1[i]<-0
    }
    
    tau.sp~dgamma(2.01,1.01)
    sigma_sp<-1/tau.sp
    phi~dunif(0.4610308, 39.28829)
    phi.inv<-1/phi
    range<-3/phi
    
    }
    
      • Dear *John Baumgartner,*

        *Thank you very much.*

        *My model works well now.*

        *I am very grateful.*

        *Thanks again.*

        *Regards*

        *Betty*

        On Mon, Oct 9, 2017 at 6:07 AM, John Baumgartner’s Research wrote:

        > John Baumgartner commented: “Hi Betty. The line part8[i] zHIVprev_preg[i] seems to be missing an operator. Hopefully that’s your > only issue. Cheers John” >

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s