Is a binomial logistic regression valid in this case, and how do I use it / interpret its results?

Question

I am facing the unusual problem that my $p$ values are too good. They are so good that I must be doing something wrong, but I don't know what.

I am working with natural language data from a text corpus. The language I am working with has forms that are, in some contexts, optionally marked for plural number. The equivalent in English would be to be able to say both two trees and two tree. I want to test the hypothesis that marking (-s or none) depends on case (nominative [subject], accusative [object], or genitive [possessive]). In particular, I predict that in the nominative case there is relatively more -s than in the other cases.

There may be a mediating effect of the type of noun (denoting humans, animals, or things), so that the effect is stronger in some types than in others. It may also be the case that there are individual differences between nouns (thus, the ratio of marked : unmarked may be different for tree vs. goat, but still the prediction would be that that ratio is higher for nominative case than for other cases).

To test this hypothesis I have collected a number of nouns of interest, so I have a table that looks like this:

Lexeme	Type	Case	Nmarked	Nunmarked
tree	other	nominative	234	123
tree	other	genitive	456	567
tree	other	accusative	678	789
goat	animal	nominative	901	12
...	...	...	...	...

I never learned any statistics, but from looking around I think this problem calls for a binomial logistic regression. So I have tried to fit such a model on my dataset using Python statsmodels:

import pandas as pd
import statsmodels.api as sm
exog_dummies = pd.get_dummies(exog)
exog_dummies = sm.add_constant(exog_dummies, prepend=False)
glm_binom = sm.GLM(endog, exog_dummies, family=sm.families.Binomial())
res = glm_binom.fit(use_t=True)
print(res.summary())

I am using get_dummies() to convert the "Lexeme", "Type", and "Case" columns to series of columns (e.g. "Case_nominative", "Case_genitive", "Case_accusative"). To be honest I am not exactly sure what add_constant() is for, but most online examples include it and my question stands regardless of it.

If I look at the data myself, it is easy to eyeball that while there is a difference between nominative and genitive, there is very little difference between nominative and accusative. Based on this I would expect that in the model only the variable for the genitive case has a significant effect.

What I find instead is that the variables for all cases have a significant effect, with a p-value of 0.000:

                          coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
Case_nominative        -0.0962      0.012     -8.330      0.000      -0.119      -0.073
Case_genitive          -0.4009      0.015    -25.916      0.000      -0.432      -0.370
Case_accusative        -0.1058      0.021     -5.138      0.000      -0.147      -0.065

(My program also reports a pseudo $R^2$ (CS) of $1$, and a Pearson $\chi^2$ of $488$. My full dataset has $18$ lexemes with $5$ cases each, i.e. $90$ observations/rows in the table, with a total of $335,464$ word forms from the corpus.)

I do see that the coefficient is much stronger for the genitive. However, I find it dodgy that the p-values are so extremely low, and I am not sure how to report on these results to discuss the hypothesis that there is a significant difference in marking between cases.

If anyone can point out where I'm going wrong, that would be much appreciated.

Answer 1

First, I think the $p = .000$ just isn't being precise enough. You likely have a very large decimal value that is cutoff by the output here (I'm not sure how this works in Python so I am making assumptions). It would certainly be suspect if it was exactly zero. Based off what I think you say above, you have an absurdly high number of observations of around $n = 335,464$ (which is a good thing of course). Looking at your standard errors, they are consequently quite small. So whatever your p-value actually is, the output at least at surface value seems logical, but perhaps not as precise as it should be. In any case, if they are indeed considerably low and this is for a scientific article, you can simply report it as $p < .001$ (as far as I know, this is the APA standard...one reports to the nearest three decimals possible).

Having said all of that, note here that your $\beta$ coefficients here are, if set to the defaults in logistic regression, on the logit scale. To get an idea of what the actual effect size is for each predictor (how strong the effect is), you should probably transform the coefficients by exponentiating them to get odds ratios that are directly interpretable (rather than relying on $p$ values alone). Or you can transform them to probabilities and generate predictions based off different inputs.

A couple more things. First:

To be honest I am not exactly sure what add_constant() is for, but most online examples include it and my question stands regardless of it.

Never use statistical functions if you don't know what they do. Check the documentation, but I think Python is sometimes known for having bad docs on stats topics. Second:

I never learned any statistics, but from looking around I think this problem calls for a binomial logistic regression.

You really ought to stop with your analysis and learn the basics of statistics before running models. Since you appear to be doing analysis on linguistics data, I recommend reading the book Statistics for Linguists by Bodo Winter. Though it is written in the R programming language, it does an excellent job of explaining regression modeling strategies that I think would be helpful in your case.

Answer 2

Shawn's answer (+1) covered most of what I would have said. Just a little bit of elaboration follows.

Fitting the model and getting the coefficients is only the first step. That doesn't by itself answer this part of your question:

I am not sure how to report on these results to discuss the hypothesis that there is a significant difference in marking between cases.

For that you need to perform post-modeling tests that evaluate differences among the coefficients for the levels of Case (including the level subsumed in the Intercept for the model, which presumably was included when you called add_constant()). That evaluation needs to go beyond "statistical significance"; as Shawn pointed out, it's important to evaluate the odds or probabilities of differences in outcome as a function of the predictor values, and then use your understanding of the subject matter to decide on the practical significance or the differences.

I don't know what tools are available for that in statsmodels. My sense is that the documentation for statsmodels tends to assume that you already know what you're doing in terms of statistics. You will be better off using a statistical package that has accessible explanations of the assumptions in different modeling choices, like those provided by this UCLA web page. Of those, R has the advantage of being available without financial cost.

In R, the emmeans package provides such tools for many different types of regression models. The rms package provides a combination of regression-modeling functions with useful ways to present their results, with a corresponding online text.

I strongly recommend that you consult with a local statistician before you try to publish these results. From what you describe, it seems that you have made an implicit assumption that the effects of Case are independent of the effects of Type. As you have a lot of observations, you might be able to relax (and test) that assumption by including an interaction between Case and Type in your model. Interpreting the coefficients in an interaction model can get tricky, however, so it would be best to work with someone who can help you go through the details and assumptions of the model and how to present your results.