Guidance on the visualization of phylolm models

[MarcRieraDominguez]

Hi! I am a frequent user of phylolm (currently I rely on version 2.6.4), and I regularly find myself stumped when it comes to visualizing the models. This is particularly so because I average PGLS models (MuMIn::dredge() + subset(delta <= threshold) + MuMIn::model.avg()).
For instance, ggeffects::gpredict() won't work with model-averaged phylolm models (strengejacke/ggeffects#387). Another example is lack of compatiblity between phlyolm() and visreg() (#64), which does not show a confidence band.
Admitedly, one can use model coefficients to plot a line (https://stackoverflow.com/questions/43441467/plot-phylogenetic-logistic-regression-with-binary-response-variable), but this does not show uncertainty around the estimate, which is particularly important when visually assessing overlap among factor levels. As far as I can tell, phylolm lacks a se.fit = TRUE option for the predict() method, which would be a very interesting feature.
Can anyone recommend a way to visualize model predictions, effects, marginal means? I am particularly interested in a way to show uncertainty around a prediction or estimated effect. If the solution also applies to model-averaged objects it would be best!
Thank you for your time, and congratulations for the great package!