2
$\begingroup$

I'm trying to figure out the proper way to run a Linear Mixed-Effects Model. I am looking at the length of a certain gene modification for 3 genes and how they change upon treatment (there's one control and one experimental condition). Measurement of this length comes from 4 replicates but, due to experimental limitations, the replicates are unbalanced (some have 1000 observations, others have 5000).

There was a similar experiment published and in the paper they used a Linear Mixed-Effects Model (R package lmerTest with the lmer function). As I understand it, this test would account for the unbalanced groups and allow me to use all of my many (MANY) data points. One option is to simply subset my data by gene and then input into lmer like this:

lmerTest::lmer(formula = log(length) ~ treatment + (1|replicate), data = df %>% filter(gene = X))

And this produces a single coefficient and p value that make sense. Then I collect the results and correct the p values for multiple comparisons afterwards.

I think it would be better to actually test all my genes of interest in the linear mixed effects model, but I'm not sure how to do it or how to interpret the results. I think I want an interaction term between gene and treatment.

Here's what I tried:

lmerTest::lmer(formula = log(length) ~ treatment*gene + (1|replicate), data = df %>% filter(gene = X))

but this produces a results table that I'm not sure I understand:

Fixed effects:
                                  Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                      4.310e+00  1.143e-01  3.489e+00  37.691 1.11e-05 ***
conditionDRUG                    2.250e-01  4.939e-02  2.256e+04   4.554 5.28e-06 ***
gene_nameA                       1.286e-01  3.296e-02  2.256e+04   3.901 9.60e-05 ***
gene_nameB                      -7.296e-02  3.556e-02  2.256e+04  -2.052   0.0402 *  
conditionDRUG:gene_nameA        -1.309e-01  5.139e-02  2.256e+04  -2.548   0.0108 *  
conditionDRUG:gene_nameB        -1.175e-01  5.449e-02  2.256e+04  -2.157   0.0310 *  

So I think:

  • The intercept is how much of the variance the model can't account for (a bucket for "other stuff")
  • conditionDRUG is telling me the effect of the drug on ALL genes
  • gene_nameA is gene A compared with the missing gene, C, which I guess is the reference?
  • gene_nameB same as above but for B
  • conditionDRUG:gene_nameA these last two I don't quite understand. Are the relative to the missing interaction term or what?

Ultimately I want to know, for each gene does the treatment change its length.

I'm happy to provide a toy dataframe if someone wants to run it, but the code (I think) runs fine, I just want to understand the output. Also apologies if this would be better on the statistics stackOverflow (I don't have any rep there to post)

$\endgroup$
3
  • 1
    $\begingroup$ You tagged me in the other thread. Sorry I don't have the time to provide a full answer, but I think you might find this link useful: marginaleffects.com/articles/experiments.html The context is simpler, but the same commands should work in analogous fashion in a mixed effects model with interactions like the one you are proposing. This may also be somewhat relevant: marginaleffects.com/articles/lme4.html $\endgroup$
    – Vincent
    Oct 7, 2023 at 12:16
  • 1
    $\begingroup$ Is there a typo in the model output ? conditionDRUG:gene_nameA appears twice ... $\endgroup$ Oct 10, 2023 at 11:06
  • 1
    $\begingroup$ yes that's a typo -- i fixed it $\endgroup$
    – Max F
    Oct 10, 2023 at 17:27

1 Answer 1

2
$\begingroup$

To summarize: condition (aka treatment) has two levels, "DRUG" and "NO DRUG". And gene (aka gene_name) has three levels, "A", "B" and "C". The "NO DRUG" and "C" levels are the reference levels for the condition and gene variables, respectively.

And the research question is:

For each gene does the treatment change its length? That is, estimate the average difference in log length due to treatment, for each of the three genes A, B and C.

Since both fixed effects are categorical, the model coefficients have straightforward interpretation:

  • The intercept is how much of the variance the model can't account for (a bucket for "other stuff") This is the population-level E(log(length) | condition = "NO DRUG", gene = "C"). That is, the average log length of gene C when no treatment is applied.
  • conditionDRUG is telling me the effect of the drug on ALL genes This is the population-level difference E(log(length) | condition = "DRUG", gene = "C") - E(log(length) | condition = "NO DRUG", gene = "C"). This is, the change in the average log length of gene C due to the treatment. It's the answer to the research question for gene C.
  • gene_nameA is gene A compared with the missing gene, C, which I guess is the reference? This is the population-level difference E(log(length) | condition = "NO DRUG", gene = "A") - E(log(length) | condition = "NO DRUG", gene = "C"). That is, the difference between the average log length of un-treated A gene and untreated C gene.
  • gene_nameB same as above but for B
  • conditionDRUG:gene_nameA these last two I don't quite understand. Are the relative to the missing interaction term or what? Yes, these terms are interactions defined as population-level differences compared to the intercept. So E(log(length) | condition = "DRUG", gene = "A") - E(log(length) | condition = "NO DRUG", gene = "C"). This interaction doesn't answer the research question directly because both condition and gene change while we want to condition on gene and estimate the difference due to treatment.

Here is where it's helpful to express the research question in terms of pairwise comparisons of marginal effects and compute those comparisons with either marginaleffects or emmeans. This approach scales effortlessly with larger and more complex models.

fit <- lmer(
  log(length) ~ treatment * gene + (1 | replicate),
  data = df
)

# With emmeans:
pairs(emmeans(
  fit, ~ treatment | gene
))

# With marginaleffects:
avg_comparisons(
  fit, variables = "treatment", by = "gene"
)
$\endgroup$

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.