**Wiekvoet**, and kindly contributed to R-bloggers)

Last week is had a look at the standard R routines for estimating models for ordinal data. This week, I want to have a look at JAGS for examining the same data. To be honest, most of it is taking an example (inhaler) and removing code. To my surprise this was almost as difficult as adding complexity in models. I did add some output from JAGS. The data, as before, is the cheese data.

### Data preparation and model

The data is same as last week. Hence I will pick up at a created dataframe cheese2. One new item is added, a dataframe difp which is used to index differences between cheese means. The model is a simplified version of inhaler. Note that loads of pre processing of data is removed from JAGS example and placed into R. On the other hand, to minimize the amount of output, some post processing is added, hence not done in R. I guess that’s a matter of taste.

Other than that is a fairly straightforward JAGS run with the usual preparation steps.

ordmodel <- function() {

for (i in 1:N) {

for (j in 1:Ncut) {

# Cumulative probability of worse response than j

logit(QQ[i,j]) <- -(a[j] – mu[cheese[i]] )

}

p[i,1] <- 1 – QQ[i,1];

for (j in 2:Ncut) {

p[i,j] <- QQ[i,j-1] – QQ[i,j]

}

p[i,(Ncut+1)] <- QQ[i,Ncut]

response[i] ~ dcat(p[i,1:(Ncut+1)])

}

# Fixed effects

for (g in 2:G) {

# logistic mean for group

mu[g] ~ dnorm(0, 1.0E-06)

}

mu[1] <- 0

# ordered cut points for underlying continuous latent variable

a <- sort(a0)

for(i in 1:(Ncut)) {

a0[i] ~ dnorm(0, 1.0E-6)

}

#top 2 box

for (g in 1:G) {

logit(ttb[g]) <- -(a[Ncut-1] – mu[g] )

}

# Difference between categories

for (i in 1:ndifp) {

dif[i] <- mu[difpa[i]]-mu[difpb[i]]

}

}

difp <- expand.grid(a=1:4,b=1:4)

difp <- difp[difp$a>difp$b,]

datain <- list(response=cheese2$Response, N=nrow(cheese2),

cheese=c(1:4)[cheese2$Cheese],G=nlevels(cheese2$Cheese),

Ncut=nlevels(cheese2$FResponse)-1,difpa=difp$a,difpb=difp$b,ndifp =nrow(difp))

params <- c(‘mu’,’ttb’,’dif’)

inits <- function() {

list(mu = c(NA,rnorm(3,0,1)),

a0= rnorm(8,0,1))

}

jagsfit <- jags(datain, model=ordmodel, inits=inits, parameters=params,

progress.bar=”gui”,n.iter=10000)

### Results

Legend for output:

jagsfit

**leave a comment**for the author, please follow the link and comment on his blog:

**Wiekvoet**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...