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)
}
Samuel H. Field
Statistician Investigator
Frank Porter Graham Child Development Institute
University of North Carolina - Chapel Hill
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
