# crossed and nested random effects differences in a specific case in lme4 in R

I have the following data:

data <- data.frame(id_pers=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
Income=c(10000, 12000, 1000, 9000, 10000,
16000, 100000, 120000, 119000, 10000, 11000),
family=c(1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5),
BIRTH_ORDER=c("firstborn",  "notfirstborn",
"firstborn", "notfirstborn", "firstborn",
"notfirstborn", "firstborn", "notfirstborn",
"notfirstborn", "firstborn", "notfirstborn"),
income_father=c(12000, 12000, 20000, 20000,
15000, 15000, 92000, 92000, 92000, 9000, 9000),
income_father_cat=c("low", "low", "middle",
"middle", "middle", "middle", "high", "high",
"high", "low", "low"),
region=c("south", "south", "north","north",
"west", "west", "south",
"south","south","north","north"))


I want to find the effect of being part of a certain family. What does that mean? How much of the variation in Income of a certain person can be explained with being part of that family? In lme4, we would write: baseline_model <- lmer(Income ~ (1 | family), data= data)

At the end I want to test, whether the Interclass-correlation-coefficient (ICC) is getting smaller after controlling for a certain family-internal or family-external factors. For most of the variables it's pretty clear what to do. For example in testing the family-external effect of the region or the family-internal (categorized) income of the father

region_model <- Income ~ (1 | family) + (1 | region) , data= data)

father_model <- Income ~ (1 | family) + (1 | income_father_cat) ,
data= data)


All of these effects I will treat as crossed-random-effects, because the "sub-classes" belonging to every family. Nice explanation I read here: Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?

As far as good!

But now it's the question how to handle to family-internal factor of being first-born or not firstborn sibling within a family. There will always be at least one first-born individual but eventually more than one not first-born siblings (as in family "4"). How do I handle that problem with using either treating BIRTH_ORDER as a crossed random effect or as a nested random effect?

crossed_model <- Income ~ (1 | family) + (1 | BIRTH_ORDER) ,
data= data)

or

nested_model <- Income ~ (1 | family) +  (1 | family:BIRTH_ORDER) ,
data= data)


What does these different models tell us about the family background?

What do you think about it?

It doesn't really make sense (IMO) to treat birth order on its own as a random effect, as the levels are not exchangeable: that is, it wouldn't really make sense given birth order "first"/"not first" to switch the labels (whereas if you have families "Jones", "Smith", "Woodward", relabeling the families won't change anything in the interpretation). Therefore, I would go with (1|family:BIRTH_ORDER), which says that birth order within family matters. In fact, you might want (BIRTH_ORDER|family) (or (BIRTH_ORDER-1|family)), which would allow for a correlation in the effects of different birth orders within families.
(Furthermore, since BIRTH_ORDER only has two levels, it wouldn't really be practical to treat it as a random effect - unless you have strong Bayesian priors it doesn't really work to fit factors with fewer than 5-6 levels as random effects.)