longitudinal cross membership multilevel model

Question

I am compiling a data of club formation in which a person (id, the basic unit of observation) can belong to more than 1 club, while a majority of the sampled population do not belong to any clubs. Also, in a given year, more than 1 club were formed, so a person may share overlapping membership in multiple clubs with other persons; hence, the data has a cross membership structure (but not exactly a nested structure unless all clubs formed in a given year do not share overlapping membership). Finally, the data I have have observations spanning across 25 years.

To conceptualize the above description, my data structure looks somewhat like this:

there are many id(s) like id_5 who does not belong to any clubs in a given year.

The outcome of interest is to study a person's degree of involvement in social activities (on an ordinal scale). One way to code the data is to have all individuals (id(s)) nested within year(s) and then have their club affiliation crossed-classified. This should look like the format suggested in this statalist thread; however, the data I have also contains many time-varying individual attributes which I would like to use as regressors so the recommended Stata xtmixed syntax may not suit my case.

Alternatively, I am also thinking about modeling the outcome variable for each club (formed in a given year) separately as one panel so as to make those time-varying individual attributes "constant" in that year. Suppose the data records j clubs formed among n persons between the 25-year time span in the data, then I should have j panels and a total of j * n rows. But I am not sure if this makes any statistical sense.

Have anyone here dealt with data of this property before and could share insights on how to code and estimate this kind of data using Stata and/or R's lmer function?

####### updates

Much thanks to Erik for commenting on this. My earlier communication with Stata statistician suggests that one should use "wide format" to code each club as separate columns and then use each "id-year" observation's club membership to fill in each cell entry below those separate club columns as shown below. (also recommended using xtmixed command)

as compared to the hierarchical format included in Erik's reply.

The wide format coding essentially treats each club affiliation as a separate level (to borrow a jargon from hierarchical model), it needs fewer rows (for each id, should it belongs to more than 1 club in a given year) but can get unwieldy as the number of higher level units (e.g., clubs) increases. Also, it's unclear how the outcome variable, y, should be coded. Some recent papers (e.g., Mo and Wellman (2014) follow Simmel (1955) and Breiger (1974) by treating such cross membership as meta level and use the share of one's involvement in a group as cell entry (for example, if a person joins three clubs in a given year, then its cell entry in the three club columns are all 0.333), but it's also unclear if researchers should take other club members' (say, j's and k's) club involvement into account when calculating person i's club involvement as well as what these cell entries add up to within an id or in the cross-section.

Erik's method makes more intuitive sense but I suspect more entries associated with a particular set of ids (because they are more connected, having more friends, etc.) in a given year will create unbalanced panel. My other proposed solution will be to treat each club formation from the full sample (all possible dyads) or apply cluster errors by club id, but I am not sure if this will violate any distributional assumptions.

Answer 1

This is indeed a complicated analysis and the resulting dataset, I'm guessing, is also a beast. My advice would be to start simple. Figure out the model and syntax for a single year. Once you have that and feel good about it, you can ramp up the complexity by adding predictors and then multiple years' worth of data for each individual/club.

It would be helpful to have a snapshot of your data (maybe just for two players) to see how you have set it up. But assuming that it is something like so...

ID Year Club Outcome  
1   1    1     3.5
1   1    3     3.8
2   1    1     2.1
2   1    2     4.2
3   1    2     1.9
3   1    4     5.7
...

The Stata code to get a variance components cross-classified model without predictors from this hypothetical dataset is mixed outcome || _all: R.Club || ID: or in lmer lmer(Outcome ~ 1 + (1 | Club) + (1 | ID)).

Supposing you get this running, the addition of year into the model would be straightforward. In Stata mixed outcome || _all: R.Club || ID: || Year: and in lmer lmer(Outcome ~ 1 + (1 | Club) + (1 + Year/ID)).* You may also want to play with estimating the year random effect separately in lmer.

To the degree that you have time-varying predictors, which will be constant within-person and within-year, then you would simply add them as "fixed effects" predictors to the above models. You may also want to consider adding the person means for each of those predictors to disentangle the within- from between-person associations between these predictors and your outcome.

Regarding your idea to aggregate all the data up to the club level for a given year, you can certainly do that to create a panel dataset. But you can easily incorporate club-level predictors into these mixed models, and to the extent that those predictors do not vary by year, they will explain club outcomes in aggregate. If they vary by year, then you would again want to compute club means for the variables across years and include those as predictors in addition to the year-varying versions of the predictors to disentangle the within- and between-club associations of those variables with your outcome.

*Edit: I am not as familiar with lmer as I am with mixed, you may want to estimate the club effect as also nested within year - lmer(Outcome ~ 1 + (1 + Year/Club) + (1 + Year/ID)). I'll let more experienced lmer users weigh in on that.