Ordinal regression docs #719

GStechschulte · 2023-09-13T16:48:19Z

This draft PR contains a notebook with an ordinal regression model for the newly added cumulative family.

First, it is explained why ordered categorical outcomes require special treatment. Secondly, the basics of ordinal regression are explained, motivated through the use of the cumulative link function. Next, an intercept only model is fit to explain how the cumulative link function "describes" an ordered distribution. Lastly, a model with predictors is developed.

I will also add a section with the sratio family.

Edit by Tomas Capretto: Closes #583

review-notebook-app · 2023-09-13T16:48:23Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

tomicapretto · 2023-09-15T07:31:13Z

Hi @GStechschulte! I recommend to rebase so it runs the latest version of the CI :)

Btw, is there anything I can do here to help you?

GStechschulte · 2023-09-15T07:48:15Z

Hi @GStechschulte! I recommend to rebase so it runs the latest version of the CI :)

Btw, is there anything I can do here to help you?

Hey, and will do! Thanks :) Not necessarily. Although, there is a problem with ordinal models and the interpret sub-package. Since an ordinal model prediction is a vector of probabilities of length $K$, we get shape errors when attempting to use any of the functions in interpret 😵‍💫.

I need to see if marginaleffects supports ordinal models, and see if we can adopt their solution.

tomicapretto · 2023-09-15T08:05:46Z

Hi @GStechschulte! I recommend to rebase so it runs the latest version of the CI :)
Btw, is there anything I can do here to help you?

Hey, and will do! Thanks :) Not necessarily. Although, there is a problem with ordinal models and the interpret sub-package. Since an ordinal model prediction is a vector of probabilities of length K, we get shape errors when attempting to use any of the functions in interpret 😵‍💫.

I need to see if marginaleffects supports ordinal models, and see if we can adopt their solution.

Isn't it similar to the case of the "categorical" family?

GStechschulte · 2023-09-15T08:20:30Z

@tomicapretto yes, it should be similar to the categorical family.

docs/notebooks/ordinal_regression.ipynb

GStechschulte · 2023-09-19T19:21:47Z

I think this PR is almost there. However, I am still questioning the interpretation of the threshold coefficients for the sratio model. I will paste my Slack message here.

From my understanding, for ordinal models with the sratio family and logit link, the interpretation of the threshold coefficients are not cumulative-logits, but rather logits. This is in part because the ordering of the $k$ thresholds does not matter since each outcome is characterised by its own latent distribution. Whereas ordinal models with family=cumulative and link=logit, the coefficients of the thresholds are cumulative-logits.

However, when I plot the logits of the coefficients for the sratio model, it appears as though they are “partially” cumulative (overall, the probability increases as category increases, but for some categories the probability decreases which is not possible under a cumulative specification).

Plotting the PyMC graph for the sratio model, there are no constraints.

docs/notebooks/ordinal_regression.ipynb

tomicapretto · 2023-09-21T09:24:53Z

Hi Gabriel! This is a fantastic example, it's not simple at all and you're making a great job. I left some comments with suggestions and thoughts. Also:

Can you silence FutureWarnings?
Is the warning about the usage of the C-API making the model(s) run slower?
I'm starting to think the prior for the threshold in the StoppingRatio family is not the most adequate. The means are sorted, but there's no reason to do that. So I'm thinking it makes more sense to have something like Normal([0, 0, 0..., 0]) instead of what we have now here

bambi/bambi/priors/scaler.py

Line 117 in 169564f

mu = np.round(np.linspace(-2, 2, num=response_level_n - 1), 2)

what do you think about this last point?

GStechschulte · 2023-09-22T05:35:13Z

Thank you! It has been a fun notebook to implement 😄

Can you silence FutureWarnings?

Done.

Is the warning about the usage of the C-API making the model(s) run slower?

Yes, it "feels" like it. I haven't ran any timed experiments though.

I'm starting to think the prior for the threshold in the StoppingRatio family is not the most adequate. The means are sorted, but there's no reason to do that. So I'm thinking it makes more sense to have something like Normal([0, 0, 0..., 0]) instead of what we have now here

bambi/bambi/priors/scaler.py

Line 117 in 169564f

mu = np.round(np.linspace(-2, 2, num=response_level_n - 1), 2)

Agreed, since the ordering of the thresholds doesn't matter, doing something like mu = np.zeros(response_level_n - 1) makes more sense. I will try it out. What does brms do?

Update:

Changing the default threshold prior:

response_level_n = len(attrition["YearsAtCompany"].unique())
mu = np.zeros(response_level_n - 1)
threshold_prior = {"threshold": bmb.Prior("Normal", mu=mu, sigma=1)}

sequence_model = bmb.Model(
    "YearsAtCompany ~ 0 + TotalWorkingYears", 
    data=attrition, 
    family="sratio", 
    priors=threshold_prior
)

sequence_model

       Formula: YearsAtCompany ~ 0 + TotalWorkingYears
        Family: sratio
          Link: p = logit
  Observations: 1233
        Priors: 
    target = p
        Common-level effects
            TotalWorkingYears ~ Normal(mu: 0.0, sigma: 0.3223)
        
        Auxiliary parameters
            threshold ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
             0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.], sigma: 1.0)

works as expected. For now, I am thinking I keep this and explain how to change the default prior for the thresholds. Then, we could open a separate PR? What do you think?

tomicapretto · 2023-09-22T11:06:01Z

@GStechschulte

For now, I am thinking I keep this and explain how to change the default prior for the thresholds. Then, we could open a separate PR? What do you think?

I know it's more "tidy" to do it in a separate PR, but I think it's OK to do it in this one too. Anyway, I leave it up to you. Whatever you want is OK!

tomicapretto · 2023-09-22T11:06:46Z

@GStechschulte

For now, I am thinking I keep this and explain how to change the default prior for the thresholds. Then, we could open a separate PR? What do you think?

I know it's more "tidy" to do it in a separate PR, but I think it's OK to do it in this one too. Anyway, I leave it up to you. Whatever you want is OK!

GStechschulte · 2023-09-22T12:47:15Z

I implemented the zero mu vector in this PR 😄 and added a section in the notebook explaining the differences in the default priors for cumulative and sratio.

codecov-commenter · 2023-09-22T13:33:10Z

Codecov Report

Merging #719 (5c57e55) into main (e53f8da) will not change coverage.
Report is 2 commits behind head on main.
The diff coverage is 100.00%.

❗ Current head 5c57e55 differs from pull request most recent head 075f537. Consider uploading reports for the commit 075f537 to get more accurate results

@@           Coverage Diff           @@
##             main     #719   +/-   ##
=======================================
  Coverage   89.56%   89.56%           
=======================================
  Files          44       44           
  Lines        3525     3525           
=======================================
  Hits         3157     3157           
  Misses        368      368

Files	Coverage Δ
bambi/priors/scaler.py	`96.70% <100.00%> (ø)`

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

GStechschulte · 2023-09-25T18:59:55Z

bambi/priors/scaler.py

@@ -114,7 +114,7 @@ def scale_threshold(self):
            threshold = self.model.components["threshold"]
            if isinstance(threshold, ConstantComponent) and threshold.prior.auto_scale:
                response_level_n = len(np.unique(self.response_component.response_term.data))
-                mu = np.round(np.linspace(-2, 2, num=response_level_n - 1), 2)
+                mu = np.zeros(response_level_n - 1)


If sequential models assume that for every response level there is a latent continuous variable $Z_k$, then wouldn't we need each response level? Thus, mu should be mu = response_level_n and not response_level_n - 1?

Ohh, I think you're right!

I'm thinking why the current approach is working and not failing. Is it because it's not considering the probability of Y being larger than the largest observed category? I think it would make sense for the years example, but I'm not sure if it would make sense for cases where there is a pre-specified set of categories.

I wrote that as I was looking at this visualization from the ordinal tutorial by Bürkner and Vuorre:

and I wonder: Do we always have a Pr(Y > K)? (as in the Y > 3 in the figure)

@GStechschulte if you make that modification and run the example, does it work?

If I remove the -1, I get ValueError: Incompatible Elemwise input shapes [(35,), (36,)].

This makes sense as I stated in the docs because the sequential model is a product of probabilities, i.e., the probability that $Y$ is equal to category $k$ is equal to the probability that it did not fall in one of the former categories $1: k-1$ multiplied by the probability that the sequential process stopped at $k$.

In the case of the attrition dataset, there are 36 response categories. Because of the statement above, this makes sense why the probability of category 36 is 1. There is no category after 36, so once you multiply all of the previous probabilities with the current category, you get 1. Thus, you don't need a parameter (threshold) for it.

docs/notebooks/ordinal_regression.ipynb

tomicapretto · 2023-09-25T21:37:40Z

@GStechschulte on top of the comment that you raised, there are two nits. After that, feel free to merge. Excellent work!!

edit: closed by accident haha

ordinal models (cumulative and sratio)

* zero inflated poisson and hurdle poisson models * grammar fix and sort imports * interpret coeff. and model comparison section * code review changes * change wording in hurdle Poisson section * change posterior predictive bins to use np.arange

ordinal models (cumulative and sratio)

GStechschulte · 2023-09-28T14:18:07Z

Git got me good on this one 😢

tomicapretto · 2023-09-28T16:15:14Z

I'm not sure if I follow what happened. Do you need help?

GStechschulte · 2023-09-28T18:41:52Z

I'm not sure if I follow what happened. Do you need help?

No haha. My local branch somehow diverged from this remote branch and I attempted to fix the conflicts manually. At the end, the easiest was to force push. By the way, thanks for the reviews!

* ordinal model with cumulative link notebook * ordinal model with cumulative link function ordinal models (cumulative and sratio) * unified explanation for cumulative and sequential models * sratio model and data * code review changes * remove intercept in models * zero mu vector prior for sratio family * code review and add section on default priors * explicit explanation of K and k and added summary section * Zero inflated docs (bambinos#725) * zero inflated poisson and hurdle poisson models * grammar fix and sort imports * interpret coeff. and model comparison section * code review changes * change wording in hurdle Poisson section * change posterior predictive bins to use np.arange * ordinal model with cumulative link function ordinal models (cumulative and sratio) * use plot_ppc_discrete for posterior predictive samples * add plots explaining the ordinal outcome of the dataset --------- Co-authored-by: Gabriel Stechschulte <[email protected]>

GStechschulte mentioned this pull request Sep 18, 2023

Ordinal models example #583

Closed

tomicapretto reviewed Sep 18, 2023

View reviewed changes

GStechschulte force-pushed the ordinal-examples branch 2 times, most recently from d06e32d to fb1724a Compare September 19, 2023 18:57

GStechschulte marked this pull request as ready for review September 20, 2023 14:02

GStechschulte requested a review from tomicapretto September 20, 2023 14:02

GStechschulte added the documentation label Sep 20, 2023

tomicapretto reviewed Sep 21, 2023

View reviewed changes

GStechschulte requested a review from tomicapretto September 22, 2023 12:49

GStechschulte force-pushed the ordinal-examples branch from 77d798f to e081f47 Compare September 25, 2023 18:55

GStechschulte commented Sep 25, 2023

View reviewed changes

tomicapretto reviewed Sep 25, 2023

View reviewed changes

docs/notebooks/ordinal_regression.ipynb Show resolved Hide resolved

docs/notebooks/ordinal_regression.ipynb Show resolved Hide resolved

tomicapretto closed this Sep 25, 2023

tomicapretto reopened this Sep 25, 2023

GStechschulte mentioned this pull request Sep 27, 2023

hsgp tests failing due to incompatible Elemwise input shapes #730

Closed

GStechschulte requested a review from tomicapretto September 28, 2023 06:42

ordinal model with cumulative link notebook

edbfa6b

GStechschulte and others added 12 commits September 28, 2023 16:15

ordinal model with cumulative link function

6d89eb4

ordinal models (cumulative and sratio)

unified explanation for cumulative and sequential models

a2167bc

sratio model and data

7f9b118

code review changes

3dffcce

remove intercept in models

f7a1826

zero mu vector prior for sratio family

e27194f

code review and add section on default priors

6f357f2

explicit explanation of K and k and added summary section

0521b01

Zero inflated docs (bambinos#725)

32e49fd

* zero inflated poisson and hurdle poisson models * grammar fix and sort imports * interpret coeff. and model comparison section * code review changes * change wording in hurdle Poisson section * change posterior predictive bins to use np.arange

ordinal model with cumulative link function

ff9c044

ordinal models (cumulative and sratio)

use plot_ppc_discrete for posterior predictive samples

d7acdeb

add plots explaining the ordinal outcome of the dataset

075f537

GStechschulte force-pushed the ordinal-examples branch from 610b9b3 to 075f537 Compare September 28, 2023 14:17

GStechschulte merged commit eaae9c5 into bambinos:main Sep 28, 2023
1 of 4 checks passed

GStechschulte deleted the ordinal-examples branch January 21, 2024 20:19

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ordinal regression docs #719

Ordinal regression docs #719

GStechschulte commented Sep 13, 2023 •

edited by tomicapretto

Loading

review-notebook-app bot commented Sep 13, 2023

tomicapretto commented Sep 15, 2023 •

edited

Loading

GStechschulte commented Sep 15, 2023 •

edited

Loading

tomicapretto commented Sep 15, 2023

GStechschulte commented Sep 15, 2023

GStechschulte commented Sep 19, 2023 •

edited

Loading

tomicapretto commented Sep 21, 2023

GStechschulte commented Sep 22, 2023 •

edited

Loading

tomicapretto commented Sep 22, 2023

tomicapretto commented Sep 22, 2023

GStechschulte commented Sep 22, 2023 •

edited

Loading

codecov-commenter commented Sep 22, 2023 •

edited

Loading

GStechschulte Sep 25, 2023 •

edited

Loading

tomicapretto Sep 25, 2023 •

edited

Loading

GStechschulte Sep 26, 2023 •

edited

Loading

tomicapretto commented Sep 25, 2023 •

edited

Loading

GStechschulte commented Sep 28, 2023

tomicapretto commented Sep 28, 2023

GStechschulte commented Sep 28, 2023 •

edited

Loading

Ordinal regression docs #719

Ordinal regression docs #719

Conversation

GStechschulte commented Sep 13, 2023 • edited by tomicapretto Loading

review-notebook-app bot commented Sep 13, 2023

tomicapretto commented Sep 15, 2023 • edited Loading

GStechschulte commented Sep 15, 2023 • edited Loading

tomicapretto commented Sep 15, 2023

GStechschulte commented Sep 15, 2023

GStechschulte commented Sep 19, 2023 • edited Loading

tomicapretto commented Sep 21, 2023

GStechschulte commented Sep 22, 2023 • edited Loading

tomicapretto commented Sep 22, 2023

tomicapretto commented Sep 22, 2023

GStechschulte commented Sep 22, 2023 • edited Loading

codecov-commenter commented Sep 22, 2023 • edited Loading

Codecov Report

GStechschulte Sep 25, 2023 • edited Loading

Choose a reason for hiding this comment

tomicapretto Sep 25, 2023 • edited Loading

Choose a reason for hiding this comment

GStechschulte Sep 26, 2023 • edited Loading

Choose a reason for hiding this comment

tomicapretto commented Sep 25, 2023 • edited Loading

GStechschulte commented Sep 28, 2023

tomicapretto commented Sep 28, 2023

GStechschulte commented Sep 28, 2023 • edited Loading

GStechschulte commented Sep 13, 2023 •

edited by tomicapretto

Loading

tomicapretto commented Sep 15, 2023 •

edited

Loading

GStechschulte commented Sep 15, 2023 •

edited

Loading

GStechschulte commented Sep 19, 2023 •

edited

Loading

GStechschulte commented Sep 22, 2023 •

edited

Loading

GStechschulte commented Sep 22, 2023 •

edited

Loading

codecov-commenter commented Sep 22, 2023 •

edited

Loading

GStechschulte Sep 25, 2023 •

edited

Loading

tomicapretto Sep 25, 2023 •

edited

Loading

GStechschulte Sep 26, 2023 •

edited

Loading

tomicapretto commented Sep 25, 2023 •

edited

Loading

GStechschulte commented Sep 28, 2023 •

edited

Loading