Hi all,
I wanted to let you know that over Christmas vacation I solved this problem.
I bought Congdon's _Bayesian Statistical Modelling_ and applied the example
code he had in Chapter 7, I think it was 7.1. The main difference, it seems,
was that all other MNL code I've seen uses y ~ dmulti, though Congdon's code
uses y ~ dcat. Another thing I like about Congdon's approach is that you
don't have to modify data set from wide to tall, as with the other code I found.
Below is a working version, and below that is the original request. Thank you
all.
############################ Directory Paths
############################
# Home
MyModelPath <- "c:/R/WinBUGS/"
# Work
#MyModelPath <- "c:/Documents and Settings/My Documents/WinBUGS/"
# Universal
MyBUGSPath <- "c:/Program Files/WinBUGS14/"
MyModelFile <- paste(MyModelPath, "model.bug", sep="")
###################### Randomly Create a Data Set
#######################
y <- c(1:400); x02 <- c(1:400); x03 <- c(1:400)
y[1:100] <- 1
y[101:200] <- 2
y[201:300] <- 3
y[301:400] <- 4
x02[1:100] <- rnorm(100, 25000, 2500)
x02[101:200] <- rnorm(100, 40000, 4000)
x02[201:300] <- rnorm(100, 35000, 3500)
x02[301:400] <- rnorm(100, 30000, 3000)
x03[1:100] <- rnorm(100, 2.51, 0.25)
x03[101:200] <- rnorm(100, 2.01, 0.20)
x03[201:300] <- rnorm(100, 2.70, 0.27)
x03[301:400] <- rnorm(100, 1.66, 0.16)
MyData <- cbind(y, x02, x03)
MyData <- as.data.frame(MyData)
y.n <- NROW(y)
y.j <- length(unique(y))
############### Center and Scale any Continuous Predictors
##############
MyData$x02.c <- (MyData$x02 - mean(MyData$x02)) / (2 * sd(MyData$x02))
MyData$x03.c <- (MyData$x03 - mean(MyData$x03)) / (2 * sd(MyData$x03))
################## Frequentist Multinomial Logit Model
##################
library(nnet)
MyFMNL <- multinom(y ~ x02.c + x03.c, data=MyData)
MyFMNL
MyProbs <- predict(MyFMNL, type="probs")
MyData$FPROB_1 <- MyProbs[,1]
MyData$FPROB_2 <- MyProbs[,2]
MyData$FPROB_3 <- MyProbs[,3]
MyData$FPROB_4 <- MyProbs[,4]
MyData$F_YHAT <- predict(MyFMNL)
with(MyData, table(y, F_YHAT))
######################## Format Data for WinBUGS
########################
MyBUGSData <- list(y=MyData$y,
x02=MyData$x02.c, x03=MyData$x03.c,
y.n=y.n, y.j=y.j)
########################## WinBUGS Model File
###########################
# This model code came from Congdon's BSM textbook (2006), example 7.2
library(R2WinBUGS)
model <- function()
{
### Likelihood Definition
for (i in 1:y.n)
{
### Random Component
y[i] ~ dcat(p[i,1:y.j])
### Deterministic Component and Link Function
for (j in 1:y.j)
{
p[i,j] <- phi[i,j] / sum(phi[i,])
log(phi[i,j]) <- b[j,1] + b[j,2]*x02[i] +
b[j,3]*x03[i]
}
}
### Priors (form = b[j, n.var])
b[1,1] <- 0 #Reference category, set to zero
b[1,2] <- 0
b[1,3] <- 0
b[2,1] ~ dnorm(b.2.01.mu, b.2.01.tau)
b[2,2] ~ dnorm(b.2.02.mu, b.2.02.tau)
b[2,3] ~ dnorm(b.2.03.mu, b.2.03.tau)
b[3,1] ~ dnorm(b.3.01.mu, b.3.01.tau)
b[3,2] ~ dnorm(b.3.02.mu, b.3.02.tau)
b[3,3] ~ dnorm(b.3.03.mu, b.3.03.tau)
b[4,1] ~ dnorm(b.4.01.mu, b.4.01.tau)
b[4,2] ~ dnorm(b.4.02.mu, b.4.02.tau)
b[4,3] ~ dnorm(b.4.03.mu, b.4.03.tau)
### Hyperpriors
b.2.01.mu ~ dnorm(0.0, 1.0E-2)
b.2.02.mu ~ dnorm(0.0, 1.0E-2)
b.2.03.mu ~ dnorm(0.0, 1.0E-2)
b.3.01.mu ~ dnorm(0.0, 1.0E-2)
b.3.02.mu ~ dnorm(0.0, 1.0E-2)
b.3.03.mu ~ dnorm(0.0, 1.0E-2)
b.4.01.mu ~ dnorm(0.0, 1.0E-2)
b.4.02.mu ~ dnorm(0.0, 1.0E-2)
b.4.03.mu ~ dnorm(0.0, 1.0E-2)
b.2.01.tau <- pow(b.2.01.sigma, -2)
b.2.02.tau <- pow(b.2.02.sigma, -2)
b.2.03.tau <- pow(b.2.03.sigma, -2)
b.3.01.tau <- pow(b.3.01.sigma, -2)
b.3.02.tau <- pow(b.3.02.sigma, -2)
b.3.03.tau <- pow(b.3.03.sigma, -2)
b.4.01.tau <- pow(b.4.01.sigma, -2)
b.4.02.tau <- pow(b.4.02.sigma, -2)
b.4.03.tau <- pow(b.4.03.sigma, -2)
b.2.01.sigma ~ dunif(0, 100)
b.2.02.sigma ~ dunif(0, 100)
b.2.03.sigma ~ dunif(0, 100)
b.3.01.sigma ~ dunif(0, 100)
b.3.02.sigma ~ dunif(0, 100)
b.3.03.sigma ~ dunif(0, 100)
b.4.01.sigma ~ dunif(0, 100)
b.4.02.sigma ~ dunif(0, 100)
b.4.03.sigma ~ dunif(0, 100)
}
write.model(model, MyModelFile)
#file.show(MyModelFile)
############################# WinBUGS Model
#############################
system.time(
MyModel <- bugs(MyBUGSData, inits=NULL,
model.file=MyModelFile,
parameters.to.save=c("b"),
n.chains=3, n.iter=2000, n.burnin=1000, n.thin=2, codaPkg=TRUE,
bugs.directory = MyBUGSPath, working.directory=MyModelPath,
debug=TRUE)
# Saved draws = ((n.iter - n.burnin) / n.thin) * n.chains
)
############################ MCMC Diagnostics
###########################
library(coda)
MyHBMNL <- read.bugs(MyModel)
summary(MyHBMNL)
gelman.diag(MyHBMNL, transform=TRUE) #Increase iterations if not below 1.10
par(ask=TRUE); gelman.plot(MyHBMNL, ylim=c(1,1.2))
par(ask=TRUE); plot(MyHBMNL, density=FALSE)
par(ask=TRUE); plot(MyHBMNL, trace=FALSE, col="red")
par(ask=TRUE); autocorr.plot(MyHBMNL, col="red")
effectiveSize(MyHBMNL) #Autocorrelation effectively reduced draws to this
par(ask=TRUE); cumuplot(MyHBMNL, col="red") #Cumulative quantile plots
########################## Posterior Data Set
###########################
Chain1 <- MyHBMNL[[1]]; Chain1 <- as.data.frame(Chain1)
Chain2 <- MyHBMNL[[2]]; Chain2 <- as.data.frame(Chain2)
Chain3 <- MyHBMNL[[3]]; Chain3 <- as.data.frame(Chain3)
PosteriorDraws <- rbind(Chain1, Chain2, Chain3)
rm(Chain1, Chain2, Chain3)
#edit(PosteriorDraws)
########################## Posterior Inference
##########################
par(ask=TRUE); boxplot(PosteriorDraws[,1:9], col="red", xlab="Parameters",
ylab="Posterior Densities"); abline(h=0, col="blue")
mean(PosteriorDraws[,2] > 0)
mean(PosteriorDraws[,2] > PosteriorDraws[,3])
############################ Score Data Set
#############################
b2_01 <- mean(PosteriorDraws[,1]); b2_02 <- mean(PosteriorDraws[,2])
b2_03 <- mean(PosteriorDraws[,3]); b3_01 <- mean(PosteriorDraws[,4])
b3_02 <- mean(PosteriorDraws[,5]); b3_03 <- mean(PosteriorDraws[,6])
b4_01 <- mean(PosteriorDraws[,7]); b4_02 <- mean(PosteriorDraws[,8])
b4_03 <- mean(PosteriorDraws[,9])
XB1 <- 0
XB2 <- (b2_01 + b2_02*MyData$x02.c + b2_03*MyData$x03.c)
XB3 <- (b3_01 + b3_02*MyData$x02.c + b3_03*MyData$x03.c)
XB4 <- (b4_01 + b4_02*MyData$x02.c + b4_03*MyData$x03.c)
MyData$phat1 <- (exp(XB1) / (exp(XB1) + exp(XB2) + exp(XB3) + exp(XB4)))
MyData$phat2 <- (exp(XB2) / (exp(XB1) + exp(XB2) + exp(XB3) + exp(XB4)))
MyData$phat3 <- (exp(XB3) / (exp(XB1) + exp(XB2) + exp(XB3) + exp(XB4)))
MyData$phat4 <- (exp(XB4) / (exp(XB1) + exp(XB2) + exp(XB3) + exp(XB4)))
rm(XB1, XB2, XB3, XB4)
MyData$yhat <- 1
MyData$yhat <- ifelse(MyData$phat2 > MyData$phat1 &
MyData$phat2 > MyData$phat3 & MyData$phat2 > MyData$phat4,
2, MyData$yhat)
MyData$yhat <- ifelse(MyData$phat3 > MyData$phat1 &
MyData$phat3 > MyData$phat2 & MyData$phat3 > MyData$phat4,
3, MyData$yhat)
MyData$yhat <- ifelse(MyData$phat4 > MyData$phat1 &
MyData$phat4 > MyData$phat2 & MyData$phat4 > MyData$phat3,
4, MyData$yhat)
Classification_Table <- with(MyData, table(y, yhat))
Pred_Accuracy <- 1:4
Pred_Accuracy[1] <- (Classification_Table[1,1] /
sum(Classification_Table[1,1:4]))
Pred_Accuracy[2] <- (Classification_Table[2,2] /
sum(Classification_Table[2,1:4]))
Pred_Accuracy[3] <- (Classification_Table[3,3] /
sum(Classification_Table[3,1:4]))
Pred_Accuracy[4] <- (Classification_Table[4,4] /
sum(Classification_Table[4,1:4]))
Pred_Accuracy_Total <- ((Classification_Table[1,1]
+ Classification_Table[2,2] + Classification_Table[3,3]
+ Classification_Table[4,4])
/ sum(Classification_Table))
Balanced_Lift <- (Pred_Accuracy_Total / 0.25) - 1
Classification_Table
Pred_Accuracy
Pred_Accuracy_Total
Balanced_Lift
>
> I am new to R, WinBUGS, and email lists, but ask for your help. I am
> using the R2WinBUGS package in R to use WinBUGS. I am building small,
> simple versions of models, trying to get them right, and then build in
complexity, of course.
> I've had great success using WinBUGS to replicate a frequentist binary
> logit model in R. Now, I'm trying to do the same thing with a
> multinomial logit, but my results are very different. I am convinced
> that the problem must be in the model specification.
>
> The data set is randomly generated (so I'm obviously just learning,
> not trying to solve a real problem with it), and the continuous
> predictors have been centered and scaled (techinically not
> standardized) according to Gelman's
> (2007) book.
>
> Here is an example of how the WinBUGS MNL model is not predicting the
> same as the frequentist MNL model from R:
>
> Frequentist MNL using multinom() in R
> Classification Table
> yhat
> y 1 2 3
> 1 80 1 19
> 2 0 98 2
> 3 23 1 76
>
> WinBUGS MNL
> Classification Table
>
> PS: Although I am only using 3,000 iterations, the problem is the same
> even with 50,000.
>
> Here is the entire code:
>
> ############################ Directory Paths
> ############################ # Home MyModelPath <-
> "c:/Compound/R/WinBUGS/"
> # Work
> MyModelPath <- "c:/Documents and Settings/My Documents/WinBUGS/"
> # Universal
> MyBUGSPath <- "c:/Program Files/WinBUGS14/"
> MyModelFile <- paste(MyModelPath, "model.bug", sep="")
>
> ###################### Randomly Create a Data Set
> ####################### y <- c(1:300); x02 <- c(1:300); x03 <-
> c(1:300); x04 <- c(1:300) y[1:100] <- 1 y[101:200] <- 2 y[201:300] <-
> 3 x02[1:100] <- round(runif(100, 0, 1)) x02[101:200] <-
> round(runif(100, 0.3, 1)) x02[201:300] <- round(runif(100, -0.5, 1))
> x03[1:100] <- rnorm(100,-40,10) x03[101:200] <- rnorm(100,10,10)
> x03[201:300] <- rnorm(100,-30,10) x04[1:100] <- rnorm(100,0,2)
> x04[101:200] <- rnorm(100,1,2) x04[201:300] <- rnorm(100,2,2) MyData
> <- cbind(y, x02, x03, x04) MyData <- as.data.frame(MyData)
>
> ############### Center and Scale any Continuous Predictors
> ############## MyData$x03.c <- (MyData$x03 - mean
(MyData$x03)) / (2 *
> sd(MyData$x03)) MyData$x04.c <- (MyData$x04 - mean(MyData$x04)) / (2
*
> sd(MyData$x04))
>
> ################## Frequentist Multinomial Logit Model
> ##################
> library(nnet)
> MyFMNL <- multinom(y ~ x02 + x03.c + x04.c
> , data=MyData)
> MyFMNL
> MyProbs <- predict(MyFMNL, type="probs")
> MyData$FPROB_1 <- MyProbs[,1]
> MyData$FPROB_2 <- MyProbs[,2]
> MyData$FPROB_3 <- MyProbs[,3]
> MyData$F_YHAT <- predict(MyFMNL)
> with(MyData, table(y, F_YHAT))
>
> ##################### Convert DV from Wide to Tall
> ###################### y <- MyData$y y.new <-
> as.numeric(as.vector(outer(sort(unique(y)), y, "=="))) y.k <-
> length(unique(y)) #Num. of Categories of y y.n <- NROW(y)
>
> ###################### Format Data for WinBUGS
> ########################## MyBUGSData <-
> list(y=structure(.Data=y.new,.Dim=c(y.n,y.k)),
> x02=MyData$x02, x03=MyData$x03.c, x04=MyData$x04.c,
> numberofobservations=y.n, numberofchoices=y.k)
>
> ############ WinBUGS Bayesian Multinomial Logit Model File
> ##############
> library(R2WinBUGS)
> model <- function()
> {
> ### Likelihood Definition
> for (i in 1:numberofobservations)
> {
> ### Random Component
> y[i, 1:numberofchoices] ~ dmulti(p[i, 1:numberofchoices], 1)
> ### Deterministic Component and Link Function
> for (j in 1:numberofchoices)
> {
> p[i,j] <- phi[i,j] / sum(phi[i,])
> log(phi[i,j]) <- (b[j,1] + b[j,2]*x02[i] + b[j,3]*x03[i] +
> b[j,4]*x04[i])
> }
> }
> ### Priors: (b[DVcat,ParamNum])
> # The first DV category is the reference category, and is set to zero.
> # Every other prior is assigned a normal distribution.
> b[1,1] <- 0
> b[2,1] ~ dnorm(0, .01)
> b[3,1] ~ dnorm(0, .01)
> b[1,2] <- 0
> b[2,2] ~ dnorm(0, .01)
> b[3,2] ~ dnorm(0, .01)
> b[1,3] <- 0
> b[2,3] ~ dnorm(0, .01)
> b[3,3] ~ dnorm(0, .01)
> b[1,4] <- 0
> b[2,4] ~ dnorm(0, .01)
> b[3,4] ~ dnorm(0, .01)
> }
> write.model(model, MyModelFile)
> #file.show(MyModelFile)
>
> ############################# WinBUGS Model
> ############################# system.time( MyModel <- bugs
(MyBUGSData,
> inits=NULL,
> model.file=MyModelFile,
> parameters.to.save=c("b"),
> n.chains=3, n.iter=3000, n.burnin=1500, n.thin=5, codaPkg=TRUE,
> bugs.directory = MyBUGSPath, working.directory=MyModelPath,
> debug=TRUE)
> # Saved draws = ((n.iter - n.burnin) / n.thin) * n.chains
> )
> ############################ MCMC Diagnostics
> ###########################
> library(coda)
> MyHBMNL <- read.bugs(MyModel)
> summary(MyHBMNL)
> gelman.diag(MyHBMNL, transform=TRUE) #Increase iterations if not below
> 1.10 par(ask=TRUE); gelman.plot(MyHBMNL, ylim=c(1,1.2)) par(ask=TRUE);
> plot(MyHBMNL, density=FALSE) par(ask=TRUE); plot(MyHBMNL,
trace=FALSE,
> col="red") par(ask=TRUE); autocorr.plot(MyHBMNL, col="red")
> effectiveSize(MyHBMNL) #Autocorrelation effectively reduced draws to
> this par(ask=TRUE); cumuplot(MyHBMNL, col="red") #Cumulative quantile
> plots
>
> ########################## Posterior Data Set
> ###########################
> Chain1 <- MyHBMNL[[1]]; Chain1 <- as.data.frame(Chain1)
> Chain2 <- MyHBMNL[[2]]; Chain2 <- as.data.frame(Chain2)
> Chain3 <- MyHBMNL[[3]]; Chain3 <- as.data.frame(Chain3) PosteriorDraws
> <- rbind(Chain1, Chain2, Chain3) rm(Chain1, Chain2, Chain3)
> #edit(PosteriorDraws)
>
> ########################## Posterior Inference
> ########################## par(ask=TRUE);
> boxplot(PosteriorDraws[,1:8], col="red", xlab="Parameters",
> ylab="Posterior Densities"); abline(h=0, col="blue")
> mean(PosteriorDraws[,2] > 0)
> mean(PosteriorDraws[,2] > PosteriorDraws[,3])
>
> ############################ Score Data Set
> #############################
> b02_1 <- mean(PosteriorDraws[,1]); b02_2 <- mean(PosteriorDraws[,2])
> b02_3 <- mean(PosteriorDraws[,3]); b02_4 <- mean(PosteriorDraws[,4])
> b03_1 <- mean(PosteriorDraws[,5]); b03_2 <- mean(PosteriorDraws[,6])
> b03_3 <- mean(PosteriorDraws[,7]); b03_4 <- mean(PosteriorDraws[,8])
> MyData$phat1 <- (exp(0) / (exp(0) +
> exp(b02_1 + b02_2*MyData$x02 + b02_3*MyData$x03.c +
> b02_4*MyData$x04.c) +
> exp(b03_1 + b03_2*MyData$x02 + b03_3*MyData$x03.c +
> b03_4*MyData$x04.c)))
> MyData$phat2 <- (exp(b02_1 + b02_2*MyData$x02 + b02_3*MyData$x03.c
+
> b02_4*MyData$x04.c) / (exp(0) +
> exp(b02_1 + b02_2*MyData$x02 + b02_3*MyData$x03.c +
> b02_4*MyData$x04.c) +
> exp(b03_1 + b03_2*MyData$x02 + b03_3*MyData$x03.c +
> b03_4*MyData$x04.c)))
> MyData$phat3 <- (exp(b03_1 + b03_2*MyData$x02 + b03_3*MyData$x03.c
+
> b03_4*MyData$x04.c) / (exp(0) +
> exp(b02_1 + b02_2*MyData$x02 + b02_3*MyData$x03.c +
> b02_4*MyData$x04.c) +
> exp(b03_1 + b03_2*MyData$x02 + b03_3*MyData$x03.c +
> b03_4*MyData$x04.c)))
> MyData$yhat <- 1
> MyData$yhat <- ifelse(MyData$phat2 > MyData$phat1 &
> MyData$phat2 > MyData$phat3, 2, MyData$yhat) MyData$yhat <-
> ifelse(MyData$phat3 > MyData$phat1 &
> MyData$phat3 > MyData$phat2, 3, MyData$yhat)
Classification_Table <-
> with(MyData, table(y, yhat)) Pred_Accuracy <- 1:3 Pred_Accuracy[1] <-
> (Classification_Table[1,1] /
> sum(Classification_Table[1,1:3]))
> Pred_Accuracy[2] <- (Classification_Table[2,2] /
> sum(Classification_Table[2,1:3]))
> Pred_Accuracy[3] <- (Classification_Table[3,3] /
> sum(Classification_Table[3,1:3]))
> Pred_Accuracy_Total <- ((Classification_Table[1,1]
> + Classification_Table[2,2] + Classification_Table[3,3])
> / sum(Classification_Table))
> Balanced_Lift <- (Pred_Accuracy_Total / 0.33) - 1 Classification_Table
> Pred_Accuracy Pred_Accuracy_Total Balanced_Lift
>
> # End
>
> Any thoughts are appreciated. Thanks.
-------------------------------------------------------------------
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
