Model specification for causal inference with longitudinal data and multilevel modeling

Question

We studied the effect of media of delivery (i.e. augmented reality vs video) of educational content on memory (“binary_memory_score”) over time (“time”). We had 20 augmented reality experiences, with three memory questions specified for each experience(each multiple-choice with four categorical not ordered options; here’s an example question/answers -- “What colour was the makeup of the character in the theatre show you saw” / ‘red’, ‘green’, ‘blue’, ‘yellow’.). Videos were screen recordings of the augmented reality experiences.

We have two conditions (augmented reality vs video) for which we recruited randomly and separately (between subjects). Note that video were actual video recordings of the augmented reality experiences (1:1 relationship). Participants were asked to engage with educational content and answer memory questions. For both groups, participants were recruited into four “cohorts”, representing different educational themes (e.g. cohort 1 = art, cohort 2 = science), meaning cohort is also a between-subjects variable (although we have no reason to test whether cohorts differ as they were only loosely linked to specific themes). Within each cohort, there were five specific experiences all cohort participants experienced on the topic of the cohort (in cohort1, we have experience1, experience2, experience3, experience4, experience5; in cohort2, we had experience6, experience7 etc).

We think experience_question should be nested within experience as each experience has 3 separate questions that only apply to it -- So experience A has questions 1, 2, 3 and experience B has questions 4, 5, 6. We think (participants) id should be nested within cohort as each cohort has unique participants. So cohort A has participants 1, 2, 3, and cohort B has participants 4, 5, 6.

After the initial memory test after doing each experience (time_0), approx 25% of participants were asked the questions again later (time_later) either after 1-month (gap1); another 25% of participants after 6-months (gap6). Using an example to flesh this out: if we recruited 1000 people at time zero for a given cohort + media, at t1 we invited a randomly selected 500 people to answer the Qs again, but stopped recruiting at approx 250. At t2 We invited the remaining 500 (not invited before) people to answer Qs, stopping again at 250. So we have a 50% retention because out of 1000 initial recruits, only 500 joined the follow-up study.

We thus have a within-subject factor of time (time_0 vs time_later) and a between-subjects factor of gap (gap1 vs gap6). We currently have memory as a binary variable (correct vs incorrect answer) as a “score” for individual memory questions (“exp_question”). Here is our proposed analysis with glmer:

binary_memory_score ~ 1 + media*time*gap   + (1 | experience / experience_question) + (1 + time | id) + (1 | cohort / id),  family = binomial("logit"), nAGQ=0, control=glmerControl(optimizer = "nloptwrap"), data=combined_data_cor_screened, contrasts = list(media = contr.sum, time=contr.sum, gap=contr.sum))

We hypothesise that people who watch the video (media) will perform more poorly at the memory task (binary_memory_score), compared to people who do the augmented reality experience. We also hypothesise that memory scores will delay more rapidly over time for people who watch the video, compared to people who do the augmented reality experience. We plan to test these hypotheses via the emmeans package by exploring the media * time * gap interaction.

Answer 1

Here is a summary of variables based on your descriptions. I updated my answer after you provided clarifications.

• Response: response. For the three memory questions created for an experience, each was a multiple-choice question with four options, out of which only one option was correct. There can be multiple ways to code this variable, as discussed below.
• Treatment: media, binary, 0 = video recording and 1 = augmented reality. We should expect positive coefficient if AR enhances memory. Thus, the hypothesis should be one-sided.
• Materials: 20 different experience categories for each media. Since videos were screen recordings of the augmented reality experiences, the 20 augmented reality experiences had a one-to-one relationship with their corresponding 20 video experiences.
• Tools: 60 unique categories of question, three for each experience. In each wave of study, a participant received 5 experience categories and answered 5 * 3 = 15 questions.
• Participants: Each participant id was assigned to either augmented reality vs video, so id and media are many-to-one relationship.
• Theme: four theme groups. Each participant id was assigned to only one theme, another many-to-one relationship. All participants in each theme were assigned five specific experience categories that were not experienced by participants in other theme groups.
• Repeats: binary time, 0 = initial test, 1 = follow-up test. However, only about 50% of participants joined the follow-up test. Each id corresponds to one or two of the time categories.
• Interval: binary gap, 0 = one month afterwards, 1 = six months afterwards. For the participants who remained in the follow-up test, each id corresponds to only one gap category. This time interval was randomly assigned by the follow-up invitation.

A few variables that I consider important but are not mentioned

• The duration elapsed between the experience and the question. Answering a question regarding the beginning of an experience may get different results than answering regarding the end.
• Sequence of the experience and question. Starting with Experience 5 may get different results than starting with Experience 2.
• Attributes of experience. For example, playtime, music speed, brightness, time of day.
• Attributes of participants. For example, age, cognitive capability, memory capability, height/weight or neck muscle strength since AR device can bother those who do not like heavy headset and affect attention.
• Attributes of options. Position in a drop down menu (e.g., top vs. bottom, left to right, or ABCD), length of prompt (number of characters), whether the sequence of questions and options is randomized. Maybe participants tend to select the first/last/shortest/longest option regardless of its content.

The grouping structure is partially clear as I outline below. What is unclear is whether the four options within each question require grouping, which depends on the presentation of options.

• theme / experience / question;
• theme / id / time;
• media / id / time.

First, I would test different ways of constructing the response variable. Depending on whether the presentation or sequence of each option carries any meaning, you may or may not need to create rows for each option.

1. The simplest case is to treat correctness of an answer to each question as a binary correct, forming 15 rows for each id in a wave of study. This is appropriate if selection of the other two options offers no additional information. This approach requires a mixed-effect binary logit model. See lme4::glmer() and glmmTMB::glmmTMB().
2. A more complex case is to treat answer to each question as four TRUE/FALSE records in the variable selected, forming 60 rows for each id in a wave of study. Even if four options in each question are mutually exclusive—selecting one option means that the other three options are not selected—there could be other dimensions to compare options, such as the relative position and length. In this case, the correct option should be marked by an additional TRUE/FALSE variable, correct, to distinguish from the selected option where selected == TRUE. This approach requires mixed-effect multinomial logit regression. See https://cran.r-project.org/web/packages/mlogit/vignettes/c5.mxl.html and https://cran.r-project.org/web/packages/mlogit/vignettes/e3mxlogit.html. The memory capability is represented by the predicted probability of the option where correct == TRUE. Another benefit of using mlogit is that we can constrain the sign of random effects. For example, time should have a negative sign even if its coefficient varies by id. This can be controlled by a lognormal distribution of the random slope of time.
3. Construct the response variable count as the number of correct answers to memory questions in each experience, forming 5 rows for each id in each wave of study while dropping the question distinction. This might be beneficial because binary responses contain the least information compared to count, ordinal, and continuous ones. For longitudinal binary response, mixed-effects modeling will discard any group entirely if the binary response stays constant within the group. For example, if a participant id == 7 answers all questions correctly, all the rows where id == 7 are removed from model fitting if id has random intercept. This results from the random intercept at id == 7 being +Inf to reach probability of one, breaking the estimating routine. Therefore, mixed-effects binary logit models may often run into numerical issues. For each id-experience-time combination, the count of correct answers is a limited integer ranging between 0 and 3. Count models (such as Poisson and negative binomial regression) that allow the response to take any natural number may not work well. The most appropriate technique is ordinal regression, particularly cumulative link models. I would use ordinal::clmm() or ordinal::clmm2() for cumulative link models with mixed effects. See https://cran.r-project.org/web/packages/ordinal/vignettes/clmm2_tutorial.pdf. Even if binary regression suffices, this package lends support to "scale effects."

Second, I recommend estimating three models separately, one for the initial test, another for whether a participant remained in the follow-up test, and the last for participants who joined both initial and follow-up test.

• The initial-test model uses data where time == 0 and each id answered 15 questions. Focusing on only the participants who joined both tests causes sample selection issues.
• The retention model uses data where each id has only one row and the binary response variable retention is TRUE if the id joined the follow-up study. It is plausible that participants who had greater memory capabilities remembered better about joining the follow-up study, or participants who performed better in the initial study would feel more comfortable to accept the invitation to the follow-up study, so a self-selection process was possible that attracted those who had higher memory. We should use the memory score at time == 0 as a predictor. We had better see that no predictor in this retention model is significant, to demonstrate a random attribution process and no sample-selection bias. See http://www.jstatsoft.org/v27/i07/.
• The follow-up model uses data where retention == TRUE and each id answered 30 questions over two waves of study. Due to the attrition in the follow-up study, setting up the before-and-after model specification will discard participants who did not remain in the follow-up test.

Once the response is transformed into a count of correct answers across three questions, each id has only one count measurement for each among the five experience categories at each time. Tentative specifications of the three models are

clmm(count ~ 1 + media * theme + (1 | theme : experience) + (1 | theme : id), link = "logit")
glm(retention ~ 1 + media * theme + count, family = binomial())
clmm(count ~ 1 + media * time * gap * theme + (1 | theme : experience) + (1 | theme : id), link = "logit")

Third, it is worth considering which variables should exhibit random intercepts or slopes.

• With only four categories of theme, there may not be enough data to estimate random intercepts by theme. It is implausible to consider these four themes as a random sample of all possible themes and the findings generalizable to other unobserved themes. On the other hand, it may be important to understand how one theme contrasts with another in memory scores and mediates effects of media, time, and gap. Therefore, I consider fixed effects only of theme. Nevertheless, random effects of other categorical variables can be grouped within theme.
• It is unclear to me whether the data support random slopes and correlation between random intercepts and random slopes, as the fixed-effect predictors are mainly binary. You can test random components (1 + time + gap | theme : id) and (1 | theme : id) + (0 + time | theme : id) + (0 + gap | theme : id). I suspect that adding random slopes may lead to estimates at boundaries, such as infinite coefficient, very large standard error, or zero variance. Watch out for these signals of numerical issues.