Conditional logit in R + JAGS

cl-jags.R

## Load libraries
library(mclogit)
library(reshape)
library(rjags)
library(R2WinBUGS) ## for write.model

## Load the data file from mclogit
data(Transport)

## Do the mclogit model
summary(mclogit.mod <- mclogit(
	cbind(resp,suburb)~distance+cost,
	data=Transport
))

## Now some reshaping
## We need to create matrices for the DV
## and for each IV.

## reshape choices to a suburb-by-mode data frame with values as frequencies
transport.c <- cast(Transport,suburb~transport,value="resp")

## Get predictors in same format
cost.c <- cast(Transport,suburb~transport,value="cost")
distance.c <- cast(Transport,suburb~transport,value="distance")

## Write our JAGS model
model <- function() {

	for (the.suburb in 1:nSuburbs) {
		for (the.mode in 1:nModes ) {
			log(phi[the.suburb,the.mode]) <- alpha[1] * cost[the.suburb,the.mode] + alpha[2] * distance[the.suburb,the.mode]  
			p[the.suburb,the.mode] <- phi[the.suburb,the.mode] / sum(phi[the.suburb,1:nModes])
		}
		y[the.suburb,1:nModes] ~ dmulti(p[the.suburb,1:nModes], nResp[the.suburb])
	}

	alpha[1] ~ dnorm(0,.01)
	alpha[2] ~ dnorm(0,.01)
}

write.model(model, "model1.bug")

## Set up data
y <- as.matrix(transport.c[,2:ncol(transport.c)])
nSuburbs <- nrow(y)
nModes <- ncol(y)
nResp <- rowSums(y)
cost <- as.matrix(cost.c[,2:ncol(cost.c)])
distance <- as.matrix(distance.c[,2:ncol(cost.c)])

## Write out list
data <- list("y"=y,
	"nResp"=nResp,
	"cost"=cost,
	"distance"=distance,
	"nSuburbs"=nSuburbs,
	"nModes"=nModes)

jags <- jags.model('model1.bug',
                   data = data,
                   n.chains = 3,
                   n.adapt = 100)
 
system.time(update(jags, 1e4))

system.time(out <- coda.samples(jags,
             c('alpha', 'p'),
             1e6,thin=1000))

holder <- summary(out)
gelman.diag(out)
## This fails, but see
## http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/28cef6e5 
## for an explanation

geweke.diag(out)
## Looks like we achieved convergence

## Compare the two sets of estimated coefficients
summary(mclogit.mod)
holder$statistics[1:2,]