Mixed effect model specification for repeated measure placebo controlled trial

Question

I am analyzing gene expression data for a repeated measure placebo controlled clinical trial. Example dataset:

df <- data.frame(
  subject = c(paste0('Healthy', 1:8), paste0('PID', rep(1:12, 2))),
  disease = c(rep(FALSE, 8), rep(TRUE, 24)),
  treatment = c(rep('none', 8), rep(c('Placebo', 'Drug'), 12)),
  visit = c(rep('BL', 8), rep(c('BL', 'Day3'), each=12)),
  score = c(rep(NA,8), rep(c(0,0,30,30,50,50,70,70,90,90,100,100), 2)),
  gene = rnorm(32, 5)
)

Variables:

subject: IDs of either healthy controls or patients (PID prefix).

disease: is TRUE if subjects have the disease and FALSE if they are healthy controls.

treatment: 'none' for healthy and either 'Placebo' or 'Drug' for patients that will or did receive the corresponding treatment.

visit: is 'BL' for samples collected prior to treatment and 'Day3' for patient samples collected three days after treatment.

score: in an outcome that represents improvement in disease outcome and is measured 15 days after baseline. It could reasonably be modeled as an ordered factor but I decided to treat it as continuous for simplicity. This might not be appropriate as the maximum score is 100? Interpreting the polynomial terms confused me.

gene: is the expression of a sample gene that I want to ask questions about.

Model and Rationale:

mod <- lmerTest::lmer(gene ~ disease + treatment:visit + score:treatment:visit + (1|subject), data=df)

disease: I expect disease to independently affect the expression of gene (main effect).

treatment:visit: I expect to see an effect of being in the 'Drug' treated group only after treatment (visit Day3). I have omitted the main effect for treatment as it merely identifies healthy and randomization to Placebo or Drug groups. I have also omitted the main effect for visit as I don't expect gene expression to change on Day3 independent of treatment. Do these model specifications reflect my assumptions correctly?

score:treatment:visit: I am specifically interested in whether the expression of a given gene at both baseline and 3 days after treatment depends on disease improvement as measured by score on day 15. So, for example, gene expression could differ between subjects with disease before and after treatment in a way that depends on eventual score. Should the interaction be score:treatment:visit:disease to capture this?

Questions I want to address:

Does the expression of gene at baseline differ based on score at day 15?

Does the change in expression of gene between baseline and Day3 in subjects receiving the drug differ based on score at day 15?

Does the expression of gene at baseline differ between individuals with the disease and those without the disease?

How do I go about asking these questions using the model above (or an alternative specification if the above is incorrect)?

Thank you!

Answer 1

I think that you're trying to do too much with a single model. Things will be easier to understand and to explain to others if you break down the problem. Let's deal with the questions in reverse order.

Does the expression of gene at baseline differ between individuals with the disease and those without the disease?

To answer this question, it's simplest just to do a direct comparison of gene values between those 2 groups at baseline.

You might get slightly refined estimates with a mixed model that combines information from the 2 measurements taken on each of the individuals with disease. That could improve precision by a factor of about $\sqrt 2$ for the Disease group, but at some cost in extra assumptions and potential difficulty in explanation.

Although there are neither Day3 gene measurements nor a Drug category available for the Normal individuals, you could try something like the following. Apply the principle in this answer to deal with the data values that are impossible in your experimental design. Use 0/1 coding for disease (1 = Disease, 0 = Normal), treatment (1 = Drug, 0 for Normal and Placebo), and visit (1 = Day 3 visit, 0 for Normal and for the baseline visits). (Or if you prefer, False/True with False as the reference category.) That could give this mixed model (omitting the score variable for reasons explained below):

gene ~ disease + treatment*visit + (1|subject).

With the default treatment contrasts in R, the intercept in that model will be the model-estimated gene value for Normals (necessarily at a baseline visit and in the Placebo group with this coding). With that coding there are no further estimates made about treatment and visit effects for Normals, as they have 0 values for those. The random effect would model a distribution of intercept-value differences from that estimate (for Disease cases, taking into account the estimated fixed effect of Disease and any Drug-Placebo group difference at baseline).

Then: the disease coefficient would be the (Disease - Normal) group gene difference at the baseline visit (which could be considered an answer to this part of your question, although strictly this is the modeled difference of the Placebo group from Normal); the treatment coefficient would be the (Drug - Placebo) group gene difference at the baseline visit (ideally 0 if there are no baseline differences arising from treatment assignment); the visit coefficient would be the (Day3 - baseline) gene difference for the Placebo group (again, ideally 0 unless there is a Placebo effect on gene); the treatment:visit interaction would be the (Day3 - baseline) difference of gene values between the Drug group and the Placebo group (another result of interest).

I could see some objections, however, to modeling the random intercepts in the same way for the Normal and Disease individuals, even though the main effect of disease is accounted for in this modeling. For more complicated designs (like allowing for random effects of disease) you would have to make some assumptions about correlations between random intercepts and slopes, assumptions implicit but not always thought through in designing mixed-model formulas.

Does the change in expression of gene between baseline and Day3 in subjects receiving the drug differ based on score at day 15?

I (and I suspect many reviewers) would have a problem with causality in this type of model. Day 15 happens after the gene-expression measurements and the treatment. Using a value measured on Day 15 as a predictor, as in your proposed model, does not seem logical.

It would make more sense to treat score as the outcome and use the gene_difference as the predictor. (If you think that baseline gene values also matter, as implied by another of your questions, then you could use both baseline gene and gene_difference as predictors.) Restricting this analysis to those who received Drug, as you seem to propose, makes that pretty straightforward. You probably can get away with treating score as a continuous variable if its values cover a reasonable range.

Then you could combine information on all of your genes into a single model for score as the outcome variable. With a large number of genes you will have to model with methods like LASSO or ridge regression, but that's standard practice in gene-expression studies.

Does the expression of gene at baseline differ based on score at day 15?

In addition to the causality issue noted above, it's not clear from this how you intend to deal with the treatment issue here. With score as the outcome, I suppose you could restrict this analysis to the Placebo group and proceed as suggested above (using baseline gene values instead of the gene_difference values proposed for the Drug group). If there are no major differences in gene values for the Placebo group between baseline and Day3 then you could consider using the average of the baseline and Day 3 values to get potentially better estimates.

Finally, as a caution for future work, it's generally important to include all lower-level coefficients along with the interaction terms you are interested in. Your proposed model in this version of the question omits all single-predictor coefficients except disease and the 2-way interactions of treatment and visit with score. The full interaction model you proposed in an earlier version of the question would have been the correct form for a model with those 3 interacting predictors. See this page for extensive discussion.