Hello,
I am seeking confirmation on manually setting contrasts for a 3-factor design in DESeq2 with multiple interaction terms
I have a 3-factor experiment, with each factor containing 2 levels:
Type: Juvenile vs Sediment
Shore: Inshore vs Offshore
Sector: North vs Central
I am interested in the following contrasts:
- Average main effects of Type, Shore and Sector across all the levels.
- Interactive effect of Type:Shore, specifically:
- Juvenile: Inshore vs Offshore
- Sediment: Inshore vs Offshore
- Interactive effect of Type: Sector, specifically:
- Juvenile: North vs Central
- Sediment: North vs Central
- Interactive effect of Type vs Shore:Sector, specifically:
- Juvenile: Inshore North vs. Inshore Central
- Juvenile: Inshore North vs. Offshore Central
- Juvenile: Inshore North vs. Offshore North
- Juvenile: Offshore Central vs. Inshore Central
- Juvenile: Offshore North vs. offshore Central
- Juvenile: Offshore North vs. Inshore Central
5. 1-6: Sediment: with same 6 comparisons above
My experimental design formula is as follows:
~Type*Shore*Sector
Base levels for all factors are: Juvenile, Inshore, Central
And my Model Matrix column names outputted with resultsNames function:
resultsNames(varstab_geomLove_T.S.S) [1] "Intercept" [2] "Type_Sediment_vs_Juvenile" [3] "Shore_Offshore_vs_Inshore" [4] "Sector_North_vs_Central" [5] "TypeSediment.ShoreOffshore" [6]"TypeSediment.SectorNorth" [7] "ShoreOffshore.SectorNorth" [8]"TypeSediment.ShoreOffshore.SectorNorth"
|
Effect |
Function |
Ok? |
|
Main effect of Type (aka Juvenile vs Sediments) |
results(varstab2, contrast=c("Type", "Juvenile", "Sediment")) |
yes |
|
Main effect of Shore (aka Inshore vs Offshore) |
results(varstab2, contrast=c("Shore", "Inshore", "Offshore")) |
yes |
|
Main effect of Sector (aka North or Central) |
results(varstab2, contrast=c("Sector", "North", "Central")) |
yes |
|
Juvenile: Inshore vs Offshore
|
Subtract column with all inshore vs offshore comparison from sediment: inshore vs offshore
[3] "Shore_Offshore_vs_Inshore" 1,0,1,0,0,0,0,0 [5] "TypeSediment.ShoreOffshore" 1,0,0,0,1,0,0,0 Subtract the two to get the following: results(varstab2, contrast=c(0,0,1,0,-1,0,0,0))
|
Not sure |
|
Sediment: Inshore vs Offshore
|
results(varstab2, name="TypeSediment.ShoreOffshore") |
yes |
|
Juvenile: North vs Central |
Subtract column with all north vs central comparison from sediment: north vs central
[4] "Sector_North_vs_Central" 1,0,0,1,0,0,0,0 [6]"TypeSediment.SectorNorth" 1,0,0,0,0,1,0,0 Subtract the two to get the following: results(varstab2, contrast=c(0,0,0,1,0,-1,0,0))
|
Not sure |
|
Sediment: North vs Central |
results(varstab2, name=”TypeSediment.SectorNorth”) |
yes |
|
Juvenile: Inshore North vs. Inshore Central |
[1] "Intercept" 1,0,0,0,0,0,0,0 which is juvenile, inshore, central [4] "Sector_North_vs_Central" 0,0,0,1,0,0,0,0
Subtract the two to get the following: results(varstab2, contrast=c(1,0,0,-1,0,0,0,0))
|
Not sure |
|
Juvenile: Inshore North vs. Offshore Central |
|
? |
|
Juvenile: Inshore North vs. Offshore North |
|
? |
|
Juvenile: Offshore Central vs. Inshore Central |
|
? |
|
Juvenile: Offshore North vs. offshore Central |
[7] "ShoreOffshore.SectorNorth" results(varstab2, name="ShoreOffshore.SectorNorth") |
Ok |
|
Juvenile: Offshore North vs. Inshore Central |
|
? |
|
Sediments: Inshore North vs. Inshore Central |
|
? |
|
Sediments: Inshore North vs. Offshore Central |
|
? |
|
Sediments: Inshore North vs. Offshore North |
|
? |
|
Sediments: Offshore Central vs. Inshore Central |
|
? |
|
Sediments: Offshore North vs. offshore Central |
[8]"TypeSediment.ShoreOffshore.SectorNorth" results(varstab2, name="TypeSediment.ShoreOffshore.SectorNorth") |
Ok |
|
Sediments: Offshore North vs. Inshore Central |
|
? |
I feel confident in the comparisons with the a "Yes" in the last column. However, could someone please check that I have specified the contrasts in the rows with "Not sure?" Especially the contrasts "Juvenile: Inshore North vs. Inshore Central." If I can get a confirmation that that contrast is set correctly, then I think I can manage the rest.
Thank you for any insights!
cheers,
k
Thank you Michael. I can keep looking through the glmm literature now that you confirm that contrasts in DESeq2 and other negative binomial models are specified the same. I just wasn't sure if there was some unique DESeq2 functionality I wasn't finding that would draw out these contrasts.
cheers,
Hello Michael:
I would like to propose a possible solution to such contrast questions and run it by you.
To estimate higher order contrasts:
1) Construct a design matrix in GLM coding. An example is given here:
https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_glm_sect030.htm
2) Each categorical column (“A1”, “AB21”, etc) in the design can be seen as a separate binary factor. I can also pretend that every column is a quantitative factor (please let me know if it makes a difference for DESeq2). Either way, an additive model that includes all of such factors is fed to DESeq2. Since interactions are actually present, shrinkage of regression coefficients is disabled (betaPrior = FALSE).
3) Obtain the contrast vector as a difference between the two corresponding lsmeans vectors. E.g. for “(A1, B1) vs (A3, B1)” contrast the vector is LSM(AB11) – LSM(AB31). How to obtain lsmeans is described here:
https://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#statug_glm_a0000000871.htm
4) Feed the contrast vector to results() function from DESeq2. Here I hope to rely on the option: “a numeric contrast vector with one element for each element in resultsNames(object) (most general case)”.
Please tell me if DESeq2 can digest an overparameterized design “as is” and whether the results() will work the way I expect.
You can pass in any full rank design matrix to the 'full' argument, and use numeric contrasts to compare coefficients using results(). The current release of DESeq2 has betaPrior=FALSE by default (see NEWS).
If the user can only feed full rank designs to DESeq2, do you know if there is an R package that could help figure out how to specify contrast coefficients for non-trivial cases? If not, it's very hard to use the design produced by model.matrix() if the goal is to estimate non-trivial contrasts.
The
contrastpackage works with 'lm' like objects. So it doesn't work with DESeq2 but you could extract the linear contrast from the output of the contrast() function from that package, for example running on fake data for 'y' but using the same model.matrix. See the example of the contrast package here.For non-trivial contrasts I find it just as hard to specify the desired contrast using that package, as working it out myself.