Using the {ggplot2} package there's quite a bit you can do with model plotting. For example
library(ggplot2)
fit <- lm(mpg ~ drat, data = mtcars)
ggplot(mtcars,aes(drat,mpg)) + geom_smooth(method = "lm") + theme_minimal()
#> `geom_smooth()` using formula 'y ~ x'

However, when the y-axis is binary, it ends up looking more like
It's possible to get to something more like
but getting there isn't simple. See this explainer.