Skip to content

Instantly share code, notes, and snippets.

@chrishanretty
Created October 27, 2012 17:50
Show Gist options
  • Save chrishanretty/3965490 to your computer and use it in GitHub Desktop.
Save chrishanretty/3965490 to your computer and use it in GitHub Desktop.
Conditional logit in R + JAGS
## 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,]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment