Ordered logistic regression

[datalorax]

I started a branch for dealing with ordered logistic regression via MASS::polr. I'm heading out of town through the end of the work week and won't have time to work on it more, but it might be worth taking a look. It works okay overall but the big issue currently is that the numerator and denominator are wrapped in \text and so I can't insert code like \geq easily. I don't think this should be a terrifically difficult fix, but it's something we should address.

Here's how it looks currently

library(magrittr)
library(rio)
library(equatiomatic)

d <- import("https://stats.idre.ucla.edu/stat/data/ologit.dta",
                 setclass = "tibble") %>%
	factorize()
m <- MASS::polr(apply ~ pared + public + gpa, data = d)
extract_eq(m)

On a side note, I was really happy with how little code actually had to be written to make this extension. I think we've set ourselves up well for future development.

[andrewheiss]

Hm, yeah, we'll need to somehow break it up so that \text{} goes around just the levels: P(\text{unlikely} \geq \text{somewhat likely})

[datalorax]

Yeah, so probably something with the wrap text function? Or maybe it’s just earlier in the lhs functions and I didn’t go deep enough.

Also I think we’ll have to make the the rhs a generic function too because the intercepts need to be different for these lines and we’ll need to do that down the road anyway, so I’ll get that started as I try to address the rest of this one too.

[datalorax]

Well I fixed the travis fail but it looks like interactions aren't working so maybe this was a bit premature...

[datalorax]

In fact, worse than that, it broke the interactions for ALL models, so this should definitely not be merged in yet. I'll look into why this is happening tonight, hopefully.

Side note, we should create some tests for interactions (not that you have to do it, Andrew. Thank you for setting up the initial testing structure).

[datalorax]

Addressed by #43.

[andrewheiss]

Reopening this with a couple issues (I just made a branch to add some tests, but then ran into these two things):

1. Right now this creates a set of equations for each level of the outcome variable. Right now by default, the second equation is added to the same line as the first unless the user enables wrap = TRUE. For instance:

housing <- MASS::housing
house.plr <- MASS::polr(Sat ~ Infl, data = housing)
preview(extract_eq(house.plr))

That generates this:

$$
\log\left[ \frac { P( \text{Low} \geq \text{Medium} ) }{ 1 - P( \text{Low} \geq \text{Medium} ) } \right] = \alpha_{2} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon\\\log\left[ \frac { P( \text{Medium} \geq \text{High} ) }{ 1 - P( \text{Medium} \geq \text{High} ) } \right] = \alpha_{1} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon
$$

Which looks like this when previewed:

If we turn on wrapping, each of the equations is placed in its own aligned environment instead of one shared one:

$$
\begin{aligned}
\log\left[ \frac { P( \text{Low} \geq \text{Medium} ) }{ 1 - P( \text{Low} \geq \text{Medium} ) } \right] &= \alpha_{2} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon
\end{aligned}\\\begin{aligned}
\log\left[ \frac { P( \text{Medium} \geq \text{High} ) }{ 1 - P( \text{Medium} \geq \text{High} ) } \right] &= \alpha_{1} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon
\end{aligned}
$$

Which leads to the same inline preview issue:

Ideally the TeX should look like this:

$$
\begin{aligned}
\log\left[ \frac { P( \text{Low} \geq \text{Medium} ) }{ 1 - P( \text{Low} \geq \text{Medium} ) } \right] &= \alpha_{2} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon \\
\log\left[ \frac { P( \text{Medium} \geq \text{High} ) }{ 1 - P( \text{Medium} \geq \text{High} ) } \right] &= \alpha_{1} + \beta_{1}(\text{Infl}_{\text{Medium}}) + \beta_{2}(\text{Infl}_{\text{High}}) + \epsilon
\end{aligned}
$$

Given the longer sets of equations, we also might want to think about setting wrap = TRUE by default for extract_eq.polr() or something

[andrewheiss]

And…

2. idk if this is a MASS thing or an equatiomatic thing, but it won't display some kinds of models. For instance, I made this data frame for testing purposes:

set.seed(1234)
df <- data.frame(outcome = rep(LETTERS[1:5], 60),
                 categorical = sample(letters[1:3], 100, replace = TRUE),
                 continuous_1 = rnorm(300, 100, 1),
                 continuous_2 = rnorm(300, 50, 5))
df$outcome <- ordered(df$outcome)

model_ologit <- MASS::polr(outcome ~ ., data = df)

tex <- extract_eq(model_ologit)

But I get this error when running extract_eq:

Error in as.data.frame.default(data, optional = TRUE) : 
  cannot coerce class ‘"function"’ to a data.frame 

[datalorax]

Hmm... Okay. Thanks. This has ended up being tougher than I initially though but I'm with you on all of this.

[datalorax]

So I'm not sure why this is happening, but for 2 I get that error even when I call summary(model_ologit), which is getting called by broom::tidy. So I think that's the issue. It can be resolved if you fit the model with MASS::polr(outcome ~ ., data = df, Hess = TRUE) but I'm not sure if this is generally the case? We could embed a check if so. When you try to call summary without Hess = TRUE it tries to refit it to get the Hessian but fails. So seems like it's something to do with the function rather than anything specific to extract_eq.

I'm open to suggestions for how to address this.

[datalorax]

For the first issue, it looks like this has something to do with preview.equation, but I can't seem to get it to work. If you just knit the file the equation looks as it should.

[andrewheiss]

Cool, yeah, the second issue is fixed with Hess = TRUE. Googling the error message shows that lots of people get the same error, and it happens with summary() and broom::tidy() too, so I think we just live with it and assume users will include Hess = TRUE on their own (I've always included it in my own stuff, but I had no idea why :) ).

Maybe we just include a note in the documentation that says Hess = TRUE should be included.

[datalorax]

Yeah that sounds good. I'd vote we do that and maybe even include an example of how it won't work without Hess = TRUE. I can probably add that later tonight and dig into the preview.equation issue a bit more too.