I have two questions regarding modeling a 3-level multilevel model in R.

I have a dataset of different variables that were assessed 4x as part of a longitudinal study. At each of the four assessments that were 1 month apart, participants filled out questionnaires (= one mean score per participant per assessment) and used daily diary/ecological momentary assessment (EMA) over the course of 4 days. My issue is that the EMA data consists of individual entries per day per assessment (e.g., 3 entries on 4 days = 12 total entries per assessment). This data is on a different level than the questionnaire data (clustered within a person in addition to the timepoint level).

I am unsure how to calculate a model where I can include the EMA data as a predictor of the questionnaire data while acknowledging the individual entries.

So far, I have used a two-level model where I aggregated the EMA data to yield a mean per timepoint so it matches the data from the questionnaire. So level 1 is the within-person variation across the three timepoints, and level 2 is the between-person variation.

However, because I aggregated the EMA data, I lost the advantages of the dense data that was collected. I don't know how to perform a 3-level MLM where I include the individual EMA entries per day nested in each assessment.

As it is now, my data is in long format (and in addition to the original long format, I have also centered it around the grand mean and within clusters (=individual person), but that is not included here in this sample data). I hope I described the data type sufficiently.

structure(list(id = structure(c(1, 1, 1, 
                                1, 2, 2, 2, 2, 
                                3, 3), label = "serial", 
format.spss = "N12", display_width = 12L), 
               time = structure(c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2), 
format.spss = "F4.0"), 
               gender = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2), 
label = "gender", format.spss = "F3.0", display_width = 5L, 
labels = c(`na` = -9,                                                                                                                                               female = 1, male = 2, other = 3), class = c("haven_labelled",                                                                                                                                                                                                  "vctrs_vctr", "double")),
               questionnaire_mean = structure(c(9, 7, NA, NA, 6, 6, 13, 8, 
               11, 7), format.spss = "F8.2", display_width = 10L), 
               ema_mean = structure(c(4.21, 5, NA, 5, 2.6, 
                                      NA, NA, 5.9, 2.95, 2.79), 
format.spss = "F8.2", display_width = 10L)), row.names = c(NA, -10L), 
class = c("tbl_df", "tbl", "data.frame"))
  • 1
    $\begingroup$ Welcome to CV. Questions about coding are off topic here, but your question also contains statistical content so I am voting to keep it open. But, if you edit it to emphasize the statistics, then it will be less likely to be closed (unless, of course, your question really is about code). $\endgroup$
    – Peter Flom
    Commented Nov 30, 2023 at 18:21

1 Answer 1


In a multilevel model, the lowest level of the data is the limiting factor in terms of frequency or repeatedness of measurement. This is the level at which your outcome is measured. If the outcome is measured at four time points, as in the questionnaire_mean variable in your example data, then any predictors could at most, be repeatedly measured four times. Hence your approach of creating a mean of the EMA data is correct. You could also consider modeling the dispersion in each individual's EMA ratings at a given time point to capture some of the within-person variability you are interested in. In your EMA data, you could do something like the following:

ema_data <- ema_data %>% 
          group_by(id, time) %>% 
          mutate(ema_sd = sd(ema_response, na.rm=TRUE) %>% 

Note that you could use any number of measures of spread or dispersion - I just went with the standard deviation b/c it is well known. Then you could merge ema_sd into your outcome data frame and enter it as another time-varying predictor in your model for the questionnaire_mean:

m1 <- lmer(questionnaire_mean ~ 1 + gender + ema_mean + ema_sd + (1|id), data = dat)

Other related approaches that might give you more reliable estimates of the ema_means would be to move to a multilevel structural equation modeling framework (lavaan in R) with a measurement model for the EMA data points at the lowest level. This gets complicated and again, will only give you better estimates of the time point specific means.


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.