Estimating random slopes with missing data and comparing within/between-subject designs

[Jon Fawcett]

I am analyzing a data set for someone in which participants make two binary decisions (related to art preference) in either one or two conditions each. Individual participants were tested across 6 experiments split between these designs. This results in some participants with 1 response in Condition A and 1 response in Condition B (within-subjects design), some participants with 2 responses in Condition A and some participants with 2 responses in Condition B (between-subjects design). In addition to assuming that participants will vary in their overall propensity to make one decision or the other (random intercept) it also seems reasonable to consider the possibility that they would vary in their response to the experimental manipulation (random slope). A similar argument might be made for variation in response at the level of the experiments as well. Therefore, the appropriate formula seems to be:

stan_glmer(decision ~ condition + (condition | participant_id) + (condition | experiment_id), data = data_frame, family = binomial, adapt_delta=.999)

This fits the model just fine, and the results look reasonable. However, it is interesting that this model *does* work so well. Trying variations of this model using glmer results in an error concerning identifiability, presumably because the model is trying to estimate random slopes for some participants that have no responses in one of the two conditions. I believe that rstanarm deals with those cases by imputing slopes based on the joint posterior: This means (roughly) that the random slopes for participants in the between-subject experiments will be based on the random slopes for participants in the within-subject experiments. This seems fine, and the specific random slopes for the between-subject experiments are not important. My question is whether this understanding is correct or whether I am missing some pernicious element of modelling the data in this manner. Further, we would ultimately like to compare performance across experimental design, for example using:

stan_glmer(decision ~ condition * design + (condition | participant_id) + (condition | experiment_id), data = data_frame, family = binomial, adapt_delta=.999)

Here, design refers to whether a given experiment employed a within- or between-subjects design. Again, the model fits and everything appears to be in order, but I am a little cautious in part due to earlier critiques of within/between-subject comparisons of this nature using more “traditional” techniques (for example, Erlebacher, 1977). I believe this sort of comparison to be a non-issue with stan_glmer, but I wanted a second opinion before committing.

Thanks in advance!

[Jonah Gabry]

Hi Jon, 

I wouldn't say that rstanarm is "imputing slopes based on the joint posterior", although I think I know what you mean. It's not really imputation but more borrowing information. But yes I think you have the basic idea correct. 

Regarding non-identifiability, the priors (even the defaults) that stan_glmer uses should in general help make it more robust to non-identifiability than regular glmer. There can still be problems, but priors help a lot. Also, how many levels (unique values) does participant_id have? 

I'm not familiar with Erlebacher (1977), so I can't really say if the concerns in that paper are or are not an issue here. Can you give a brief summary?

Jonah

[Jon Fawcett]

Thanks, Jonah!

You are correct, this is precisely what I meant (albeit expressed in an imprecise manner, it would seem). As for the number of levels, I have somewhere in the order of ~600 participants and 6 experiments. But again, for approximately half of those participants there is 1 observation for each condition and for the remaining half are 2 observations in one of the the two conditions (meaning that those participants lack entirely data in the remaining condition).

As for Erlebacher (1977), the precise concerns voiced in that article are not relevant because he is talking about trying to compare designs using Frequentist ANOVAs. However, the thesis is that psychologists of that era had been drawing comparisons across within- and between-subject studies for some time in a misleading or wasteful manner. Erlebacher proposed to correct this by providing a variation of a standard ANOVA capable of making across design comparisons. I really just want to ensure that I am not analogously blind to an error in the current model. In short, if by running the aforementioned model:

stan_glmer(decision ~ condition * design + (condition | participant_id) + (condition | experiment_id), data = data_frame, family = binomial, adapt_delta=.999)

I am able to conclude that the effect of condition is (for example) larger for studies in which this variable is manipulated within-subjects relative to between-subjects (the condition x design interaction), I am all set!

[Jonah Gabry]

Hi Jon, can you share a bit of (potentially fake) data? For example, a small R data frame with the same structure as you have but with completely made up values would be fine (I don't need to see the real data if it is private). I just want to be 100% sure that I understand the distinction you're making between "design", "condition", and "experiment" and seeing an example would make that clear. I'm pretty sure I'm following you and that what you describe seems reasonable, but I don't want to inadvertently give you bad advice. 

Thanks, 

Jonah

[Jon Fawcett]

Thanks again, Jonah.

I really appreciate your interest here. I have appended an RDS file I generated depicting data that looks like the data I am using (although I have vastly more data than this, it seemed sufficient for demonstration purposes). I have modified the variable names to match the example I provided in my original post. Participant is nested within experiment. As a refresher:

experiment_id: Unique ID for specific experiments

participant_id: Unique ID for specific participants

condition: conditions within the experiment (you could think of these as treatment groups, coded as 0 or 1)

design: whether condition is manipulated between- or within-subjects (coded within or between)

decision: Binary outcome variable

Using these data, I am fitting:

dat = readRDS("sample_dat.rds")

glmer_model = stan_glmer(decision ~ condition * design + (condition | participant_id) + (condition | experiment_id), data = dat, family = binomial)

I would like to be able to determine whether the difference between the two conditions is larger for within-subject designs than between-subject designs (so the interaction).

[Jonah Sol Gabry]

Thanks Jon, that's helpful. Ok, so you have a logistic regression with a model formula of

decision ~ condition * design + (condition |participant_id) + (condition | experiment_id)

which corresponds to a linear predictor that is a function of:

- a global intercept ((Intercept))

- a global "main effect" of condition (condition)

- a global "main effect" of design (designwithin) 

- a global interaction term (condition:designwithin)

- shifts to the intercept and coefficient on condition specific to each participant_id (e.g. b[(Intercept) participant_id:1], b[condition participant_id:1])

- shifts to the intercept and coefficient on condition specific to each experiment_id (e.g. b[(Intercept) experiment_id:1], b[condition experiment_id:1])

and there will also be estimated covariance matrices for the multivariate normal prior on the participant-specific varying intercept/slope and for the experiment-specific varying intercept/slope.

That corresponds to what you're intending to do, correct?

Regarding the comparison of conditions across the two designs, I think preferable to focusing directly on the sign/magnitude of the interaction term would be to 

(1) use your posterior distribution to predict a bunch of outcomes from the posterior predictive distribution under the various possible scenarios you have  (you can use posterior_predict to get the 0/1 outcomes themselves or posterior_linpred with transform=TRUE to get the posterior probabilities of the outcomes)

and then 

(2) make the comparison of interest on the scale of the outcome via summaries of those predictions (rather than on the scale of the interaction coefficient). 

I think that's generally much more intuitive than directly interpreting interaction coefficients (especially when the linear predictor is transformed by a non-linear inverse-link function) and also the implications for the outcome tend be what's of practical interest. What do you think?

Best,

Jonah

P.S. It might also be interesting to do this both with the model that has the interaction term explicitly included and for the same model but without the interaction term so that you can compare that too. 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/mBMwLjFGB6k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jon Fawcett]

Thanks, Jonah!

Sorry for not replying immediately, I was caught up in the end-of-term rush and had fallen behind. Your post was very helpful. I have always undertaken comparisons in-text by using the coefficients to calculate an estimate for each level of a given condition (or range of values for continuous predictors), used the inverse transformation on the resulting posterior for each of those conditions and then subtract them to get an estimated difference in probability space (in the case of logistic regression). This was cumbersome, and from what I understand, posterior_linpred with transform=TRUE does the same thing in a single function call. Fantastic!

I have one more question for you, if I may. In addition to the analyses detailed above, I am also conducting analyses of the individual experiments, each with their own unique predictors. In one example, I have 96 participants in a within-subject design with 7 predictors. I have appended simulated data that matches the properties of my actual data very closely. The DV is still decision and the participant identifier is still participant_id, but the other variables are slightly different. condition and order are categorical (each with two levels) and they have been centred by defining these levels as -0.5 and 0.5; the remaining variables (choices, relprefc1, ici, bore, fat) are either continuous or ordinal and have been centred at their respective means. I analyze them using:

simdat2 = readRDS('sample_dat2.rds')
simmod2 = stan_glmer(decision ~ condition + order + choices + relprefc1 + ici + bore + fat + (condition|participant_id), data = simdat2, family = binomial, adapt_delta=.999, prior_covariance=decov(scale=.5))

I tweaked the prior on the RE covariance structure to be more constrained, because otherwise the random slope by condition takes on some pretty unreasonable values. In any event, using posterior_linpred as you suggested allows me to estimate the mean marginal probability for each level of condition – at the mean of all other variables – using the function call:

sim_conditions = posterior_linpred(simmod2, re.form=NA, transform=TRUE, newdata=expand.grid(condition=c(-.5,.5), order=0, choices=0, relprefc1=0, ici=0, bore=0, fat=0))
apply(sim_conditions, 2, function(x){return(quantile(x, probs=c(.025,.5,.975)))})

#                  1          2
#  2.5%  0.004359592 0.02508615
#  50%   0.032214613 0.10240261
#  97.5% 0.079338786 0.19666528

This is great! However, my question pertains to the fact that these estimates (M = .03 and .10) are quite different from the simple proportions obtained from the same data by calling (M = .08 and M = .22):

library(dplyr)
summarize(group_by(simdat2, condition), mean=mean(decision))

#  condition       mean
#      <dbl>      <dbl>
#1      -0.5 0.08333333
#2       0.5 0.21875000

I believe the reason for this is a consequence of modelling random effects so close to a boundary (i.e., 0). If I were to simulate the dataset as a whole, accounting for the random effects, as such:

y_rep <- posterior_predict(simmod2)
summary(apply(y_rep, 1, mean))

#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.05729 0.13020 0.15100 0.15120 0.17190 0.27080

The resulting distribution of means is spot on ([0.08 + 0.22] / 2 = 15). So, the question here is how to best interpret the output from posterior_linpred. Could this be viewed as the most credible median value for each condition – at the mean of every other variable – for an “average” participant? The implication would be that the difference between these estimates and the simple proportions would then be attributed to heterogeneity across subject. Or is there an even better way to compare the two conditions?

Cheers!

Jon

To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.