Binomial GLM probability and cdplot probability for same data are wildly different - why? Code using glm:
model <- glm(vs ~ hp, data=mtcars, family=binomial)
newdata <- data.frame(hp=seq(min(mtcars$hp), max(mtcars$hp),len=100))
newdata$vs = predict(model, newdata, type="response")
plot(vs ~ hp, data=mtcars)
lines(vs ~ hp, newdata, lwd=2)
plogis(predict(model, newdata = data.frame(hp=150)))

Code using cdplot:
mtcarscdplot <- cdplot(as.factor(vs) ~ hp, data = mtcars)
mtcarscdplot$`1`(150)

The probabilities for 150 should at least be close, but from the glm it is 0.1294217 while from the cdplot it is 0.3168849 - a massive difference. Could anyone shed light on this, please? Which probability is the correct (and reportable) one and where is the other one going wrong?
 A: It helps to think about what lies behind these / how they are generated.
The conditional density plot (cdplot) is an exploratory data analysis approach for the relationship between a binary and a continuous variable.  It is not a model.  Basically, two kernel density plots are made, and you get the first as a proportion of the sum of both densities.  I illustrate the process here: Any necessary EDA before logistic?  I make them occasionally, mostly because it's quick and easy.  There are other EDA methods as well.
The logistic regression model is a statistical model.  However, your model only fits a straight line (in the log odds transformed space).  It's quite possible that a straight line is not appropriate.  We can try adding a squared term, but: 1) with N=32 data, 14 v-shaped and 18 non v-shaped engines, we just don't have much information to fit a flexible line; and 2) this still isn't as flexible as the conditional density plot.  Notice the probability is now lower!
m2 = glm(vs~poly(hp, 2), mtcars, family=binomial)
lines(50:335, predict(m2, data.frame(hp=50:335), type="response"),
      col="red")
predict(m2, data.frame(hp=150), type="response")  # 0.1041799


Assuming the straight-line fit is correct, you can get a 95% confidence interval for that predicted probability.  It definitely includes .317:
y.hat.150 = predict(model, data.frame(hp=150), type="link", se.fit=TRUE)
with(y.hat.150, plogis(c(fit-(1.96*se.fit), fit+(1.96*se.fit)))) 
# 0.01580684 0.57913339

If you're married to the CD-plot, and you're sure the kernel and default bandwidth estimation are correct, you could try to resample to get a sense of the uncertainty there.  While admittedly, I should put more into polishing this, it's... kind of a mess:
set.seed(7454)
mat = matrix(NA, nrow=512, ncol=100)
for(i in 1:100){
  dat = mtcars[sample(1:32, size=32, replace=T), c("hp", "vs")]
  dy  = with(dat, density(hp[vs==1], from=50, to=335))
  dn  = with(dat, density(hp[vs==0], from=50, to=335))
  sm  = dy$y + dn$y
  mat[,i] = dy$y/sm
}
windows()
  plot(1,1, xlim=c(50,335), ylim=c(0,1), ylab="CD estimated proportion", 
       xlab="hp")
  for(i in 1:100){
    lines(seq(50,335,length.out=512), mat[,i], 
          col=rgb(.25, .25, .25, alpha=.25), lwd=2)
  }


Oddly, there's a big dip in the maximum at 150hp, but we can get a crude approximation of a 95% CI for this (it'd be better to do more iterations to get better precision at the limits):
dy$x[180]  # [1] 149.8337
round(sort(mat[180,])[c(3,98)], 3)  # [1] 0.000 0.289

