Exploring shifts in response to dichotomous dependent variable

Question

I have one dichotomous dependent variable (buried with grave goods or not) and a series of categorical and continuous independent variables (age at death, year of death, sex, socioeconomic status) for a dataset of c.5000 individuals. Therefore, I have used binary logistic regression in SPSS to explore how each of the independent variables impact the likelihood of the individual being buried with grave goods.

However, in my preliminary analysis, I found that the proportion of individuals buried with grave goods seemed to increase until around the age of 10-12 and then plateaued. I have two questions related to how I can explore this trend:

1. I performed a binary logistic regression for individuals aged 1-12 and another for individuals 13-30 and found age was a significant contributing factor in the first regression and not in the second. While this does demonstrate that the pattern is likely there, I am unhappy with this approach and have not been able to find good examples of how these patterns are normally tested for significance (while controlling for the other variables). Any suggestions or examples of how this has been done would be greatly appreciated.

2. I am interested in identifying the age at which this trend shifts so that I can then see if there are differences between the socioeconomic groups however I do not have a statistics background and am unsure how best to go about this.

Answer 1

It is usually best to use a model with all the data, and then include interactions. I discussed this at Separate Models vs Flags in the same model, see also How to test a curvilinear relationship in a logistic regression.

For your data I would try a generalized additive model (see tag gam), a very good and flexible implementations is mgcv R (there might be something in SPSS, but I do not know about that). In reality, gam's as implemented in mgcv can be seen as a modern, flexible implementation of glm's (generalized linear models), and what I show here could be done, more or less, with glm plus some spline packages. But the modern, flexible way is shown below. So, something like

library(mgcv)
yourdf <- data.frame(grave_goods=, age_at_death=, 
                     year_of_death=, sex=, 
                     ses=, ...)
mod0 <- mgcv::gam( grave_goods ~ sex + ses + 
                    s(age_at_death) + 
                    s(year_of_death), 
                    family=binomial, data=yourdf, 
                     ... )
summary(mod0)

Here I assume that grave_goods, sex, ses are factors, and that ses has few levels (probably 2?). The above is then a base model, with smooth effects, but without interactions. Interactions can be modeled using terms such as

s(age_at_death, by=sex)

or s(age_at_death, by=ses)

or even

s(age_at_death, by=interaction(sex,ses))

(same with year_of_death). For interactions between the two numerical variables we have to be careful, since they will have very different ranges ... but something like the following:

s(age_at_death) + s(year_of_death) + 
ti(age_at_death, year_of_death)  

Here is some useful material on mgcv made available by ith author Simon Wood.