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:
model <- gam (y ~ s(age, by = as.numeric(group=="green")) + s(age, by = as.numeric(group=="red")) + factor(group), family = binomial)
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.
I made a simulation based on your graph and I assumed there were 30 subjects per age/group pair, n = 360 in total.
n <- 30
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)
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:
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:
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:
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.
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!