adjustCovariate for a repeated measures two way ANOVA-like setup

[lars-work-sund]

I have an experimental setup where each animal has been assigned to an experimental group and then measured twice, once before treatment and once after a treatment, I have an example of the pData table below. I am interested in comparing after treatment to before treatment within the groups, while adjusting for the effect of animal. Ideally I would also like to find an interaction, but that does not seem possible.

When running the dmrseq function, I get the following error:
Detecting candidate regions with coefficient larger than 0.1 in magnitude. ...Chromosome chr1: Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 't': invalid 'k' argument Calls: test_covariate ... t -> tcrossprod -> backsolve -> .handleSimpleError -> h

Is what I'm trying to do possible with dmr-seq?

SampleID	Treatment	Timepoint	Animal
S1	Ctrl	Before	animal1
S2	Ctrl	After	animal1
S3	Ctrl	Before	animal2
S4	Ctrl	After	animal2
S5	Treat	Before	animal3
S6	Treat	After	animal3
S7	Treat	Before	animal4
S8	Treat	After	animal4

Best

[kdkorthauer]

Hi @lars-work-sund,

Can you provide me with the call to the dmrseq function that throws the error?

Looking at your covariate table, it's possible you're trying to estimate effects that you can't estimate with that experimental design. If that's the case, then adding a more informative error message would be helpful. But I can't be sure until I know what model you're trying to fit. Thanks!

[lars-work-sund]

Of course, I created a column called group by pasting Treatment and Timepoint together, then I ran the command dmrseq(bs, testCovariate = "group", adjustCovariate = "Animal")

Looking at it now I realize it only looks for a difference across all groups. Are there any plans to support more complicated experimental setups? Also the ability to test for interactions or specific contrasts in dmrseq would be really nice. It is by far the best tool for analyzing bisulfite data I have found.

[kdkorthauer]

Hi @lars-work-sund ,

Thanks for the additional information. When you create a new variable by pasting treatment and timepoint together, this is effectively creating a 4 level factor, with 2 replicates (animals) in each group. This would be perfectly fine to pass into testCovariate on its own, and would look for any effect of treatment or time. However, when you try to add animal as the adjustCovariate, this is not an effect that is possible to estimate with your experimental design, because there is only one observation per combination of treatment*timepoint x animal. Mathematically, there are columns of the design matrix that you can obtain via a linear combination of other columns. Intuitively, the effect of animal is confounded with the effect of treatment (each animal is only in control or treatment, not both).

In a linear modeling sense, typically the animal effect would be better represented by a random effect than a fixed effect, but linear mixed models are outside the scope of dmrseq (as is testing specific contrasts, since each one needs its own permutation scheme - instead of specific contrasts, you can set up your input variables so that the testCovariate represents the effect of interest and rerun).

I can suggest a couple of alternative strategies for your situation, but I'd advise you consult with a statistician who has knowledge of your experimental system and biological questions (my advice is only mathematical, because this is all the information I have). One option is to model a new response variable that reflects your paired design - difference in DNAm level before and after for each animal, then input treatment as your testCovariate. Another option is to do as you did before, but not adjust for animal. As you point out, this tests for any difference amongst the groups, so if you're interested in looking at effect of time controlling for treatment (or the reverse), just input time as the testCovariate and treatment as the adjustCovariate (or the reverse).

Hope that helps!

[kdkorthauer]

Pushed an addition to the devel branch that should give a more informative error message for putting in a design that is not full rank.