# Replicating ologit results from Econometrics.jl

**URL:** <https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819>\
**Category:** Data\
**Created:** [September 7, 2021, 2:51pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819 "2021-09-07T14:51:30Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 7, 2021, 2:51pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/1 "2021-09-07T14:51:30Z")

</div>

I’m having a hard time replicating the results of an ordered logit regression in `Econometrics.jl`. The code for the estimation is done [here](https://github.com/Nosferican/Econometrics.jl/blob/32659ae549100488a00e4d016870fea877c915e4/src/solvers.jl#L228). The objective function is [here](https://github.com/Nosferican/Econometrics.jl/blob/32659ae549100488a00e4d016870fea877c915e4/src/solvers.jl#L257).

Here is my implementation (you can download the data from [here](https://github.com/Nosferican/Econometrics.jl/blob/master/data/auto.csv))

```julia
using ComponentArrays
using UnPack
using Econometrics
using Distributions
using Optim
using CategoricalArrays
using DataFrames
using CSV

datadf = DataFrame(CSV.File("auto.csv"))
datadf.foreign = categorical(datadf.foreign)
datadf.foreign = recode(datadf.foreign, "Domestic" => 0, "Foreign" => 1)
datadf.rep77 = categorical(datadf.rep77, ordered=true)
datadf.rep77 = recode(datadf.rep77, "Poor" => 1, "Fair" => 2, "Average" => 3, "Good" => 4, "Excellent" => 5)

datadfnm = dropmissing(datadf[!, [:rep77, :foreign, :length, :mpg]])

# Econometrics.jl results
ecn = fit(EconometricModel, @formula(rep77 ~ foreign + length + mpg), datadfnm)

# my implementation
θ = ComponentArray{Float64}(β=[0.1, 0.1, 0.1], K=(1:4)./4)

function Λ(x)
    cdf(Logistic(), x)
end

function mycumsum(K)
    n = length(K)
    s = similar(K)
    s[1] = K[1]
    for i in 2:n
        s[i] = abs(K[i]) + s[i-1]
    end
    return s
end

function Pi(j, θ, x)
    @unpack β, K = θ
    KK = mycumsum(K)
    kk = length(K)
    n = length(x)
    xb = sum(x[i]*β[i] for i in 1:n)
    if 2 <= j <= kk
        return Λ(KK[j] - xb) - Λ(KK[j-1] - xb)
    elseif j == 1
        return Λ(KK[1] - xb)
    elseif j == kk+1
        return 1 - Λ(KK[kk] - xb)
    end
end

function boole(x, pred)
    x == pred
end

function logℒ(θ; Yvec=Yvec, Xmat=Xmat)
    lev = unique(Yvec)
    s = zero(eltype(θ))
    for i in 1:size(Xmat, 1), j in 1:length(lev)
        xi = @view(Xmat[i, :])
        s += log(Pi(j, θ, xi)) * boole(Yvec[i], lev[j])
    end
    return s
end

Xmat = Matrix(datadfnm[!, [:foreign, :length, :mpg]])
Yvec = datadfnm[!, :rep77]

θ0 = similar(θ)
θ0.β .= 0.0
θ0.K .= log.(1:4)

opt = optimize(θ -> -logℒ(θ), θ, Newton(), Optim.Options(g_tol = 1e-10, iterations=10_000); autodiff=:forward)

# my coefficients
opt.minimizer.β
# my thresholds
opt.minimizer.K |> mycumsum

# Econometrics.jl coefficients and thresholds:
ecn.β

```

Can any one give some advice?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [September 7, 2021, 2:57pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/2 "2021-09-07T14:57:36Z")

</div>

I think your best bet is to summon @Nosferican. In these types of situation I usually use RCall to try and get another opinion as well.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [September 7, 2021, 3:08pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/3 "2021-09-07T15:08:46Z")

</div>

I have been summoned. 🙂  
Did you take a look at the test suite? The results are verified against Stata 13 `ologit` (should be the same or quite close to Stata 15) and R’s `MASS:polr`. Did you check the your implementation against other software?

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 7, 2021, 9:34pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/4 "2021-09-07T21:34:30Z")

</div>

Stata 16.1 and `Econometrics.jl` output the same results. I still don’t see how my implementation is different. Where in the algorithm do I do something fundamentally different than `Econometrics.jl`.

I need to run a ologit many times, hence why I’m trying to implement a barebones algorithm.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [September 7, 2021, 9:46pm UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/5 "2021-09-07T21:46:30Z")

</div>

The easiest implementation to “borrow from” would be  
[https://github.com/Nosferican/Econometrics.jl/blob/32659ae549100488a00e4d016870fea877c915e4/src/solvers.jl#L228-L295](https://github.com/Nosferican/Econometrics.jl/blob/32659ae549100488a00e4d016870fea877c915e4/src/solvers.jl#L228-L295)  
You can adapt the “barebones” routine and re-use arrays and scratch objects as needed for a speed-up as applicable. If you have some guarantees such that you don’t have to check for linearly separable cases or perfect multicollinearity or other checks, you could get rid of those to make it less “safe” but potentially faster. It depends on the specifics of your use case.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 8, 2021, 2:09am UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/6 "2021-09-08T02:09:58Z")

</div>

I finally figured it out!  
The problem was that `unique(Yvec)` does not yield a sorted vector. Once I changed the line `lev = unique(Yvec)` to `lev = levels(Yvec)` I got the same result as `Econometrics.jl` (and Stata).

Aside, I use `abs` to compute the increasing thresholds, while @nosferican you use `exp`. For values away from 0, it should be fine I think. I have used `abs` in other cases to ensure that some quantity remains positive, but is there any concern in using `abs` vs `exp`?

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [September 8, 2021, 9:51am UTC](https://discourse.julialang.org/t/replicating-ologit-results-from-econometrics-jl/67819/7 "2021-09-08T09:51:42Z")

</div>

That, I recall, the `exp` method is definitely better. I think that was a numerical trick I learned from the MASS authors.