# Anyone developing multinomial logistic regression?

**URL:** <https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222>\
**Category:** Optimization (Mathematical)\
**Tags:** proposal\
**Created:** [April 16, 2019, 9:39pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222 "2019-04-16T21:39:58Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [April 16, 2019, 9:39pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/1 "2019-04-16T21:39:59Z")

</div>

Hi all,

I’m looking to write a multinomial logistic regression function (including conditional logit). From what I can see, there is no extant package that can do this. Given the popularity of such a model, I’m surprised that it doesn’t exist yet in Julia, but would also be interested in collaborating. A good model for such a package would be the R [mlogit](https://cran.r-project.org/web/packages/mlogit/index.html) package.

@Nosferican @mcreel

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [April 16, 2019, 10:16pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/2 "2019-04-16T22:16:58Z")

</div>

Not sure what mlogit is but I have one in my package, Haven’t fully tested it yet I’d love feedback

> **[GitHub - caseykneale/ChemometricsTools.jl: A collection of tools for...](https://github.com/caseykneale/ChemometricsTools.jl)**
>
> A collection of tools for chemometrics and machine learning written in Julia. - GitHub - caseykneale/ChemometricsTools.jl: A collection of tools for chemometrics and machine learning written in Julia.

I also have classification measures for multiclass classifications, etc.

If there’s anything missing you’d like to see let me know and I probably have the wherewithall to add it.

Edit - Ah sorry I have a multinomial softmax regression. Sorry thats a different bird.

---

<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:** [April 17, 2019, 12:08am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/3 "2019-04-17T00:08:47Z")

</div>

I will be registering my package in the next few days which has ologit and mlogit. The implementation for mlogit is similar to Stata’s mlogit rather than the R package. The R package allows for outcome specific features (e.g., price of y0, price of y1, etc) while Stata’s and mine only handle features for the unit of observation (e.g., age, education, etc.).

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [April 17, 2019, 1:52am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/4 "2019-04-17T01:52:49Z")

</div>

Awesome! That’s definitely a great help. I am basically looking for the equivalent of Stata’s `asclogit` command but I suspect that I can implement that using your package as a starting point.

I have implemented this in Matlab but will be curious to see how you handle the optimization.

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [April 17, 2019, 1:55am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/5 "2019-04-17T01:55:14Z")

</div>

I think softmax regression is the same as multinomial logistic regression. This package is really cool and I will take a look! Thanks for responding.

---

<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:** [April 17, 2019, 2:57pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/6 "2019-04-17T14:57:51Z")

</div>

For nominal response models, it would be nice to add the generalization for conditional multinomial logistic (PR welcome). For ordinal logistic regression, I implemented the proportional odds logistic regression (POLR), but the generalization would also be welcomed as a PR. I am waiting for StatsModels to tag a new release for the Terms 2.0 era and for me to finish up the last final touches.

In terms of the optimization approach, (1) the random effects (Swamy-Arora) is fitted through a partial demeaning framework, (2) absorption of categorical fixed effects is done trough the method of alternating projections, (3) multinomial logistic regression is optimized through Fisher Scoring in a Vector Generalized Linear Model (VGLM) framework, and (4) ordinal logistic regression uses Optim.jl with a Newtonian with AD for the Hessian because I lost so much time trying to write up the analytical Hessian and failed.

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [April 17, 2019, 3:11pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/7 "2019-04-17T15:11:28Z")

</div>

Awesome, thanks! My Matlab implementation of Stata’s `asclogit` command is basically the equivalent of calling `Optim.jl`. I suspect Fisher scoring is faster, but I’ve never programmed that up before.

I had a peek at your source code this morning and will be in touch with a possible PR in the coming weeks.

---

<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:** [April 17, 2019, 3:22pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/8 "2019-04-17T15:22:11Z")

</div>

It might still need a tweak or two (e.g., StatsModels changes in master)… I want to get the beta test out in the coming weeks and then start working on optimizing the code and more robustness checks. That was a chapter in my defense (two weeks ago) so probably would be good to release it soon and have it stable by JuliaCon.

---

<div class="post-metadata">

**Author:** ![infinitylabs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/infinitylabs/32/10563_2.png) [@infinitylabs](https://discourse.julialang.org/u/infinitylabs)\
**Post date:** [October 1, 2019, 9:31am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/10 "2019-10-01T09:31:38Z")

</div>

Can you provide the link to your package?

---

<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:** [October 1, 2019, 1:16pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/11 "2019-10-01T13:16:55Z")

</div>

[https://github.com/Nosferican/Econometrics.jl](https://github.com/Nosferican/Econometrics.jl)

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [October 2, 2019, 6:52am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/12 "2019-10-02T06:52:16Z")

</div>

For those interested in regularised regression models including multinomial regressions, I’ve just released [MLJLinearModels.jl](https://github.com/alan-turing-institute/MLJLinearModels.jl) which will soon be integrated with the MLJ environment; it’s competitive with Sklearn and similar R packages though rough around the edge and without docs (which is why it hasn’t been officially announced yet)

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 2, 2019, 10:57am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/13 "2019-10-02T10:57:07Z")

</div>

@tlienart - When this matures I’ll likely interface my package with it. Looks like a great start.

---

<div class="post-metadata">

**Author:** ![alequa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alequa/32/12338_2.png) [@alequa](https://discourse.julialang.org/u/alequa)\
**Post date:** [April 3, 2020, 12:51pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/14 "2020-04-03T12:51:44Z")

</div>

Hi,  
I am trying to use your package, but I really have no clue how!  
Is there any documentation or example I can follow?

Best,  
Alessio

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [April 3, 2020, 1:09pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/15 "2020-04-03T13:09:50Z")

</div>

Yeah this is my fault really… there’s a docs branch but since starting a new job I’m struggling a bit to find the time to push it through

The tests provide quite good examples including comparisons with sklearn: [MLJLinearModels.jl/logistic-multinomial.jl at master · JuliaAI/MLJLinearModels.jl · GitHub](https://github.com/alan-turing-institute/MLJLinearModels.jl/blob/master/test/fit/logistic-multinomial.jl)

Long story short if you have a n x p matrix X and targets y encoded as 1,2,3,… (indices of the class) then you can just do the following:

```julia
mnr = MultinomialRegression(0.5) # 0.5 is the strength of the L2 regularisation
theta = fit(mnr, X, y)

```

not very dissimilar to how sklearn is called (and actually you can see the comparison code [MLJLinearModels.jl/logistic-multinomial.jl at b5d91ac9d13bb98f78eda9b9f347e8a6ddb39aaa · JuliaAI/MLJLinearModels.jl · GitHub](https://github.com/alan-turing-institute/MLJLinearModels.jl/blob/b5d91ac9d13bb98f78eda9b9f347e8a6ddb39aaa/test/fit/logistic-multinomial.jl#L90-L101))

To do the predictions you can apply X to theta then use a softmax:

```julia
preds = apply_X(Xmatrix, θ, c) # c is the number of classes
preds .= softmax(preds)

```

for binary, same process but with ±1 labels and using sigmoid instead of softmax.

If you go the MLJ route, you can use the fit/predict machinery, see [this example](https://alan-turing-institute.github.io/MLJTutorials/end-to-end/wine/index.html)

---

<div class="post-metadata">

**Author:** ![alequa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alequa/32/12338_2.png) [@alequa](https://discourse.julialang.org/u/alequa)\
**Post date:** [April 3, 2020, 3:37pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/16 "2020-04-03T15:37:22Z")

</div>

Ok, Cool.  
Indeed I looked at the tests, but still I was quite lost!

And how the matrix theta is structured?

```julia
X = rand(100,3)
y = rand([1,2,3],100)
λ = 0.5
mnr = MultinomialRegression(λ; fit_intercept=true)
θ = fit(mnr, X, y)
X = rand(100,3)
y = rand([1,2,3],100)

parameters = reshape(θ,(3,4))

```

where parameters[1, 1:3] are the coefficient of the first class and parameters[1, 4] is the intercept?  
Maybe I am dumb and not understanding how Logit works, but I am referring to this:

> **[Multinomial logistic regression | Linear predictor](https://en.wikipedia.org/wiki/Multinomial_logistic_regression#Linear_predictor)**
>
> As in other forms of linear regression, multinomial logistic regression uses a linear predictor function 
>   
>     
>       
> f
> (
> k
> ,
> i
> )
>       
>     
> {\\displaystyle f(k,i)}
>   
> to predict the probability that observation i has outcome k, of the following form:

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [April 3, 2020, 4:06pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/17 "2020-04-03T16:06:35Z")

</div>

So I think the way the prediction work should help you there, basically (assuming no intercept first) what it does is:

1. apply\_X, this basically means (X \* theta) = z
2. apply a softmax on z

If you check the dimensions you will see that X is (n x p ) where n is the number of records, p is the dimensionality (number of features), and theta effectively corresponds to a (p x c) matrix where c is the number of classes (it’s stored as a single vector though, which is why you see this reshape; each column is one after the other)

Now if s = softmax z, it’s n x c (like z) with the rows summing to 1, basically each row contains the score for each class. For instance for 3 classes you would get predictions with rows like

[0.2, 0.5, 0.3]

meaning that the second class is the “most likely” for that point (most highly scored according to the model)

I hope this clarifies it a bit!

**Edit** : more specifically, in the no intercept case, each class gets `p` parameters, these are stacked one after the other in theta so `1:p` then `p+1:2p` etc and theta has total length `p*c` (change `p` for `p+1` for the case with intercept) if you want to compare it with what sklearn gives, you indeed have to do `reshape(theta, (p, c))` or `reshape(theta, (p+1, c))` (with intercept)

**Edit2** : code for apply might be useful too: [MLJLinearModels.jl/utils.jl at b5d91ac9d13bb98f78eda9b9f347e8a6ddb39aaa · JuliaAI/MLJLinearModels.jl · GitHub](https://github.com/alan-turing-institute/MLJLinearModels.jl/blob/b5d91ac9d13bb98f78eda9b9f347e8a6ddb39aaa/src/utils.jl#L35-L46)

---

<div class="post-metadata">

**Author:** ![alequa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alequa/32/12338_2.png) [@alequa](https://discourse.julialang.org/u/alequa)\
**Post date:** [April 3, 2020, 4:23pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/18 "2020-04-03T16:23:28Z")

</div>

Thank you!  
The Mnist database got classified pretty well 🙂

There was just something ‘odd’

I downloaded MNIST from

```julia
using Flux, Flux.Data.MNIST, Statistics
imgs = MNIST.images()
X = Array(transpose(hcat(float.(reshape.(imgs, :))...) ))
labels = MNIST.labels() .+1

```

The labels go from 0 to 9, but I had to translate them to [1,10] otherwise I got this error.  
Looking at the number I had the intuition it was scoring over 1:9, like the class with 0 was ‘not welcome’

For the rest it worked smoothly,

Thanks a lot

```julia
DimensionMismatch("new dimensions (785, 10) must be consistent with array size 7056")
(::Base.var"#throw_dmrsa#197")(::Tuple{Int64,Int64}, ::Int64) at reshapedarray.jl:41
reshape at reshapedarray.jl:45 [inlined]
reshape at reshapedarray.jl:116 [inlined]
apply_X!(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64) at utils.jl:66
(::MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64})(::Float64, ::Array{Float64,1}, ::Array{Float64,1}) at d_logistic.jl:149
(::NLSolversBase.var"#61#62"{NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing},Float64})(::Array{Float64,1}, ::Array{Float64,1}) at incomplete.jl:45
value_gradient!!(::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at interface.jl:82
initial_state(::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}, ::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at l_bfgs.jl:158
optimize(::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at optimize.jl:33
#optimize#93 at interface.jl:116 [inlined]
optimize(::NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at interface.jl:115
optimize(::NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}) at interface.jl:115
_fit(::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::LBFGS, ::Array{Float64,2}, ::Array{Int64,1}) at newton.jl:114
#fit#144(::LBFGS, ::typeof(fit), ::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::Array{Float64,2}, ::Array{Int64,1}) at default.jl:48
fit(::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::Array{Float64,2}, ::Array{Int64,1}) at default.jl:38
top-level scope at test_LR.jl:159

```

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [April 3, 2020, 4:39pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/19 "2020-04-03T16:39:15Z")

</div>

Yes that’s actually a requirements to have the encoding to be 1…c (partly related to if you have a categorical vector and want to encode it to integer)

Cool thanks for reporting this! Would you mind opening an issue with your reproducing code? I might add it as an example in these much needed docs 😅

---

<div class="post-metadata">

**Author:** ![alequa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alequa/32/12338_2.png) [@alequa](https://discourse.julialang.org/u/alequa)\
**Post date:** [April 3, 2020, 5:07pm UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/20 "2020-04-03T17:07:49Z")

</div>

Here you go.

[https://github.com/alan-turing-institute/MLJLinearModels.jl/issues/58](https://github.com/alan-turing-institute/MLJLinearModels.jl/issues/58)

I would say that this topic got solved 🙂  
Thanks a lot for the help  
Alessio

The example code:

```julia
# Classify MNIST digits with a simple multi-layer-perceptron
using Flux.Data.MNIST
using MLJLinearModels
# Get MNIST dataset and transpose for (records, features)
imgs = MNIST.images()
X = Array(transpose(hcat(float.(reshape.(imgs, :))...) )
# MNIST labels: Categorical labels must be 1...c, hence add .+1 to each label
labels = MNIST.labels() .+1
# and the number of classes
n_classes = length(Set(labels))
n_features = size(X,2)
# The MNIST database does not need the intercept
intercept = false
# deploy MultinomialRegression from MLJLinearModels, λ being the strenght of the reguliser
mnr = MultinomialRegression(λ; fit_intercept=intercept)
# Fit the model
θ = fit(mnr, X, labels)
# The model parameters are organized such we can apply X⋅θ, the following is only to clarify
params = reshape(θ, n_features +Int(intercept), n_classes)
# Get the predictions X⋅θ 
preds = MLJLinearModels.softmax(MLJLinearModels.apply_X(X,θ,n_classes))
# map each vector to its maximal element 
targets = map(x->argmax(x),eachrow(preds))
#and evaluate the model over the labels
scores = 1- sum(targets-labels)/length(preds)

```

---

<div class="post-metadata">

**Author:** ![Sourish](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sourish/32/34156_2.png) [@Sourish](https://discourse.julialang.org/u/Sourish)\
**Post date:** [October 6, 2022, 9:30am UTC](https://discourse.julialang.org/t/anyone-developing-multinomial-logistic-regression/23222/21 "2022-10-06T09:30:02Z")

</div>

Here is an implementation of **[Bayesian Multinomial Logistic Regression](https://turing.ml/dev/tutorials/08-multinomial-logistic-regression/)**