Visualising developmental binary data: plot model fit or data means? In my disciple (developmental psychology) it is becoming common to visualise data as the fitted model with CIs, rather the actual data (means with CIs). The typical example is where we are interested in the development with age of the tendency to choose 1 rather than 0 in a simple binary choice task. Data is visualised as the fit of the logistic model (e.g. http://www.pnas.org/content/110/36/14586.figures-only).
With my own data I am comparing these methods. I am initially attracted to plotting the fitted model, because (a) I think it looks cool, (b) everyone is doing it, and (b) more seriously, because the method seems like it focusses on signal rather than noise. The problem is, with my data set, it just doesn’t seem to do the data justice. Have a look. One of these graphs is model fits (separate for two conditions), the other is mean at each age (each +-95% CI). The most salient problem is that the model implies that the two conditions are quite likely to have different starting points at 4 years, which is clearly not the case from the data. The models are the most simple logistic regressions possible (R glm default settings for family=”binomial”).  I have also tried included a quadratic age term. The fit then looks a bit more like the actual data, but the conditions still don’t converge at 4 years, and above all the quadratic term is not significant and the model has worse AIC, so the quadratic term does not feel justified.
Any advice?

 A: Regarding the nominal question of showing the data or the model, sometimes it works out to show both. That is, put the main message, the model, in the foreground, but also include data in the background to give some sense of how well the model represents the data.
Here is an example showing a 3-parameter logistic model on @Jonas Berge's simulated data. I have no idea if this model is appropriate--it's only to illustrate the point of layering.

A: As you recognize, the real issue isn't whether to plot the data or the model, as the title puts it: it's that the model doesn't do justice to the data.
The data suggest an interaction between the treatment group and the shape of the relation of age to the outcome variable, which your original model didn't allow for. If the GAM doesn't work, you might consider treating age as an ordinal predictor variable and include an interaction with treatment group in a logistic regression. With 237 participants and about half making each outcome choice, you might have enough data for this to work.
A: I thought I'd answer my own question because I ended up with a solution I am happy with and is a bit simpler than the other suggestions. I have this rather nice looking figure:

The figure is simply the predictions +/- 95% CI (and raw means as dots) from this model with two spline knots (condition, age with knots, and the interaction):
knotLoc <- quantile(d$Age,probs=c(1/3,2/3))
mod <- glm( Response ~ Condition + ns(Age,knots=knotLoc) +
    Condition:ns(Age,knots=knotLoc),
    data=d, family="binomial"
)

As you can see the model works really nicely to generate the confidence intervals showing at what ages the conditions differ. However, I also note that in this model, neither the main effect of condition, nor any of its interactions with the age splines, are significant. This puzzles me a bit, because a really simple linear logistic regression with no splines show a very clear significant effect of condition. If anyone has any comment on that, I'd be interested. But for my purposes I think it may do to show the significance with the non-spline model, and the prettiness with spline model.
