[BioC] Defining the model matrix for a multi-factorial design with interactions - edgeR

Ryan rct at thompsonclan.org
Wed Jan 29 19:12:56 CET 2014


Hi Marco,

Note that your models 2 & 3 are exactly the same model. The expression 
in your formula for model 3:

Genotype:Load:Time

is more or less equivalent to:

factor(paste(design$Genotype, design$Load, design$Time, sep=":"))

except that the factor levels might possibly be in a different order. As 
for your question about doing the same test with Models 1 & 2, I believe 
that both models are equivalent parametrizations of the same design, so 
you should get (nearly) identical dispersion estimates, and performing 
equivalent tests with the two models should give you the same p-values. 
If you are not getting the same p-values, then you have not properly set 
up the contrast for Model 1. This isn't surprising, because a 3-way 
interaction model is really confusing to work with, and I can't even 
tell you what the correct contrast would be. I would recommend that you 
stick with Model 2 (or 3), where each coefficient has a direct 
biological interpretation as the estimated expression level of a group. 
Designing contrasts for this parametrization is much easier, in my 
opinion. I think you would agree, since you got this contrast correct in 
your example code.

-Ryan

On 1/29/14, 7:11 AM, Marco [guest] wrote:
> Hi all!
>
> I have a question about the way to define the model matrix for a multi-factorial design in which I want to model all the interactions. In particular, I have a 3x3x2 factorial design (3 factors) and three biological replicates, so I have 3x3x2x3 = 54 sample.
> The three factors are:
> - GENOTYPE with the three levels F, G and P
> - LOAD with levels HIGH, MEDIUM and LOW
> - TIME with levels H and PH
>
> So this is the design matrix:
>
>> design
>         Time   Load Genotype
> FHH_1     H   High        F
> FHPH_1   PH   High        F
> FHH_2     H   High        F
> FHPH_2   PH   High        F
> FHH_3     H   High        F
> FHPH_3   PH   High        F
> FMH_1     H Medium        F
> FMPH_1   PH Medium        F
> FMH_2     H Medium        F
> FMPH_2   PH Medium        F
> FMH_3     H Medium        F
> FMPH_3   PH Medium        F
> FLH_1     H    Low        F
> FLPH_1   PH    Low        F
> FLH_2     H    Low        F
> FLPH_2   PH    Low        F
> FLH_3     H    Low        F
> FLPH_3   PH    Low        F
> GHH_1     H   High        G
> GHPH_1   PH   High        G
> GHH_2     H   High        G
> GHPH_2   PH   High        G
> GHH_3     H   High        G
> GHPH_3   PH   High        G
> GMH_1     H Medium        G
> GMPH_1   PH Medium        G
> GMH_2     H Medium        G
> GMPH_2   PH Medium        G
> GMH_3     H Medium        G
> GMPH_3   PH Medium        G
> GLH_1     H    Low        G
> GLPH_1   PH    Low        G
> GLH_2     H    Low        G
> GLPH_2   PH    Low        G
> GLH_3     H    Low        G
> GLPH_3   PH    Low        G
> PHH_1     H   High        P
> PHPH_1   PH   High        P
> PHH_2     H   High        P
> PHPH_2   PH   High        P
> PHH_3     H   High        P
> PHPH_3   PH   High        P
> PMH_1     H Medium        P
> PMPH_1   PH Medium        P
> PMH_2     H Medium        P
> PMPH_2   PH Medium        P
> PMH_3     H Medium        P
> PMPH_3   PH Medium        P
> PLH_1     H    Low        P
> PLPH_1   PH    Low        P
> PLH_2     H    Low        P
> PLPH_2   PH    Low        P
> PLH_3     H    Low        P
> PLPH_3   PH    Low        P
>
>
> I am interested in finding all possible interactions between the three factors, so I want to model also the three way interactions GENOTYPE:LOAD:TIME. I have already read the edgeR vignette so I was wondering if the best way to define the model matrix is by defining it with a formula or by defining each treatment combination as a group.
> I'm trying to understand which are the differences between the two alternatives, so I tried to estimate both the models:
>
> Model 1)
>    mydesign1<-model.matrix(~Genotype * Load * Time, data=design)
>    D1 <- estimateGLMCommonDisp(D, mydesign1) # D is the DGEList of the normalized counts
> 	D1 <- estimateGLMTrendedDisp(D1, mydesign1)
> 	fit1 <- glmFit(D1, mydesign1)
>
> The estimated coefficients are:
> 	> colnames(fit1)
> 	 [1] "(Intercept)"                 "GenotypeG"
> 	 [3] "GenotypeP"                   "LoadLow"
> 	 [5] "LoadMedium"                  "TimePH"
> 	 [7] "GenotypeG:LoadLow"           "GenotypeP:LoadLow"
> 	 [9] "GenotypeG:LoadMedium"        "GenotypeP:LoadMedium"
> 	[11] "GenotypeG:TimePH"            "GenotypeP:TimePH"
> 	[13] "LoadLow:TimePH"              "LoadMedium:TimePH"
> 	[15] "GenotypeG:LoadLow:TimePH"    "GenotypeP:LoadLow:TimePH"
> 	[17] "GenotypeG:LoadMedium:TimePH" "GenotypeP:LoadMedium:TimePH"
>    
>    
> Model 2)
>    Group <- factor(paste(design$Genotype, design$Load, design$Time, sep="."))
>    # I have 18 unique levels in Group (54 samples / 3 biological replicates)
>
>    mydesign2<-model.matrix(~0+Group)
>    colnames(mydesign2) <- levels(Group)
> 	
> 	D2 <- estimateGLMCommonDisp(D, mydesign2)
> 	D2 <- estimateGLMTrendedDisp(D2, mydesign2)
> 	
> 	fit2<-glmFit(D2, mydesign2)
> 	> colnames(fit2)
> 	 [1] "H.High.F"    "H.High.G"    "H.High.P"    "H.Low.F"     "H.Low.G"
> 	 [6] "H.Low.P"     "H.Medium.F"  "H.Medium.G"  "H.Medium.P"  "PH.High.F"
> 	[11] "PH.High.G"   "PH.High.P"   "PH.Low.F"    "PH.Low.G"    "PH.Low.P"
> 	[16] "PH.Medium.F" "PH.Medium.G" "PH.Medium.P"
>    
> Now, to understand the differences between the two models, my question was this:
> if I define the contrast:
>
>    contrast1<-c(rep(0,14), 1, 0, -1, 0) # "GenotypeG:LoadLow:TimePH" vs "GenotypeG:LoadMedium:TimePH"
>
> on the first model and I make a LR test, will I obtain the same result as making a LR test on the contrast:
>
>    contrast2<-makeContrasts(G.LowvsMedium.PH = PH.Low.G - PH.Medium.G, levels=mydesign2)
>
> on the second model?
> This is the code I used to make the tests:
>
>    lrt1.G.LowvsMedium.PH<-glmLRT(fit1, contrast=contrast1)
>    summary(decideTestsDGE(lrt1.G.LowvsMedium.PH))
>    
> and
>
>    lrt2.G.LowvsMedium.PH<-glmLRT(fit2, contrast=contrast2[," G.LowvsMedium.PH"])
>    summary(decideTestsDGE(lrt2.G.LowvsMedium.PH))
>
> The answer to the question is no, because I've tried to do it on my data and the results are different between the two tests. But why? Am I doing two completely different tests? The two parameters of the two models aren't saying me the same thing?
>
> In a moment of madness :) I tried to estimate another model, with only the three-level interactions between the factors.
>
> Model 3)
>    mydesign3<-model.matrix(~0+Genotype:Load:Time, data=design)
>
>    D3 <- estimateGLMCommonDisp(D, mydesign3)
>    D3 <- estimateGLMTrendedDisp(D3, mydesign3)
>
>    fit3<-glmFit(D3, mydesign3)
>    > colnames(fit3)
>    [1] "GenotypeF:LoadHigh:TimeH"    "GenotypeG:LoadHigh:TimeH"
>    [3] "GenotypeP:LoadHigh:TimeH"    "GenotypeF:LoadLow:TimeH"
>    [5] "GenotypeG:LoadLow:TimeH"     "GenotypeP:LoadLow:TimeH"
>    [7] "GenotypeF:LoadMedium:TimeH"  "GenotypeG:LoadMedium:TimeH"
>    [9] "GenotypeP:LoadMedium:TimeH"  "GenotypeF:LoadHigh:TimePH"
>    [11] "GenotypeG:LoadHigh:TimePH"   "GenotypeP:LoadHigh:TimePH"
>    [13] "GenotypeF:LoadLow:TimePH"    "GenotypeG:LoadLow:TimePH"
>    [15] "GenotypeP:LoadLow:TimePH"    "GenotypeF:LoadMedium:TimePH"
>    [17] "GenotypeG:LoadMedium:TimePH" "GenotypeP:LoadMedium:TimePH"
>    
> Well, this model gave me the same parameter estimates as the model 2 (each treatment combination as a group) and, obviously, the same results in terms of DEG. But is it statistically correct to estimate a model with three-level interactions without any principal effect?
>
> Thanks in advance for any replies!
>
> Marco
>
>   -- output of sessionInfo():
>
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=Italian_Italy.1252  LC_CTYPE=Italian_Italy.1252
> [3] LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C
> [5] LC_TIME=Italian_Italy.1252
>
> attached base packages:
> [1] splines   parallel  stats     graphics  grDevices utils     datasets
> [8] methods   base
>
> other attached packages:
>   [1] EDASeq_1.8.0         ShortRead_1.20.0     Rsamtools_1.14.2
>   [4] lattice_0.20-24      Biostrings_2.30.1    GenomicRanges_1.14.4
>   [7] XVector_0.2.0        IRanges_1.20.6       Biobase_2.22.0
> [10] BiocGenerics_0.8.0   edgeR_3.4.2          limma_3.18.9
> [13] aroma.light_1.32.0   matrixStats_0.8.14   yeastRNASeq_0.0.10
>
> loaded via a namespace (and not attached):
>   [1] annotate_1.40.0      AnnotationDbi_1.24.0 bitops_1.0-6
>   [4] DBI_0.2-7            DESeq_1.14.0         genefilter_1.44.0
>   [7] geneplotter_1.40.0   grid_3.0.2           hwriter_1.3
> [10] latticeExtra_0.6-26  R.methodsS3_1.6.1    R.oo_1.17.0
> [13] RColorBrewer_1.0-5   RSQLite_0.11.4       stats4_3.0.2
> [16] survival_2.37-6      tools_3.0.2          XML_3.98-1.1
> [19] xtable_1.7-1         zlibbioc_1.8.0
>
> --
> Sent via the guest posting facility at bioconductor.org.
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor



More information about the Bioconductor mailing list