Multinomial logistic regression pull request

[svteichman]

This is an edited version of PR #10. I couldn't figure out how to push directly to that PR, so I'm creating a new one here that takes the code from PR #10 from @shirmath and updates based on @adw96's comments on PR #10. This implements robust score tests for multinomial logistic regression, without a Firth penalty.

[adw96]

This is looking great, @svteichman !

Can you tell me why the following doesn't give me a p-value?

multinom_test(X = structure(c(0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1), dim = c(6L, 2L)), 
              Y = structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                              0, 1, 1, 2, 1, 1, 2), dim = c(6L, 4L)), 
              strong=TRUE)

Thanks so much!

(Update: I just added this as an test that currently fails -- thought it'd be more helpful)

[adw96]

This eg may shed some light...

set.seed(1) ## runs
multinom_test(X = structure(c(0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1), dim = c(6L, 2L)), 
              Y = rpois(24, 1) + structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                              0, 1, 1, 2, 1, 1, 2), dim = c(6L, 4L)), 
              strong=TRUE)
set.seed(11) # fails
multinom_test(X = structure(c(0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1), dim = c(6L, 2L)), 
              Y = rpois(24, 1) + structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                             0, 1, 1, 2, 1, 1, 2), dim = c(6L, 4L)), 
              strong=TRUE)

And when it runs, why is the test stat always 6...?

[svteichman]

The p-value is NA when the inverse of the matrix used in the generalized score test stat is computationally singular, so it looks like it's hitting a try-catch block, muting the error from solve(), and returning an NA test stat and p-value. I'll look into the case where it works and is equal to 6.

[adw96]

Hmmmmm... do you think it should return a warning? Silently failing isn't good... Are we sure that matrix is truly singular in this case? (intuitively this wouldn't be too surprising since it is noiseless...)

Do you think Firth regularization fix singularity issues? It's not good to have a test that just sometimes can't be computed, unless we can really clearly characterize why...

[svteichman]

Yes, I'll have it report the matrix inverse error as a warning. I'll look into whether the Firth penalty will fix singularity issues.

[svteichman]

I checked your example and the matrix H2 %*% solve(I2) %*% D2 %*% solve(I2) %*% t(H2), which is the inner part of the robust score test, is indeed singular.

[codecov]

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered ☂️

[adw96]

should multinomial be included in the function title? Could be good as we plan to expand the package

[adw96]

(this may apply to other functions as well)