I have collected a data-set and run into some problems on how to correctly use a linear mixed model to analyse my data. The data collection (simplified): I have 24 family groups (“family”) of 2 individuals each that were treated with A, B, C or D (“treat1”, e.g. group 1-6 was treated with A, group 7-12 was treated with B etc.; total n = 48). Then, I separated all individuals and one of the two individuals per family was treated with E, while the other was not as a control (“treat2”). As dependent variable, I measured various genes’ expression levels (which I log-transformed subsequently). The main question I'd like to address is: does “treat1” (A,B,C,D) affect the individuals’ response (gene expression) to “treat2” (E)? Assuming log-transformed gene expression is normally distributed, I’d use the following model for each gene individually (from the nlme package): lme(gene_expression~treat1*treat2, random = 1|family, method=”REML”)

I conceptually assume that “family” should not affect the interaction of treat1 and treat2, which is why I only have a random intercept here, but I wonder if the model can actually deal with the experimental design, specifically the fact that each family only reflects one level of treat1, but two levels of treat2. While this model mostly appears to produce reasonable results, I occasionally get a “singular convergence (7)” error message, and I am not sure why that is (the response data looks fine as genes with many 0 expression levels are not considered).

It would be wonderful to get somebody’s opinion on this, specifically if the random intercept is fine or maybe a random slope has to be added, and on what might cause the error message. Thank you so much!

  • 1
    $\begingroup$ Have you tried using lmer() from the lme4 package instead? Sometimes that can avoid this error; see this thread. lme() allows for some correlation structures that aren't available in lmer(), but a simple random intercept isn't a problem. Also, how many genes are involved and how are you handling multiple-comparison issues? $\endgroup$
    – EdM
    Dec 9, 2021 at 21:15
  • $\begingroup$ Thank you and sorry for the delayed reply - will try lme4. I have ~16k genes and what I would do is note the p-value (etc.) per gene and variable and then, per variable, do a fdr correction. However, I have noticed that this is probably too conservative (?), because I have variables that are not significantly affecting y in any gene after fdr correction, but when I plot histograms of raw p-values, I observed an over-representation of small p-values. Assuming I test 16k genes, I'd expect approx. 0.05*16k=800 false positives, but before fdr correction, I have 2k genes <0.05, but none after fdr. $\endgroup$ Feb 3, 2022 at 7:49

1 Answer 1


With 16000 genes you will be better off using software packages that are designed to handle large-scale gene-expression data. The DESeq2 and edgeR Bioconductor packages were designed for working with RNAseq data, and the venerable limma package originally designed for spotted microarrays can work with such data, too.

Your gene-by-gene analysis does not model variance as a function of gene expression level, unlike the packages cited above. For example, the DESeq2 package performs negative binomial modeling of RNA-seq counts as a function of both sample-specific and gene-specific factors. That provides better pooled error estimates of variances to use for differential expression analysis than does the constant variance on a log scale implicit in your approach, potentially improving power to detect true changes. Those packages can help deal with outliers and handle the multiple-comparisons problem. They accept design matrices in ways that should accommodate your experimental design.

Model matrix example

The main question I'd like to address is: does “treat1” (A,B,C,D) affect the individuals’ response (gene expression) to “treat2” (E)?

Section 3.5 of the edgeR Users Guide has a design addressing essentially the same question as yours. Each individual Patient, having one of 3 Disease types, received both a control and a hormone Treatment. The question there is whether the Disease affects the response to Treatment, like your interest in whether "treat1" affects the response to "treat2"; it has pairing like yours.

To get a corresponding design matrix for your study, replace the User Guide's 9-level Patient with your 24-level "family"; its 3-level Disease with your 4-level "treat1"; its 2-level Treatment with your 2-level "treat2":

fam <- factor(rep(1:24,each=2))
trt1 <-factor(rep(LETTERS[1:4],each=12))
trt2 <- relevel(factor(rep(c("None","E"),24)),ref="None")
AE <- trt1=="A" & trt2 =="E"
BE <- trt1=="B" & trt2 =="E"
CE <- trt1=="C" & trt2 =="E"
DE <- trt1=="D" & trt2 =="E"
design <- model.matrix(~fam)
design <- cbind(design,AE,BE,CE,DE)

That accomplishes with pairing the critical part of what your mixed model was doing. That section of the User Guide then shows how to use the resulting model to find genes that respond differentially to combinations of conditions.

The section of the DESeq2 vignette on "Group-specific condition effects, individuals nested within groups" suggests a more efficient model-matrix coding inheriting from that edgeR method. As the particular "family" names aren't themselves important and no "family" is involved in more than 1 level of "treat1", you can set up an "ind.n" factor that just annotates the 6 separate families within each level of "treat1." Then your model matrix could be based on the formula ~treat1 + treat1:ind.id + treat1:treat2. The vignette goes on to illustrate how to get comparisons of interest.

I haven't carefully thought through the differences between those two suggestions. The point is that these standard packages should be able to answer your fundamental question.

A comment on an earlier version of this answer suggests that you might have additional covariates. If so, the edgeR-recommended model matrix only has 28 columns and the DESeq2 recommendation only has 12, allowing you to add columns for some additional covariates (if they aren't linearly dependent on the columns already included).

If you do need to use a mixed model, you might need to consider a two-step approach as in Trabzuni et al., Bioinformatics, Volume 30, Issue 11, 1 June 2014, Pages 1555–1561, which combined mixed modeling with subsequent finite mixture modeling to separate differentially expressed from non-differential transcripts.

  • $\begingroup$ Thanks for the suggestions. DeSeq2 was my first choice, but it seems my experimental design is too complex for what the package can handle. I also want to stay away from reducing the complexity of gene expression to "DE" and "nonDE", and rather be able to interpret slopes and meaningful interactions among variables and their effects on slopes, which is information I could not get out of DeSeq2. I also played around with limma before, but I also could not really find a suitable model. I will look at Trabzuni et al., which looks interesting but probably I am not smart enough to implement it :-( $\endgroup$ Feb 3, 2022 at 14:28
  • $\begingroup$ @RalfSchneider I'm a bit surprised that DESeq2 didn't work. Rather than coding interactions in a model matrix it can be less error prone to specify a single multi-level categorical predictor and then do post-hoc contrasts; see the vignette. Also, it's not clear what you mean by "interactions among variables and their effects on slopes"; I didn't see those in the question. If there are additional predictors in your model, those might have been the source of original problems with mixed models. $\endgroup$
    – EdM
    Feb 3, 2022 at 17:04
  • $\begingroup$ @RalfSchneider the edgeR Users Guide has an example almost identical to your "main question." I've edited the answer to show how to apply that approach to construct a model matrix for the situation described in your question. The pairing accomplishes pretty much what you get from a mixed model; using DESeq2 or edgeR avoids the downside of using mixed models that assume constant variance on the log scale. $\endgroup$
    – EdM
    Feb 3, 2022 at 20:02
  • $\begingroup$ Thanks so much! I'm looking into matrix modeling and the DeSeq2 vignette now. However, the actual experiment has another factor (treat3, 2 levels) and treat1 may actually be better treated as 2 separate factors with 2 levels each (i.e. the parental animals were treated or not: father (treat1.1): yes/no; mother (treat1.2): yes/no), n=94 in total. Also, interactions are central to the question, e.g., which genes respond to treat2 and/or 3 differently depending on treat1.1 and 1.2? My mixed model gives me coefs and p for all of this, but seems very conservative --> I'd love to use DeSeq2... $\endgroup$ Feb 10, 2022 at 8:43
  • $\begingroup$ @RalfSchneider the model matrix used by DESeq2 is as in standard R modeling. Play with small, maybe simulated, data sets. If a model matrix works for a single gene, it should work for the whole set. Look carefully at the DESEq2 vignette about interactions, because much can be obtained more simply by contrasts. You might use DESeq2 to get variance-stabilized transformations of counts for other analyses like mixed modeling, although the authors recommend using DESeq2 directly for differential expression. $\endgroup$
    – EdM
    Feb 10, 2022 at 20:52

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.