I know this kind of question has been asked before, but I can't find anything that clearly elucidates the issue. What is the 'right' way to model the follow situation:

Let's say I pair two people up randomly (on the basis of their sex) such that there are four possible groups: MM, FF, MF, and FM. I have them play some behavioral economics game, and I am interested in their before-and-after differences on some variable, let's say their happiness.

In R, I can think of several potential ways to model this relationship. Assume treatment is a categorical variable containing the above four categories:

lm(happy_score ~ treatment + education + income)

However, this doesn't account for the fact that, within every 2-person team, there will be correlations unaccounted for in the above specification.

So, we can add a fixed effect for team (a unique identifier for each 2-person team)

lm(happy_score ~ treatment + education + income + team)

This seems to not be getting at what I want, as it will compute a separate slope for each team, when I really want to just 'control' for the correlation that is likely to be shared between two individuals on the same team.

I've also been told to address this issue we can cluster standard errors at the team level, so:

  happy_score ~ treatment + education + income,
  data = data,
  clusters = team,
  se = "stata"

But I'm not sure what this is doing that is different from adding a fixed effect. Should I have both fixed effects and clustered standard errors?

I've also been pushed in the direction of mixed models, but I don't understand what they are substantively doing that is different from some combination of fixed effects and interactions:

 lmer(happy_score ~ treatment +  education + income 
                (1|team), data = data)    

Are there intuitive ways of thinking about each of these options?


2 Answers 2


The approaches you've outlined differ in terms of how you're using the team-specific data to measure the baseline happy score (happy score in absence of intervention). These approaches correspond to a "pooled" approach, an "unpooled" approach, and a "partially pooled" approach.

The partially pooled approach is the most reasonable of the approaches since it balances the team-specific information with information from the grand mean across teams (see http://www.stat.columbia.edu/~gelman/research/published/multi2.pdf). I explain why this is, below.

Pooled model: lm(happy_score ~ treatment + education + income)

  • This approach ignores variability in the baseline. The variation in the baseline happy score would be reflected in the standard error of the estimated treatment effect.

Unpooled model: lm(happy_score ~ treatment + education + income + team)

  • This approach will estimate a separate underlying baseline for each team. This approach is thought of as the unpooled approach since no information on baseline happy score is shared across the teams. The problem is that there is probably a fair bit of error in measurement of the baselines, since you only have 4 observations per team. Some teams with extremely high or extremely low initial scores will probably experience regression to the grand mean, irrespective of the intervention. So in this unpooled approach, the treatment effect is measured on top of this noisily measured baseline. The variation due to the poorly measured baseline happy score would be reflected in the estimated treatment effect, which would have a relatively high standard error.

Partially pooled model: lmer(happy_score ~ treatment + education + income + (1|team), data = data)

  • The multilevel model approach will provide estimates of each team's baseline that is "partially pooled" across teams. This means that team-specific means are "shrunken" towards the grand mean. The optimal degree of shrinkage/pooling is determined by the data. If there is lots of real variation (that lasts across time points) between teams, then the degree of pooling will be small. If most of the variation between teams is noise (i.e. it's due to measurement error or something that is not consistent across time points), then the degree of pooling will be large. With this partial pooling approach, your estimated treatment effect will benefit from the more accurate measurement of the baseline and should have a relatively low standard error.
  • The intuition for the partial pooling approach is nicely explained in Efron's classic paper - "Stein's Paradox in Statistics" (http://statweb.stanford.edu/~ckirby/brad/other/Article1977.pdf)

The robust approach should yield a similar answer to the multilevel model approach in this scenario, except you won't get the "shrunken" team-specific estimates that the multilevel model would provide. If this were a Poisson or logistic multilevel model, then estimated effects would differ, since the robust approach would yield "population average" effects while the multilevel model would yield "subject specific" effects.

Also for all the above approaches, I'd suggest to change the coding of the treatment variable to a 5 category variable [before (ref), after-FF, after-FM, after-MF, after-MM]. This coding is the constrained baseline approach to the analysis of cluster RCTs (see Hooper R et al. Analysis of cluster randomised trials with an assessment of outcome at baseline. BMJ 2018).

  • $\begingroup$ Thank you for the comprehensive answer, Kevin! A couple follow-ups: given that each team only has two people in it, does it make sense to add team-level fixed effects? From what I understand, adding fixed effects restricts us to only examining variability WITHIN teams, which is not what we want to do (we want to compare across teams with different treatment conditions, while controlling for the fact that members of the same team will have some shared characteristics --e.g. the same joint score on the game). Second, could you explain is the multilevel specification is actually doing that? $\endgroup$ Oct 26, 2018 at 1:38
  • $\begingroup$ I'm not 100% certain what is meant by "this means that team-specific means are "shrunken" towards the grand mean." $\endgroup$ Oct 26, 2018 at 1:39
  • 1
    $\begingroup$ In my understanding the treatment variable was: before=0 vs after=1 (not the MF, FM, FF, MM groups) / the outcome was the raw happy score (not the differences). The approach that I've outlined works on the understanding that you would have 4 observations per team (before/after x 2). You'd lose important info if you pre-calc'd the before/after differences. To measure group-specific effects you could create a 5 category treatment var (before, after-FF, after-FM, after-MF, after-MM). This assumes that all groups are the same in the "before" period, which is reasonable if they're randomly assigned $\endgroup$ Oct 26, 2018 at 2:55
  • $\begingroup$ Thank you for your help! Herein lies my problem: the 'treatment' in my experiment is whether you are randomly assigned a female or a male partner. Thus, there are four 'conditions' (MM, FF, FM, and MF). All four groups play the same game, and their happiness is measured before and after. I'm trying to quantify if individuals in one of these groups have different before and after changes than individuals in other groups. The problem is that we can't assume that 'all groups are the same in the 'before' period, because they'll crucially vary by sex. What's the right way to deal with this? $\endgroup$ Oct 26, 2018 at 3:01
  • 1
    $\begingroup$ The shrinkage estimate is shrunken towards the grand mean of all the teams by a certain shrinkage factor. So the shrinkage estimator lies somewhere between the team mean and the grand mean. z = y_bar + c(y - y_bar) where y_bar is the grand mean, y is the team mean, and 0 <= c <= 1 is the shrinkage factor. c = 0 is complete shrinkage (i.e. the pooled model), while c=1 is no shrinkage (i.e the unpooled model). Usually the c will lie somewhere in between. $\endgroup$ Oct 26, 2018 at 3:04

The first approach (i.e., lm()) postulates that there are differences in the average happy_score per team. However, it still assumes that measurements within a team are uncorrelated. If teams are truly random (i.e., you have a sample of teams), theoretically this approach violates one of the assumptions of standard asymptotics that requires that the dimension of the parameter spaces does not change with the sample size.

I am not totally familiar with lm_robust() but what I think it does is that it fits the model assuming independent error terms, and corrects the standard errors using a sandwich-type estimator. Two points regarding this approach: (1) if you have teams that miss the responses of some members (i.e., you have missing data), the derived results will only be valid under the missing completely at random assumption. (2) The sandwich estimator is known to protect against misspecification (i.e., here you have correlated data but you fit a misspecified model that ignores the correlations) but at the expense of power.

The last approach explicitly accounts for the correlations by including the random effect per team. Under the assumption that the postulated correlation structure is correct, this will provide you with correct inferences, and also be valid under the missing at random assumption (that is less stringent than the missing completely at random).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.