#############################################################################
# Ordered Logit model
# Alan Agresti "Categorical Data
Analysis"
#
Dependent variable :
# Mental Impairment (mental,
1: well, 2:mild, 3: moderate, 4:impaired)
# Independent
variables:
# Socioeconomic status: (ses,
0: low, 1:high)
# Life event index (a measure
of the number and severity of
#
important life events within
the past three years)
#############################################################################
model
{
for(
i in 1:n)
{
mental[i]
~ dcat( p[i,] )
p[i,
1] <- 1/( 1 + exp( -c[1] + mu[i] ) )
p[i,
2] <- 1/( 1 + exp( -c[2] + mu[i] ) ) - 1/( 1 + exp( -c[1] + mu[i] ) )
p[i,
3] <- 1/( 1 + exp( -c[3] + mu[i] ) ) - 1/( 1 + exp( -c[2] + mu[i] ) )
p[i,
4] <- 1- 1/( 1 + exp( -c[3] + mu[i] ) )
mu[i] <- b[1]*ses[i] +
b[2]*life[i]/sd( life[] )
}
c[1]
<- g[1]
c[2]
<- g[1] + g[2]
c[3]
<- g[1] + g[2] + g[3]
g[1]
~ dnorm(0, 1.0E-6)
g[2]
~ dgamma(0.001, 0.001)
g[3]
~ dgamma(0.001, 0.001)
b.ses <- b[1]
b.life <- b[2]/sd( life[] )
b[1]
~ dnorm(0, 1.0E-6)
b[2]
~ dnorm(0, 1.0E-6)
# Prediction Pr( mental>=2 | ses, life, y)
for(
i in 1:9)
{
life.pred[i] <- i
mental.pred.0[i]
<- 1-1/( 1 + exp( -c[1] + b.life*i ) )
mental.pred.1[i]
<- 1-1/( 1 + exp( -c[1] + b.ses + b.life*i ) )
}
}
# Init
list( b = c(0, 0), g = c( 1, 1, 1) )
# Data 1
list( n=40)
# Data 2
mental[] ses[] life[]
1 1 1
1 1 9
1 1 4
1 1 3
1 0 2
1 1 0
1 0 1
1 1 3
1 1 3
1 1 7
1 0 1
1 0 2
2 1 5
2 0 6
2 1 3
2 0 1
2 1 8
2 1 2
2 0 5
2 1 5
2 1 9
2 0 3
2 1 3
2 1 1
3 0 0
3 1 4
3 0 3
3 0 9
3 1 6
3 0 4
3 0 3
4 1 8
4 1 2
4 1 7
4 0 5
4 0 4
4 0 4
4 1 8
4 0 8
4 0 9
END