# 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. You seem to have fitted a linear model when the association between age and the dependent variable may not be linear, especially for the green group. You can try this:

library(mgcv)
model <- gam (y ~ s(age, by = as.numeric(group=="green")) + s(age, by = as.numeric(group=="red")) + factor(group), family = binomial)
plot(model)
summary(model)


This should give you information whether the effect of age on the dependent variable is different between groups, and it enables the effect of age to be non-linear, which might be the case with the green group.

EDIT:

I made a simulation based on your graph and I assumed there were 30 subjects per age/group pair, n = 360 in total.

library(mgcv)
n <- 30
set.seed(1)
age40 <- rep(4, n)
group40 <- rep(0, n)
outcome40 <- rbinom(n,1,0.2)
age41 <- rep(4, n)
group41 <- rep(1, n)
outcome41 <- rbinom(n,1,0.2)
age50 <- rep(5, n)
group50 <- rep(0, n)
outcome50 <- rbinom(n,1,0.23)
age51 <- rep(5, n)
group51 <- rep(1, n)
outcome51 <- rbinom(n,1,0.52)
age60 <- rep(6, n)
group60 <- rep(0, n)
outcome60 <- rbinom(n,1,0.28)
age61 <- rep(6, n)
group61 <- rep(1, n)
outcome61 <- rbinom(n,1,0.67)
age70 <- rep(7, n)
group70 <- rep(0, n)
outcome70 <- rbinom(n,1,0.48)
age71 <- rep(7, n)
group71 <- rep(1, n)
outcome71 <- rbinom(n,1,0.63)
age80 <- rep(8, n)
group80 <- rep(0, n)
outcome80 <- rbinom(n,1,0.61)
age81 <- rep(8, n)
group81 <- rep(1, n)
outcome81 <- rbinom(n,1,0.67)
age90 <- rep(9, n)
group90 <- rep(0, n)
outcome90 <- rbinom(n,1,0.67)
age91 <- rep(9, n)
group91 <- rep(1, n)
outcome91 <- rbinom(n,1,0.64)
age <- c(age40, age41, age50, age51, age60, age61, age70, age71, age80, age81, age90, age91)
group <- c(group40, group41, group50, group51, group60, group61, group70, group71, group80, group81, group90, group91)
outcome <- c(outcome40, outcome41, outcome50, outcome51, outcome60, outcome61, outcome70, outcome71, outcome80, outcome81, outcome90, outcome91)
group <- factor(group)
red <- c(sum(outcome40)/n, sum(outcome50)/n, sum(outcome60)/n, sum(outcome70)/n, sum(outcome80)/n, sum(outcome90)/n)
green <- c(sum(outcome41)/n, sum(outcome51)/n, sum(outcome61)/n,     sum(outcome71)/n,     sum(outcome81)/n, sum(outcome91)/n)
lines((4:9), green, col="green")
age <- age + rnorm(length(age), 0, 0.01)
M1 <- glm(outcome ~ age + group, family=binomial)
M2 <- glm(outcome ~ age + I(age*age) + group, family=binomial)
M3 <- glm(outcome ~ age * group + group * I(age * age), family=binomial)
M4 <- glm(outcome ~ age + group * I(age * age), family=binomial)
M5 <- gam(outcome ~ group + s(age), family=binomial, method="ML")
M6 <- gam(outcome ~ s(age, by=as.numeric(group==0)) + s(age,  by=as.numeric(group==1)), family=binomial, method="ML")
plot((4:9), red, type="l", col="red", ylim=c(0,1),
main="Sample means",xlab="Age (years)", ylab="Proportion choosing 1")
AIC(M1, M2, M3, M4, M5, M6)
Summary(M6) I don't know how to create those confidence interval markers..

You may notice that I added a small random component to the age variable. This was done because GAM wouldn't run otherwise, so it seems you were right about GAM being picky about discrete variables. But by adding a small random component, it should be okay (I'm talking about adding a random normal variable with mean 0 and standard deviation 0.01 so it should not influence the results).

I haven't seen a reference to this procedure or seen it used before, I just came up with it as a solution right now. The way I see it, we just add a minimal amount of measurement error to what is already there (because I assume that most of the individuals in the study weren't studied on their birthdays). So adding this small measurment error should be okay, because if we had access to their real age, we wouldn't have this problem.

If you have R, you can see that the graph looks similar to yours. The AIC command gives the following output:

         df      AIC
M1 3.000000 445.1460
M2 4.000000 442.0030
M3 6.000000 439.7423
M4 5.000000 443.7428
M5 3.938145 443.6351
M6 5.670207 439.2681


The models M3 and M6 seems best. Indeed, when comparing these models to the rest using anova(Mx, My, test="LR"), M3 and M6 are superior to the rest, but there is no significant different between the models:

anova(M3, M6, test="LR")

Analysis of Deviance Table

Model 1: outcome ~ age * group + group * I(age * age)
Model 2: outcome ~ s(age, by = as.numeric(group == 0)) + s(age, by = as.numeric(group == 1))
Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
1    354.00     427.74
2    354.84     427.93 -0.84363 -0.18533   0.5974


The GAM model (M6) is more easily understood visually, so I'd pick that one:

plot(M6) This is the fitted model for group 0 (red) And this is for group 1 (green).

So we see that the relationship for group 0 is linear and non-linear for group 1. Bear in mind that the variable on the y axis is not the logit but the unexponentiated estimate, so that why the curve for group 0 is linear and does not assume the classical s-curve.

You can see that both curves starts at about -1 at age 4, which together with the intercept of -0.47 gives an odds of exp(-1.47) = 0.23, which corresponds to a probability of 0.23/(1+0.23) = 0.19. At age 9, the curves are at about +1, so the odds are about exp(1-0.47) = 1.70, which gives a probability of 1.7/(1+1.7)= 0.63. As you can see these figures closely match the probabilities I specified (p 0.2 for 4-year-olds and 0.67/0.64 for 9-year-olds in groups 0 and 1, respectively). So using this approach, the curves will indeed match at ages 4 and 9, just like in your data.

I think you might be able to create better graphs if you could extract the smoothers with confidence intervals from the M6 object, and you could then add the intercept to each smoother and exponentiate and convert to probabilities, and plot both curves in the same graph. Then it would be really easily understood. Unfortunately, I don't know how to extract the smoothers so I can't do that yet.

Finally, we can see that the smoothers are significant:

summary(M6)

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4673     0.1651   -2.83  0.00466 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df Chi.sq  p-value
s(age):as.numeric(group == 0) 1.000   1.00  21.15 4.24e-06 ***
s(age):as.numeric(group == 1) 3.156   3.67  40.20 8.09e-08 ***


So it seems that using this non-linear approach might be a good idea and worth a try.

EDIT 2:

Actually.. now that I'm thinking about it, it is also possible to plot like this:

plot(age[which(group==0)], M3$fit[which(group==0)], type="l", col="red", ylim=c(0,1)) lines(age[which(group==1)], M3$fit[which(group==1)], col="green") This is similar to the curves above. I don't know how to add the 95% CI lines here. So in this case, the quadratic terms with interactions between age variables and groups give similar results and the likelihood ratio test shows no significant difference between this model and the GAM model.

Either way, provided the results in the model you choose for your data are significant as well, it seems that while there is a clear positive (linear) effect of age, it seems that the grouping (red/green) only makes a difference in a certain age range, and that the difference seems to disappear when the children grow older. I think you should show outputs from a GAM or quadratic model rather than your linear model, since it's obvious upon visual inspection that the linear model isn't appropriate. I'm looking forward to seeing your results!

• From the graph it seems like age has only 6 possible outcomes (age 4-9). Will gam, or other smoothers, be reasonable with such a discrete variable? – Maarten Buis Sep 3 '15 at 9:59
• I'm not sure.. it should be better than assuming a linear effect of age though. Perhaps age should be used as a factor instead, but the data may be too small for this since the interaction terms between age and group should be tried then (in the the first graph it looks like there is a different effect of group on the outcome depending on age). In any case, the models could be compared to see which model provides the best fit and lowest AIC/BIC, and if the GAM is better than the others it should be acceptable I think. – JonB Sep 3 '15 at 10:11
• Thanks for the suggestion, although I also worry because age is indeed somewhat discrete, I only have years-old (not exact age). I like the idea of automatically fitting the right sort of polynomial, as I read here that gam does, but the idea of splitting the data into different age sections seems excessive. I will have a play and see how it looks. – Amorphia Sep 3 '15 at 10:12
• The thing about a linear effect of age: because it's a logistic regression, I was hoping it would fit anyway. A linear but logistic model is not a straight line plotted as here: I can see in my mind's eye both the red and green curves as being part of the typical s-shaped logistic curve, and was a bit disappointed that the models did not really end up fitting like that - I'm not sure why not. – Amorphia Sep 3 '15 at 10:19
• I don't think that the green curve looks like a typical curve in a logistic regression model, but the red one does. It looks like there is some kind of interaction between age and group, so try taking this into account in your quadratic model and see if it helps: y ~ age * group + age^2 * group, or perhaps just y ~ age + age^2 * group. – JonB Sep 3 '15 at 10:22

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. 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.

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.