# Implementing an ordered probit model in pymc

I'm trying to implement an ordered probit model in pymc, and I'm stuck. The model is similar to Welinder's "multidimensional wisdom of crowds", with coders (indexed by i) and documents (indexed by j). Coders assign codes to documents, but the coding process is noisy.

We wish to estimate two things. First, $z_j$, the true, underlying value of each document along some latent dimension. Second, $\beta_i$, bias terms revealing how accurately each coder's assessments line up with the group average.

This would be pretty easy if codes were on a continuous scale, but the data I have is ordinal. Codes fall in the range 1,2,...,5 .

Formally, we can represent this as follows:

$x_{ij} \sim Normal( \alpha_i + \beta_iz_j, 1)$

$code_{ij} = cut(x_{ij}, w)$

Where $cut$ is a cutoff function that assigns values based on the index of the highest cutpoints, $w$, exceeded by $x_{ij}$.

So far, so good. The problem with this model is that the cutpoint function is deterministic, and codes are observed. But in pymc (and in other MCMC programs, e.g. JAGS), a deterministic node cannot also be observed. So this model can't be built directly in pymc.

It seems that there's probably a way to treat $x$ as deterministic, and $codes$ as a random function of $x$. This would probably make $code$ involve a Categorical node. But I'm not sure how to specify the probability function, and I'm a little worried my whole approach may be off. Can anyone set me straight?

Also -- long shot -- if there's anyone who codes in pymc, it'd be great to see source code for this. Ordinal probit/logit is a pretty standard model, but I can't find a pymc example anywhere online.

-
What you have there looks like an ordered probit, not logit, model. –  Macro May 22 '12 at 18:07
Ah. Good point. I've changed the title to "probit," since the focus of the question was really on the "ordered" part. –  Abe May 23 '12 at 16:02

This took quite a bit of work, but I got it in the end. Note that I used the development version (pymc 2.2grad) from github, not the older version available on pypi.

Also, this runs pretty slowly, since it doesn't make good use of numpy's array manipulation. For data sets of reasonable size, some smart preprocessing and a rewrite of the cutpoints node function could probably fix this.

Here's full source code:

from scipy.stats import norm, pearsonr
from pymc import Normal, Lambda, Uniform, Exponential, stochastic, deterministic, observed, MCMC, Matplot
from numpy import mean, std, log
import numpy as np

#Set array dimensions
(I, J, K, M, N) = (5, 20, 3, 4, 1000)

#Set simulation parameters
alpha_star = np.random.normal(0, 1, size=(I,K))
beta_star = np.random.normal(1, 1, size=(I,K))
z_star = np.random.normal(0, 1, size=(J,1))
w_star = np.array([0,1,3])

#Generate data
coder = np.random.randint(I, size=(N))
doc = np.random.randint(J, size=(N))
item = np.random.randint(K, size=(N))

code = np.zeros(shape=(N))
for n in range(N):
i, j, k = coder[n], doc[n], item[n]
m = alpha_star[i,k] + beta_star[i,k] * z_star[j] + np.random.normal(0, 1)
code[n] = 1+sum(m > w_star)
#    print "\t".join([str(x) for x in [i, j, k, m, code[n]]])

#Set GLM parameters
alpha = Normal('alpha', mu=0.0, tau=0.01, value=np.zeros(I*K))
beta = Normal('beta', mu=1.0, tau=0.01, value=np.ones(I*K))

z = Normal('z', mu=0.0, tau=0.01, value=np.random.normal(0,1,J))

w = Exponential('w', .1, value=np.ones(M-3))

mu = Lambda('mu', lambda alpha=alpha, beta=beta, z=z, i=coder, j=doc, k=item: alpha[i+I*k]+beta[i+I*k]*z[j])

@deterministic(plot=False)
def cutpoints(w=w):
w2 = [-np.inf, 0.0, 1.0]
v = 1
for i in w:
v += i
w2.append(v)
w2.append(np.inf)

cp = np.array( w2 )
return cp

@stochastic(dtype=int, observed=True)
def y(value=code, mu=mu, cp=cutpoints):

def logp(value, mu, cp):
d = norm.cdf(cp[value]-mu)-norm.cdf(cp[value-1]-mu)
lp = sum(log(d))
return lp

#Run chain
M = MCMC([alpha, beta, z, mu, w, cutpoints, y])
M.isample(10000, 5000, thin=5, verbose=0)

#Summarize results
Matplot.summary_plot([alpha], name="alpha", path="./graphs/")
Matplot.summary_plot([beta], name="beta", path="./graphs/")
Matplot.summary_plot([z], name="z", path="./graphs/")
Matplot.summary_plot([w], name="w", path="./graphs/")

print pearsonr( alpha_star.transpose().reshape((I*K,)), alpha.stats()['mean'])
print pearsonr( beta_star.transpose().reshape((I*K,)), beta.stats()['mean'])
print pearsonr( z_star.reshape((J,)), z.stats()['mean'])

-