The sleuth example workflows generally have experimental designs of unpaired data with only 2 experimental conditions.
so <- sleuth_prep(s2c, ~ treatment) so <- sleuth_fit(so) so <- sleuth_fit(so, ~1, 'reduced') so <- sleuth_lrt(so, 'reduced', 'full')
Ensuring that fitted models are aware of paired data seems relatively straightforward, we simply include a biologicalreplicate variable in the model and retain it in the reduced model.
so <- sleuth_prep(s2c, ~ biologicalreplicate + treatment) so <- sleuth_fit(so) so <- sleuth_fit(so, ~ biologicalreplicate, 'reduced') so <- sleuth_lrt(so, 'reduced', 'full')
But what if our treatment column can have three different states (A, B and C), instead of just two? The code stays the same and dummy variables are automatically generated. But the likelihood ratio test now tests which transcripts are affected by any of the two conditions B and C. If we want to examine the effects of a specific condition, we can use the Wald test.
so <- sleuth_prep(s2c, ~ biologicalreplicate + treatment) so <- sleuth_fit(so) so <- sleuth_fit(so, ~ biologicalreplicate, 'reduced') so <- sleuth_wt(so, which_beta = 'treatmentB', 'full') so <- sleuth_wt(so, which_beta = 'treatmentC', 'full')
If the treatment factor in the s2c object is properly set up, this lets us compare treatment A to B and A to C.
But how do we compare B to C, short of redoing the analysis and reordering the factor to take B as the baseline, instead of A?
In the same vein, is it possible to make these comparisons using the LRT, instead of the Wald test considering the reservations regarding the Wald test expressed by Harold in this Biostars question? The only way I can think of to do this would be to run three different analyses, using restricted datasets that only contain the data relevant to the conditions of interest used in the three different models. We would build a 'full' and 'reduced' model for a dataset consisting of samples treated with A and B, then again for B and C and again for A and C and perform the LRT on each pair of models in turn. But there's no way that can be accurate, surely...
I appreciate any help. If my google-fu has missed workflows where this is explained, I apologize and will gladly settle for a link to those workflows.