Look into stability of brglm2, then add as an option (replacing glm.fit with brglm.fit as needed) for running glm to adding Firth penalty functionality. This will require a bit of craftiness to go from a design matrix (what raoBust uses to drop a column to estimate under the null) to a formula and data object (what brglm expects).