I am working with some biological data where I am trying to compare 5 categories of genes by expression which occur at different loci (locations on the chromosomes). I have samples from different tissue categories and 3 technical replicates per sample. Since I am interested in comparing gene categories by locus, I considered using 2-way repeated measures ANOVA, with the categories as repeated measurements, the locus as the subject ID, and the tissues as the second factor. However, I ran into two problems:

  1. the issue of how to handle the technical replicates
  2. some loci have fewer than 5 categories of genes i.e. non-random missing values for the repeated measures

Therefore I used a linear mixed model, treating the locus as a random effect using the lme function in R as follows:

mod = lme(expression~(tissue_type)*(gene_category), random = ~1|locus, data=df)

followed by the Anova() function from the car package to check significance of the interaction.

Regarding the posthoc analysis, I was interested in comparing gene categories pairwise for each tissue type so as to ascertain if categories differ significantly (for each tissue) and if so, which one dominates in gene expression. This is where I am unsure of how to proceed. The usual recommendation is to use the emmeans package, but I am unsure of how this handles the missing categories. Suppose I run

mod.emm <- emmeans(mod, ~ tissue_type*gene_category)

pairwise = pairs(mod.emm, simple="gene_category",adjust="tukey")

I assume this gives me a pairwise comparison of gene categories using a pooled t-test. However, I want paired comparisons by locus ID, even if that means filtering down to only those loci that have existing pairs for the two selected gene categories. Is is inappropriate to do a paired t-test this way? The following post has me confused. Any advice would be greatly appreciated!

link to post

  • $\begingroup$ Your pairwise object gives you pairwise comparisons of gene categories separately for each tissue type. Meanwhile, locus is a random effect, which means it is out of the picture as far as emmeans goes (or rather, the variations among loci are part of the SE calculations). A random effect means that you are considering a sample from a population of loci, and your inferences apply to that whole population. If you want estimates for particular loci, you should have locus in the model as a fixed effect. $\endgroup$
    – Russ Lenth
    Nov 1, 2023 at 2:48
  • $\begingroup$ @RussLenth I am not interested in estimates for particular loci. I am only interested in comparing genes that occur together so I get an idea of how expression varies by gene category for co-occuring genes - particularly if one gene category dominates another. I initially considered just filtering down the data to the loci that have all 5 gene categories present but this seems to miss out on info in much of the data. Would it be misleading to do pairwise gene category comparisons using a paired t-test - say filtering to those loci that have non missing values for those categories? $\endgroup$
    – Merry
    Nov 1, 2023 at 2:55
  • $\begingroup$ With the usual medical analogy - let's say I have patients whose blood pressure is measured after undergoing five different treatments (A,B,C,D,E) each under 3 different conditions. Some patients don't show up for some of the treatments, but I am interested in comparing A vs B, C vs D etc. under each of the three conditions. Is it incorrect use a paired t-test to make patient-wise comparisons for A vs B by only considering patients who showed up for those two treatments? $\endgroup$
    – Merry
    Nov 1, 2023 at 3:14
  • 2
    $\begingroup$ I do not recommend subsetting the data. The mixed model you fitted has the information already in it as to which subjects receive which treatments and which don't., and the estimates and their standard errors are derived accordingly. Trying to micromanage how it does that only gets you into trouble and reduces your power. Just fit the model, make sure it's appropriate and fits well;, then test the comparisons that are of interest. $\endgroup$
    – Russ Lenth
    Nov 6, 2023 at 21:28
  • 1
    $\begingroup$ And don't use pooled t tests. Use the emmeans() results. $\endgroup$
    – Russ Lenth
    Nov 6, 2023 at 21:31


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.