I am running a logistic regression model where the outcome is relatively rare (250 out of 5000). My main interest is to see if there are differences between age groups, sex, educational levels, income levels...with regards to the outcome and i am running a model with just the main effects. It appears that since the model has so few "positive" outcomes, these are not fitted well by the model. This becomes apparent since all the standardized residuals for these cases are between 2-3.1 and that all these are seen as influential points based on Cook's distance (see graph), based on the common rule of thumb 4/N(=0.0008). This is not surprising since the highest predicted probability from the model is p=0.28 for one of the "positive" outcomes.
My question is if such a model is trustworthy. I am only interested in the obtained odds ratios, not how well the model can discriminate between positive/negative outcomes. I have used the Le Cessie-van Houwelingen-Copas-Hosmer global goodness of fit test to evaluate the model which gives a p-value of p-value=0.3, so from that prespective looks good but on the other hand this is expected since the vast majority of the sample is well predicted. I have also tried Firth's logistic regression to correct for bias and both regressions give almost identical results. I have finally checked what happens if i remove the 10 more influential points from the model which resulted in almost identical estimated coefficients as the model with all data.
My question is if one can use a model like this for inference, for example to identify differences between age groups and being able to rely on these estimates.
I should also add that the model has 6 categorical predictors which leads to 20 model parameters. One variable is binary, 2 variables have four levels, 2 variables have 5 levels and one has 6 levels. No separation problems were observed and no multicollinearity based on VIF.