How to plot survival probabilities of cox regression with firth's penalized likelihood? I am using the coxphf() function (in the coxphf package) in R and am not sure how to plot survival probabilities. Does anyone have any advice or experience with this?

If that function and package don't directly return predicted probabilities, then you can estimate the "baseline" cumulative hazard $\hat \Lambda_0(t)$ from your data and the model regression coefficients, then use the relationship $\hat S_0(t)=\exp -\hat \Lambda_0(t)$ to get the "baseline" survival curve. Those can be modified to provide probabilities for any particular set of covariate values and the corresponding hazard ratio.
Several answers on this site show how. For example, here and here and here.
You might run into problems with convergence if you have continuous covariate values like age that are far from their usual reference values of 0. Software for Cox models typically chooses a set of "baseline" covariate values for such calculations at their means or medians, or perhaps the most frequent levels of categorical predictors, to get around that problem. You might need to adapt your predictor coding similarly.
