Hi,

 

I am having a bit of trouble estimating a 3 category random effects logistic regression model in WinBugs.   I started with SAS, using NLMIXED on artificial data.  NLMIXED works perfectly.  However, when I try to estimate an identically specified model (or so I believe) in WinBugs, I get highly auto-correlated chains, poor convergence, and quite different estimates for model parameters.   I have pasted first the SAS code that generates the data and estimates the model.  Then I pasted the R code that reads the data and passes the data to Winbugs.  The actual model in Winbugs appears last.  I am new to Winbugs, so any offered help would be greatly appreciated.  I am particularly in the dark about specifying a prior for the precision matrix of the random effects. 

 

Thanks!

 

Sam

 

 

 

############################################################################################

 

SAS Code:

 

############################################################################################

 

 

data a;

 

do i = 1 to 100;

 

 

      b01 = .75;

      b02 = -.5;

 

   eta1 = 0;     

   eta2 = b01 + .5*rannor(0);

   eta3 = b02 + 2*rannor(0);

 

   pi1 = exp(eta1) / (exp(eta1) + exp(eta2) + exp(eta3));

   pi2 = exp(eta2) / (exp(eta1) + exp(eta2) + exp(eta3));

   pi3 = exp(eta3) / (exp(eta1) + exp(eta2) + exp(eta3));

 

 

      do j = 1 to 20;

 

      Y = rantbl(0,pi1,pi2,pi3);

 

      y1 = 0; y2 = 0; y3 = 0;

 

      if Y = 1 then y1 = 1;

      if Y = 2 then y2 = 1;

      if Y = 3 then y3 = 1;

 

      output;

      end;

end;

 

run;

 

 

 

proc nlmixed data=a(drop = b01 b02);

 

      parms s1 =  .5  s2 = 2 b01 = 0 b02 = -.5 cov12 = 0;   

 

   /*Code linear predictors*/

 

   eta1 = 0;     

   eta2 = b01 + u1;

   eta3 = b02 + u2;

 

   pi1 = exp(eta1) / (exp(eta1) + exp(eta2) + exp(eta3));

   pi2 = exp(eta2) / (exp(eta1) + exp(eta2) + exp(eta3));

   pi3 = exp(eta3) / (exp(eta1) + exp(eta2) + exp(eta3));

 

 

   /*Define likelihood*/

   z = (pi1**y1)*(pi2**y2)*(pi3**y3);

   ll = log(z);

  model y2 ~ general(ll);

 

   /*Specify random effect distribution*/

      random  u1 u2 ~ normal([0,0], [s1*s1, cov12 ,s2*s2]) subject = i;

    ods output ParameterEstimates = parms;

 

run;

 

 

 

 

PROC EXPORT DATA= WORK.A

            OUTFILE= "H:\bayes\multi.csv"

            DBMS=CSV REPLACE;

     PUTNAMES=YES;

RUN;

 

 

############################################################################################

 

R Code:

 

############################################################################################

 

 

rm(list=ls())

 

 

library(arm)

library(mvtnorm)

library(mc2d)

 

setwd('H:/bayes')

 

 

mult.dat <- read.csv('H:/bayes/multi.csv')

N <- 100

P <- 20

Ycov <- matrix(c(.5,0,0,2),nrow=2,ncol=2)

 

 

sel.dat <- list(N=N,P=P,C=3,Y = as.matrix(mult.dat[,c("y1","y2","y3")]),id = as.numeric(mult.dat[,"i"]),zeros=c(0,0),T=solve(Ycov))

 

inits1 <- list(b = c(NA,.75,-.5),R = solve(Ycov))

 

inits <- list(inits1)

 

testb <- bugs(sel.dat,inits,c("b","R"),"PS_mod3b.bug", n.chains = 1, n.iter = 50000,n.thin = 5,n.burnin = 5000, codaPkg=T,debug = T)        

 

               testb2 <- read.bugs(testb)

 

               raftery.diag(testb2)

 

summary(testb2)

 

 

 

############################################################################################

 

WinBugs Model: PS_mod3b.bug

 

############################################################################################

model

               {

                              for (i in 1:N*P) {

              

                                             Y[i, 1:C] ~ dmulti(p[i,1:C],1)

              

                                             p[i,1] <- exp(phi[i,1]) / (exp(phi[i,1]) + exp(phi[i,2]) + exp(phi[i,3]))

                              p[i,2] <- exp(phi[i,2]) / (exp(phi[i,1]) + exp(phi[i,2]) + exp(phi[i,3]))

                                             p[i,3] <- exp(phi[i,3]) / (exp(phi[i,1]) + exp(phi[i,2]) + exp(phi[i,3]))

 

                                             phi[i,1] <- b[1]

                                             phi[i,2] <- b[2] + u[id[i],1]

                                             phi[i,3] <- b[3] + u[id[i],2]             

                                            

                              }

 

                              for(j in 1:N) {

 

                              u[j,1:2] ~ dmnorm(zeros[1:2],R[ , ])

                             

                              }

                             

 

                              b[1] <- 0                             

 

                              for (j in 2:C) {

                                             b[j] ~ dnorm(0, .0001)

 

                              }                                           

 

                             

                              R[1:2,1:2] ~ dwish(T[1:2,1:2],2)

              

               }