Plotting a categorical response as a function of a continuous predictor using R I would like to plot the relationship between a binary categorical response variable and a continuous predictor to study its shape. The goal is to prep a logistic regression.
This image may clarify:

I have access to Minitab and R and would greatly appreciate any insight on how to recreate this histogram or alternatives that may do just as well.
 A: This is exploration: we should feel free to be creative and to look in many different ways at the data to develop insight.
In this spirit, an attractive approach eschews binning the independent variable.  Instead, compute and smooth a running summary of the dependent variable (proportion of incomes less than 50,000 per annum).  Choose suitable windows for the summary and smooth depending on how detailed or how general a picture of the relationship is needed.
Here is an example of a synthetic dataset designed to look like the illustration:

The top row presents univariate summaries of the independent variable (hours) and dependent variable (income in thousands of dollars per year).  The latter makes it evident that information about high incomes will be relatively uncertain in this example.
The bottom row presents bivariate information.  On the left is the plot of income against hours.  (In many circumstances it would be better to model this relationship rather than splitting the income into just two groups.  But sometimes we really do just have a binary dependent variable or our analytical objective truly concerns comparing the two groups.  Let's proceed...)
The bottom right illustrates the suggested solution:


*

*The wiggly blue line is the smoothed running mean of proportions of income below 50K against hours per week.

*The surrounding gray lines are separated from the blue line by one standard error of estimate.

*The red line--which in practice would not be available--is the underlying relationship used to synthesize these data.  Ideally, it would lie entirely within (or at least close to) the region enclosed by the gray lines.
The appearance of this last plot can be controlled by varying the moving window size and by applying more or less amounts of smoothing to the moving summaries.  Here are some variations from the default width of $16$:

The smaller width (upper left) provides too much detail.  The larger widths give simpler representations; at the lower right, the relationship is nearly reduced to a straight line.  In practice--were this a real dataset (and the reference [red] line unavailable), one might decide to start with a simple linear term in hours alone and test whether introducing (say) a small cubic spline would improve the model, in effect comparing the depiction in the lower right corner to that in the upper right corner.
If you wish to preserve the ability to compare models formally--that is, to trust the p-values--then it is essential to hold out some data before conducting the exploration and test the final model against the held-out data.  If it holds up, then the model can finally be fit using all the data in order to improve the estimates of its coefficients.
R code follows.
#
# Synthesize a data set.
#
set.seed(17)
n <- 300                                # Amount of data
means <- c(15, 44) / 168                # Typical hours per week
sds <- c(5, 15) / 168                   # Dispersion around those values
ab <- means * (1-means) / sds^2 - 1     # Corresponding Beta parameters
alphas <- means * ab
betas <- ab - alphas
hours <- sort(rbeta(n, alphas, betas) * 168) # Generate a mixture of Betas

par(mfrow=c(2,2))
hist(hours)

f <- function(h, m1=-2, m2=0.3) {       # Prescribe the income-hour relationship
  x <- h/100
  0.4 + m2*x + m1*(pmin(x,0.5)-0.5) -(m1+m2)*(pmin(x,0.3)-0.3)
}

# Incomes are lognormally distributed conditional on hours
sd <- 0.4                               # CV of incomes (geometric SD)
income <- exp(rnorm(n, mean=qnorm(1-f(hours))*sd + log(50), sd=sd))
hist(income)

plot(hours, income, xlab="Hours per week", ylab="Income (K$)",
     main="Income vs. Hours") # $
#
# Compute moving summaries.
#
require(zoo)
width <- floor(sqrt(n))                        # Size of moving window
smooth.width <- min(n, 3*width) / n            # Strength of the smoother
fill <- list("extend", "extend", "extend")
x.window <- rollmean(zoo(hours), width, fill=fill)
y.window <- rollmean(zoo(income <= 50), width, fill=fill)
y.window <- zoo(lowess(y.window, f=smooth.width)$y) # $
plot(c(min(x.window), max(x.window)), c(0,1), type="n", bg="#f8f8f8",
     xlab="Hours per week", ylab="Proportion Income <= 50K",
     main="Proportion Below 50K vs. Hours")
curve(f(x), add=TRUE, col="Red")
lines(x.window, y.window, col="Blue")
lines(x.window, y.window + sqrt(y.window * (1-y.window) / width), col="Gray")
lines(x.window, y.window - sqrt(y.window * (1-y.window) / width), col="Gray")

A: The plot you highlight in your question reminds of using a loess (or lowess) curve to visualise a continuous variables against a binary response:-

Of course, the line corresponds with the histogram example at where the two colours meet.  I can't see in the your example if the data is raw or modelled (as my example is). 
This SO question shows another example.
