You have a serious problem with several issues:

1 P-value selection works very poorly.

2 ordinarily I would recommend adaptive lasso for you, but that is not yet available for firth regression.

you might look at this link:

https://stats.stackexchange.com/questions/88398/model-selection-with-firth-logistic-regression

If you post your question there someone may be aware of work in development.

If not I would recommend the AIC alternative. For more on this problem in standard logistic regression see the attached. In addition this link may be helpful:

https://www4.stat.ncsu.edu/~boos/var.select/lasso.adaptive.html

Good luck, David Booth

David Eugene Booth thanks for your response, I really appreciate it.

I mentioned the variables had a P

Traditional p-value methods don't work. We cover that with citations in our paper that I attached. At this point I would consider running all interesting subsets of your IVs and and using an AIC criterion. You would have to do the regressions yourself. No package is available for firth using this method either. If you run traditional p-value selection the result is often not repeatable. See our paper and the Austin and Tu paper which I believe is available here on ResearchGate.

Good luck, David Booth

As your response variables are dichotomous, I suggest fitting both of multi-logistic and multi-nominal regressions and then using AIC criterion to find the best subsets.