Bayesian inference for rank-ordered logit model in Julia's Turing.jl

rol_model.jl

using Turing
using LogExpFunctions: logsumexp
using DataFrames

# The rank ordered logit model in Turing
# The rank-ordered logit likelihood
function rank_ordered_logit(ordered_skills::Vector{<:Real})
    ll = 0.0
    for m in 1:(length(ordered_skills) - 1)
        ll += ordered_skills[m] - logsumexp(ordered_skills[m:end])
    end
    return ll
end

# Instantiate the model
@model function rol_model(results::Vector)
    K = maximum(maximum(results))
    theta ~ MvNormal(ones(K))
    for result in results
        Turing.@addlogprob! rank_ordered_logit(theta[result])
    end
end

# here are some competition results
# each number belongs to a certain competitor
# sometimes people don't compete, that's fine
results = [
    [3, 2, 5, 4, 1, 6],
    [2, 3, 6, 5, 4, 1],
    [3, 4, 2, 6],
    [2, 4, 1, 5, 6]
]

# Sample from the posterior
chain = sample(rol_model(results), NUTS(), 5_000)

# Create a ranking with estimated abilities
theta_posterior = DataFrame(chain)[!,3:8]
df = DataFrame(
    :id  => 1:maximum(maximum(results)),
    :est => [mean(v) for v in eachcol(theta_posterior)],
    :lwr => [quantile(v, 0.025) for v in eachcol(theta_posterior)],
    :upr => [quantile(v, 0.975) for v in eachcol(theta_posterior)]
);
sort(df, :est, rev=true)

The results are reasonable 😄

6×4 DataFrame
 Row │ id     est          lwr        upr      
     │ Int64  Float64      Float64    Float64
─────┼─────────────────────────────────────────
   1 │     3   1.21856     -0.285547  2.67353
   2 │     2   0.818461    -0.545962  2.19325
   3 │     4   0.00541904  -1.25964   1.18965
   4 │     5  -0.188238    -1.53195   1.10542
   5 │     1  -0.716967    -2.05064   0.57822
   6 │     6  -1.11853     -2.52127   0.170459