Entering edit mode
3.7 years ago
sindrelee
▴
10
Cross-posted: https://support.bioconductor.org/p/133284/
Please see the example below.
set.seed(1)
pheno=cbind.data.frame(IDs=as.factor(c(1:20)),
gr=as.factor(c(rep("A",10),rep("B",10))),
oc=rnorm(20)
)
dat=rbind.data.frame(varA=rnorm(20),
varB=rnorm(20),
varC=rnorm(20),
varD=rnorm(20)
)
colnames(dat)=IDs
rownames(pheno)=IDs
form=as.formula(paste0("~oc*gr"))
design=model.matrix(form,data=pheno)
colnames(design)=make.names(colnames(design))
colnames(design)
[1] "X.Intercept." "oc" "grB" "oc.grB"
library(limma)
fit=lmFit(dat,design)
fit2=eBayes(fit)
topTable(fit2,coef=2) # Pooled oc results across groups
topTable(fit2,coef=3) # Group A vs B
topTable(fit2,coef=4) # Interaction results
The question is then: How to I get within group A and within group B results to further examine the nature of the interaction? Eg. 'Post-Hoc' testing: ~oc+grA and ~oc+grB. I know I can re-run the analyses stratifying by group, but it seems unnecessary.
fit3 <- contrasts.fit(fit, c(1, 0, 0, -1)) # Not sure how to specify this correctly
fit3 <- eBayes(fit3)
topTable(fit3)