Logistic regression yields very small p-values for any/every explanatory variable in my dataset

Question

I'm investigating whether a federal district court judge's ABA rating (rating given to the judge by the American Bar Association when he/she is nominated) significantly correlates to the rate that his/her opinions are reversed on appeal. (I'm in the USA.) There are 645 judges total. My data looks like so:

judge	gender	president	president_party	aba_rating	cases_tot	aff_tot	rev_tot	rev_rate
Judge1	Male	Clinton	Democrat	Qualified	272	211	61	0.22426471
Judge2	Male	Obama	Democrat	Not Qualified	279	215	64	0.22939068
Judge3	Female	Obama	Democrat	Well Qualified	348	310	38	0.1091954
Judge4	Female	Bush II	Republican	Well Qualified	129	97	32	0.24806202
Judge5	Male	Trump	Republican	Not Qualified	6	6	0	0
...	...	...	...	...	...	...

From the very helpful Cross Validated community, I believe I should be doing logistic regression using the glm() function. I am using the following model:

model_X <- glm(cbind(aff_tot, rev_tot) ~ [explanatory variable(s)], data = dct_data, family = binomial)

(I'm using cbind(aff_tot, rev_tot) as my dependent variable instead of a simple rev_rate because I want to account for the fact that judges have had different numbers of cases appealed.)

I'm sure this is just because I'm new to R/new to stats, but no matter what explanatory variable(s) I used, I almost always get really small p-values (examples below), even when it obviously doesn't make sense to have such small p-values. Am I setting my dataset up incorrectly? Am I using the glm() function incorrectly? I've asked several stats friends, and they couldn't figure it out.

...

[If it's helpful, I'm here are some examples...]

Example 1: Gender

All judges are either Female or Male. The average reversal rate for Females is 12.2%, and the average for Males is 13.5%. But when I use the below model, I get a p-value of 2x10^-16 for the coefficient on Females and 2.5x10^-8 for Males.

model_gender <- glm(cbind(aff_tot, rev_tot) ~ gender, data = dct_data, family = binomial)

Example 2: President Party

All judges were nominated either by a Democrat or a Republican. The average reversal rate for Democrat-appointed judges is 13.5%, and the average for Republican-appointed judges is 13.0%. But when I use the below model, I get a p-value of 2x10^-16 for the coefficient on Democrats and 0.0203 for Republicans.

model_president_party <- glm(cbind(aff_tot, rev_tot) ~ president_party, data = dct_data, family = binomial)

With 645 judges, my instinct is that it's impossible that a ~1% difference in reversal rates for Female vs. Male should yield such a low p-value. (Ditto for the 0.5% reversal rate difference for Democrats vs. Republicans.) But whenever I look for a correlation between reversal rate and any of the explanatory variables, or whenever I combine multiple variables in the same model, I almost always get really low p-values for each explanatory variable. Again, I'm sure I'm making a dumb mistake, but I'd really appreciate it if someone could show me what I'm doing wrong.

Answer 1

Your instinct that something isn't quite right here is correct in my opinion. However I would caution against the direct conclusion that the p value is "too small" given the data. If the null hypothesis is true, ie. if there is really no effect of party of gender, then the p value follows a uniform distribution, so a very small p value is just as likely as a very large one.

That said, the issue here is that the sample size is not 645. It is 113,007, and your glm model is treating the data as if they are independent. With such a sample size, even a very small effect can lead to a small p value. Maybe the data are independent, however we must not rule out the possibility that the probability of a reversal for a particular decision made by one judge is more likely to be similar to the probability of a reversal for a anothe decision made by the same judge, than that of a different judge. I would therefore consider fitting either a Generalised Estimating Equations model, or a mixed effects model, with random intercepts for judge. If there is considerably intra-class correlation, then the p values will probably be higher.