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)


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!

  • $\begingroup$ Could you please say some more about what the "score" variable represents in your model? The way you describe it, its (post-pre) difference seems more like an outcome than a predictor. And is this just one gene whose expression you're evaluating, or is this a multi-gene study? $\endgroup$
    – EdM
    Commented Apr 24, 2020 at 22:02
  • $\begingroup$ Are the gene measurements only done post, or do you have post-pre differences for those too? $\endgroup$
    – EdM
    Commented Apr 24, 2020 at 22:52
  • $\begingroup$ I added some clarifications to the original question. Thank you! $\endgroup$ Commented Apr 25, 2020 at 18:41

1 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.

  • $\begingroup$ Thank you @EdM - a couple of questions. 1) For comparing baseline disease vs healthy doesn't it make sense to also model day3 disease so that I can include the subject random effect to model the within-individual correlation? This is partly an exercise for me to understand linear models so thank you for the comments/suggestions! $\endgroup$ Commented Apr 26, 2020 at 12:09
  • 1
    $\begingroup$ @alexvpickering for comparing gene for baseline disease vs healthy, you won't have repeated measures on the healthy cases so it's not clear how much is gained by adding mixed modeling of the disease cases. I suppose a model that included treatment and visit, with their interaction, could provide somewhat more precise estimates of true "baseline" gene values for the disease cases, but I suspect that the advantage will be minimal and you will have to explain a more difficult model to your audience. Try both ways and see. Modeling is a process, not a one-off effort. $\endgroup$
    – EdM
    Commented Apr 26, 2020 at 15:41
  • $\begingroup$ @alexvpickering I've added a potential mixed-effects model to the answer, which at least should help you think about some of the issues with mixed models and interaction terms. Note that even this relatively simple model estimates an intercept, a distribution of intercept values among subjects, three single-predictor coefficients, and an interaction coefficient. As I count it, you have only (2*2*12 + 8) = 56 gene observations in your example, so if that's the case in your study that might already be close to over-fitting. $\endgroup$
    – EdM
    Commented Apr 27, 2020 at 15:33

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.