The resources that I have seen feature graphs such as the following

This is fine if the predictor $x$ is continuous, but if the predictor is categorical and just has a few levels it's not clear to me how to justify the logistic model / curve.

I have seen this post, this is not a question about whether or not binary logistic regression can be carried out using categorical predictors.

What I'm interested in is how to explain the use of the logistic curve in this model, as there doesn't seem to be a clear way like there is for a continuous predictor.

edit

data

data that has been used for this simulation

library(vcd)
set.seed(2019)
n = 1000
y = rbinom(2*n, 1, 0.6)
x = rbinom(2*n, 1, 0.6)

crosstabulation

> table(df)
 y
x 0 1
 0 293 523
 1 461 723


> prop.table(table(df))
 y
x 0 1
 0 0.1465 0.2615
 1 0.2305 0.3615

mosaic plot

First, you could make a graph like that with a categorical x. It's true that that curve would not make much sense, but ...so? You could say similar things about curves used in evaluating linear regression.

Second, you can look at crosstabulations, this is especially useful for comparing the DV to a single categorical IV (which is what your plot above does, for a continuous IV). A more graphical way to look at this is a mosaic plot.

Third, it gets more interesting when you look at multiple IVs. A mosaic plot can handle two IVs pretty easily, but they get messy with more. If there are not a great many variables or levels, you can get the predicted probablity for every combination.

Overview of Binary Logistic Regression Using a Continuous Predictor

A binary logistic regression model with continuous predictor variable $x$ has the form:

log(odds that y = 1) = beta0 + beta1 * x (A)

According to this model, the continuous predictor variable $x$ has a linear effect on the log odds that the binary response variable $y$ is equal to 1 (rather than 0).

One can easily show this model to be equivalent to the following model:

 (probability that y = 1) =

 exp(beta0 + beta1 * x)/[1 + exp(beta0 + beta1 * x)] (B)

In the equivalent model, the continuous predictor $x$ has a nonlinear effect on the probability that $y = 1$.

In the plot that you shared, the S-shaped blue curve is obtained by plotting the right hand side of equation (B) above as a function of $x$ and shows how the probability that $y = 1$ increases (nonlinearly) as the values of $x$ increase.

BinaryLogistic Regression Using a Categorical Predictor with Two Categories, Whose Effect is Encoded Using a Dummy Variable

If your x variable were a categorical predictor with, say, $2$ categories, then it would be coded via a dummy variable $x$ in your model, such that x = 0 for the first (or reference) category and $x = 1$ for the second (or non-reference) category. In that case, your binary logistic regression model would still be expressed as in equation (B). However, since $x$ is a dummy variable, the model would be simplified as:

log(odds that y = 1) = beta0 for the reference category of x (C1) 

and

log(odds that y = 1) = beta0 + beta1 for the non-reference category of x (C2)

The equations (C1) and (C2) can be further manipulated and re-expressed as:

(probability that y = 1) = exp(beta0)/[1 + exp(beta0)] for the reference category of x (D1)

and

(probability that y = 1) =
 exp(beta0 + beta1)/[1 + exp(beta0 + beta1)] for the non-reference category of x (D2)

What Is the Utility of a Binary Logistic Regression with a Categorical Predictor with Two Categories?

So what is the utility of the binary logistic regression when $x$ is a dummy variable?

The model allows you to estimate two different probabilities that $y = 1$: one for $x = 0$ (as per equation (D1)) and one for $x = 1$ (as per equation (D2)).

You could create a plot to visualize these two probabilities as a function of $x$ and superimpose the observed values of $y$ for $x = 0$ (i.e., a whole bunch of zeroes and ones sitting on top of $x = 0$) and for $x = 1$ (i.e., a whole bunch of zeroes and ones sitting on top of $x = 1$). The plot would look like this:

 ^
 |
y = 1 | 1 1 
 |
 | *
 |
 | *
 |
y = 0 | 0 0
 |
 |------------------>
 x = 0 x = 1
 x-axis

In this plot, you can see the zero values (i.e., $y = 0$) stacked atop $x = 0$ and $x = 1$, as well as the one values (i.e., $y = 1$) stacked atop $x = 0$ and $x = 1$. The * symbols denote the estimated values of the probability that $y = 1$.

There are no more curves in this plot as you are just estimating two distinct probabilities.

If you wanted to, you could connect these estimated probabilities with a straight line to indicate whether the estimated probability that $y = 1$ increases or decreases when you move from $x = 0$ to $x = 1$. Of course, you could also jitter the zeroes and ones shown in the plot to avoid plotting them right on top of each other.

What Is the Utility of a Binary Logistic Regression with a Categorical Predictor with More Than Two Categories?

If your categorical predictor variable $x$ has $k$ categories, where $k > 2$, then your model would include $k - 1$ dummy variables and could be written to make it clear that it estimates $k$ distinct probabilities that $y = 1$ (one for each category of $x$). You could visualize the estimated probabilities by extending the plot shown above to incorporate $k$ categories for $x$. For example, if $k = 3$, the plot would look like this:

 ^
 |
y = 1| 1 1 1
 | *
 | *
 |
 | *
 |
y = 0| 0 0 0
 |
 |---------------------------------->
 x = 1st x = 2nd x = 3rd
 x-axis

where 1st, 2nd and 3rd refer to the first, second and third category of the categorical predictor variable $x$.

Creating the suggested plots using R and simulated data

Note that the effects package in R will create plots similar to what I suggested here, except that the plots will NOT show the observed values of $y$ corresponding to each category of $x$ and will display uncertainty intervals (i.e., 95% confidence intervals) around the plotted (estimated) probabilities. Simply use these commands:

install.packages("effects")
library(effects)

model <- glm(y ~ x, data = data, family = "binomial") 

plot(allEffects(model))

For the data in your post, this would be:

library(effects)

set.seed(2019)

n = 1000

y = as.factor(rbinom(2*n, 1, 0.6))

x = as.factor(rbinom(2*n, 1, 0.6))

df = data.frame(x=x,y=y)

model <- glm(y ~ x, data = df, family = "binomial")

plot(allEffects(model),ylim=c(0,1))

The resulting plot can be seen at https://m.imgur.com/klwame5.