Plots with error bars using logistf object

Question

I initially did a logistic model using the glm package but wanted to correct for separation so I have used the logistf function and I'm now trying to redo my plots. I'm unsure how to make a plot like the one below with a logistf object. A lot of packages don't seem to support it, I've tried using sjPlot packages' plot_model() function which plots a dot for the predicted probability but doesn't add the error bars as it does automatically with a glm object. How can I get round this? Is there perhaps another package that would make this easier or is there a way to manually add the error bars?

The code for the plot I wish to add error bars to is:

sjPlot::plot_model(lr3, type="int", mdrt.values = "meansd", show.values = TRUE, value.offset = .3)

The output of my model lr3 is:

logistf(formula = foodbank_cv ~ wave + ff_country + relevel(race_grp, 
    ref = "White") + sex_cv + age_r + relevel(numchildren, 
    ref = "None") + wave * ff_hcondhas + relevel(carer, 
    ref = "Not") + sempderived + wave * cd_ff_furlough + 
    log(ff_hours) + qual + num + relevel(keyworksector, ref = "Not keyworker") + 
    ca_clinvuln_dv + freemeals + ca_blbenefits1 + log(hhincome_week), 
    data = data, firth = TRUE, family = binomial(link = "logit"))

Model fitted by Penalized ML
Coefficients:
                                                                         coef   se(coef)  lower 0.95  upper 0.95        Chisq            p method
(Intercept)                                                      -5.237542354 0.46736532 -6.23016284 -4.30807241          Inf 0.000000e+00      2
wave5                                                            -0.377956413 0.32598420 -1.07410577  0.28545651 1.232122e+00 2.669947e-01      2
wave7                                                            -0.929934987 0.40813067 -1.84652632 -0.12926473 5.260388e+00 2.181615e-02      2
ff_country2                                                      -0.118780142 0.33317501 -0.86893024  0.51197342 1.196576e-01 7.294061e-01      2
ff_country3                                                       0.393456771 0.25097814 -0.15010616  0.88210537 2.077828e+00 1.494527e-01      2
ff_country4                                                      -0.219066153 0.43493435 -1.23008781  0.57774984 2.481153e-01 6.184053e-01      2
relevel(race_grp, ref = "White")Asian or Asian British            0.882833792 0.22906054  0.39628625  1.33641305 1.183859e+01 5.801581e-04      2
relevel(race_grp, ref = "White")Black or Black British            1.759374627 0.27942672  1.16321835  2.29702048 2.678592e+01 2.272869e-07      2
relevel(race_grp, ref = "White")Mixed                             1.786978145 0.27773294  1.19285979  2.32350705 2.763841e+01 1.462461e-07      2
relevel(race_grp, ref = "White")Other                            -0.345106379 1.38712570 -5.19048868  1.62733736 6.509258e-02 7.986208e-01      2
ff_hcondhas                                                       0.691244774 0.26776923  0.14697164  1.25269746 6.228205e+00 1.257311e-02      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

The code that I used to make the hunger and race plot. I did some manual editing to make it look nicer but this is what I ideally want my plot to look like:

plot_model(model12, type = "pred", terms = c("race_grp"), mdrt.values = "meansd", axis.textsize = .3, wrap.labels = 5)+ theme_sjplot2() + scale_color_sjplot("simply") + ggplot2::labs(title= "Predicted probabilities of Hunger", x= "Race", y="Percentage")

Answer 1

I have found a way to get around this issue, however, not with the logistf package. In case anyone in the future wants to know the answer to this question, my suggestion is that you use the brglm package. I have checked and the results from the brglm package are exactly the same as the logistf package. This is how I recreated the Hunger plot posted above:

hi2<- brglm(formula= hungry_cv~ wave + ff_country + race_grp + sex_cv + age_r + numchildren + wave*ff_hcondhas + carer + sempderived + wave*cd_ff_furlough + log(ff_hours) + qual + num + keyworksector + ca_clinvuln_dv + freemeals + ca_blbenefits1 + log(hhincome_week), data=data, family=binomial(logit), method = "brglm.fit", pl = TRUE)
racehunger<- plot_model(hi2, type = "pred", terms = c("race_grp"), mdrt.values = "meansd", axis.textsize = .3, wrap.labels = 5, show.values = TRUE)+ theme_sjplot2() + ggplot2::labs(title= "Predicted probabilities of Hunger", x= "Race", y="Percentage")
racehunger

png(file="racehunger.png", units="in", width=11, height=8.5, res=300)
print(racehunger)
dev.off()

The output of the code is:

I am personally very happy with the result.