I have a question about specifying a mixed model with crossed data across multiple levels. Consider the following situation:

Let’s say I collected data on Americans’ perceptions about the importance of 85 discrete and well-validated managerial skills at four different points in time (Wave 1, Wave 2, Wave 3, Wave 4). In Wave 1, a batch of 200 respondents rated each of the 85 skills on a scale of 0 to 100 based on how important they thought that skill is for successful managers to possess. The same protocol was used in each of the following waves, but new respondents were recruited for each wave.

An exogenous shock happened between Waves 1 and 2. The shock randomly assigned 10 of these skills to some treatment. I am attempting to measure whether the shock had a long-term effect on the mean ratings of the treated managerial skills. The unit of analysis related to my question is the managerial skill – which is repeatedly measured across four waves, but with a new cross-sectional sample of 200 respondents for each wave.

Here are the first few observations of the analysis file (called anfile) to show the structure of the data and the names of the variables. This version of the analysis file anfile is at the respondent-skill level (implicitly the respondent-skill-wave level because unique respondent ids will never appear in more than one wave).

enter image description here

In my first attempt at analyzing the data, I computed mean ratings for each managerial skill in each wave, discarded the variance produced from individual respondents (I call this data frame anfile_means).

Then I estimate a simple linear mixed effects model in R using the lmerTest package (note: wave and treatment are formatted as factors):

m1 <- lmer(rating ~ wave + treatment + wave*treatment + (1|skill), data = anfile_means)

The results are consistent with the hypothesis: the effect of Wave2 is large and positive, but there is a strong interaction between Wave2:treatment, Wave3:treatment, and Wave4:treatment. In other words, the shock induced a difference in mean importance rank for the treated cases that persists well after the shock’s occurrence.

However, I am concerned that my results are artificially significant because I discard a large amount of noise in the estimate and only use the mean rating for each managerial skill. I suspect that the p-values associated with the interactions (currently .0000000001) would be far less significant if the estimation strategy incorporate uncertainty into each mean managerial skill.

I am unsure how to decompose my data structure into a mixed modeling strategy. Here is my attempt to do so: 

  • I have three levels of data: Wave-Level, Skill-Level, and Respondent-Level.
  • Wave and Skill are crossed (i.e., contain “multiple memberships”)
  • Skill and Respondent are also crossed (i.e., each skill was rated by all 200 respondents in a wave)
  • However, it seems like Waves and Respondents are nested (?) because respondents only appear in one wave.

Below is a schematic of the nesting/crossing in the design. Note that i'’s are all unique, I just didn’t have room to provide them with unique identifiers in the illustration.

enter image description here

I am not sure how to specify the model with the appropriate syntax given this crossed design. I am also unsure if my last bullet point about respondents nested within waves is correct – or if I should be thinking of respondents as the second level, and the skills as the third, which would allow me to draw the diagram as nested as Wave-Respondent, but crossed at Respondent-Skill. Here is my attempt, using the anfile data (i.e., respondent-skill level):

m2 <- lmer(rating ~ wave + treatment + wave*treatment + (1|skill/wave), data = anfile)

Which yields results equivalent to my hypothesis, but which I worry I have misspecified.

How can I appropriately decompose the within-skill-by-wave variance in a mixed model syntax? And are my data crossed as I visualized them, or nested and crossed based on my alternate description in the paragraph above?  


1 Answer 1


It is correct to think of respondents nested within waves, but it is useful only if you model wave effects as random effects just like classrooms nested within schools. However, you have only four waves which are not a random sample of infinite waves of surveys but rather a limited experiments with treatments. So, it is appropriate to estimate wave effects as fixed effects as you did in both m1 and m2, and it’s great that you include interaction terms. Depending on the variable types, you may also want to include nonlinear functions of treatment if you measured the treatment dosage more precisely or simplify wave effects with some functional form instead of categorical levels.

Throwing away waves out of the question, you have respondents and skills crossed in clustered errors: Each skill was measured 200 * 4 = 800 times by 800 unique respondents, while each respondent was measured 85 times across 85 different skill items, making 800 * 85 = 68,000 observations in total. So, your diagram of hierarchical structure is misleading. Crossed and nested structures are only relevant for factor variables that group clustered or correlated error terms through either random or fixed estimators, usually counting into the intercept and slopes.

Your data structure and variable coding look good, as shown in the spreadsheet. Your model m1 is okay. It can be used as part of sensitivity analysis to inform audience whether the substantive conclusions change after you aggregate data. It is unclear when you claim “The shock randomly assigned 10 of these skills to some treatment.” Were (1) 10 out of 85 skill items set aside in the treatment groups for all respondents and all waves while the other 75 were always in the control group, or did (2) each wave or each respondent receive a different set of skills to be treated? This makes a difference in terms of whether it is possible to model fixed or random slopes of the treatment indicator by skill.

If (1) is true, it is best to have skill indicators as a set of 84 dummy variables and assess the treatment effect through the ten dummy variables representing the treated skills and their interactions with wave indicators. This allows you to see whether the treatment effect varies by skill, so called heterogeneous treatment effects. Do you consider 85 skills approximately a complete set of all skills required in managerial positions or a sample of a larger set of skills? If the former, estimate treatment variation by skill through fixed effects, where you only have random effects from respondents that contribute to the intercept:

m3 <- lmer(rating ~ skill * wave + (1 | respondent), data = anfile)

If the latter, you may have random effects of skills as intercepts, but you won’t be able to estimate skill-specific treatment effects, as there are at most only 10 different values of treatment coefficients through its random slope by skill. Juxtapose one random intercept by skill and another by respondent, since they are crossed.

m4 <- lmer(rating ~ wave * treatment + (1 | skill) + (1 | responsent), data = anfile)

You can also allow the treatment coefficient to vary by respondent for individual heterogeneous treatment

m5 <- lmer(rating ~ wave * treatment + (1 | skill) + (1 + treatment | responsent), data = anfile)

Or if you want both respondent-specific and skill-specific treatment effects, recode treatment into a factor with 11 levels with the base level representing the control group.

m6 <- lmer(rating ~ wave * treatment10 + (1 | skill) + (1 + treatment10 | responsent), data = anfile)

If (2) is true, opportunities may exist to assess heterogeneous treatment effects by skill via random slopes. So you have random deviations from the mean treatment effect both by skill and by respondent as follows,

m7 <- lmer(rating ~ wave * treatment + (1 + treatment | skill) + (1 + treatment | responsent), data = anfile)

Your model m2 looks incorrect because it’s missing random intercepts by respondent. You can test and include all other models for sensitivity analysis, and see if the increasing complexity is necessary. You may also want to check whether the error variance, not just the mean, is changed by treatment. To do that, you need to use the nlme package and function lme(). Refer to modelling random structure in lmer and nlme:lme for how to specify crossed random effects. After you create a new factor variable that takes a constant level across all observations perhaps through the dplyr package,

m8 <- lme(
  rating ~ wave * treatment10,
  random = list(
    one = pdIdent(~ 0 + skill), 
    one = pdLogChol(~ 0 + respondent + treatment10)), 
  weights = varIdent(form = ~ 1 | wave : treatment10), 
  data = anfile %>% mutate(one = factor(1)))

Your dependent variable as a rating ranging between 0 and 100 might suffer from limited-response problems that make the residuals deviate from a normal distribution although a moderate deviation is okay. You want to check the residual distribution carefully, especially when you need to test hypothesis and make inference about the treatment coefficient. To combat the error distribution assumption, you can consider using ordinal regression. The R package ordinal has functions “clmm()” and “clmm2()” that allow random intercepts and slopes. You could also consider Tobit regression or survival analysis for censored response, if ratings occurred at boundaries but actually meant more important than 100 or less important than 0. You could also try beta regression since both lower and upper limits of the dependent variable are known. There is also zero and one inflated beta regression if 0 and 100 ratings appear often, otherwise you may need to manipulate boundary values to 0.5% and 99.5% or 1/n% and (100 - 1/n)% for example. Beta regression with random effects can be done in package glmmTMB. If your 85 skill items are a reflection of a few intrinsic capabilities, then you may use factor analysis and structural equation modeling. If your treatment was administered on individual skill items but not on holistic capability, such dimension reduction techniques may not address your research question at all.

A good reference of mixed modeling techniques is Galwey (2014) Introduction to Mixed Modelling - Beyond Regression and Analysis of Variance.


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.