Question: random effects + interaction term + limma
0
5 months ago by
flippy230
flippy230 wrote:

Hi,

I realize that I made a similar post, but didn't realize I was editing that, and didn't know how to go back, so this is a similar post, but different question about my interpretation, rather than the code

I am following section 9.5.2 - 9.5.4 of the limma user guide. I understand that each of the output of each of these three sections should yield identical results, but my results differ in p-value and adjusted p-value, so I suspect that there are issues with my code but can't seem to find out what it is. This is a pre/post design where every patient has received a treatment. Certain patients are responders, others are not. We're interested in this specific interaction effect: the effect of treatment in responders versus non-responders

make sure treatment is a factor variables with proper reference

>phenotype$treatment <- as.factor(phenotype$treatment)
>phenotype$treatment <- relevel(phenotype$treatment, ref = "0")


make sure responder is a factor variables with proper reference

>phenotype$responder <- as.factor(phenotype$responder)
>phenotype$responder <- relevel(phenotype$responder, ref = "non")


convert to cpm --> between sample, controlling for sequencing depth

>countsPerMillion <- cpm(dgList) ##divide each count per million for comparison


we retain only those genes that are represented at least 1cpm reads in at least two samples (cpm=counts per million)

>summary(countsPerMillion)
>countCheck <- countsPerMillion > 1
>keep <- which(rowSums(countCheck) >= 2)
>dgList <- dgList[keep,]
>summary(cpm(dgList))


TMM normalization controls for library size

>dge <- calcNormFactors(dgList, method="TMM")


design matrix

>phenotype$treat <- factor(paste(phenotype$treatment,phenotype$responder,sep=".")) >treat <- paste(phenotype$responder,phenotype$treatment,sep=".") >treat <- factor(treat, levels=c("non.0","non.1","responder.0","responder.1")) >design <-model.matrix(~0 + treat + cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + cov7,data=phenotype) >colnames(design)[1:4] <- levels(treat)  estimate within subject correlation >counts <- voom(dge,design) >corfit <- duplicateCorrelation(counts,design,block=phenotype$study_id)
>corfit$consensus >counts_voom <- voom(dge,design,block=phenotype$study_id,correlation=corfit$consensus)  account for correlation in a model with random effect >fit <- lmFit(counts_voom,design,block=phenotype$study_id,correlation=corfit$consensus) >fit <- eBayes(fit)  make contrasts >cm <- makeContrasts( DOD = (responder.1 - responder.0) - (non.1 - non.0), RTE = responder.1 - responder.0, NRTE = non.1 - non.0, levels=design )  fit model >fit2 <- contrasts.fit(fit,cm) >fit2 <- eBayes(fit2) >results <- topTable(fit2,coef="DOD")  when I follow section 9.5.3: create design matrix >design <-model.matrix(~responder + responder:treatment + cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + cov7,data=phenotype)  estimate within subject correlation >counts <- voom(dge,design) >corfit <- duplicateCorrelation(counts,design,block=phenotype$study_id)
>corfit$consensus >counts_voom <- voom(dge,design,block=phenotype$study_id,correlation=corfit$consensus)  account for correlation in a model with random effect >fit <- lmFit(counts_voom,design,block=phenotype$study_id,correlation=corfit$consensus) >fit <- eBayes(fit)  make contrasts corresponding to the respondernon:treatment1 and responderresponder:treatment1 to extract interaction of the effect of treatment in responders versus nonresponders >fit2 <- contrasts.fit(fit,c(0,0,0,0,0,0,0,0,0,-1,1) >fit2 <- eBayes(fit2) >results <- topTable(fit2)  The output from colnames(fit$design) is:

> [1] "Intercept"  "responderresponder"
> [3]  "cov1" "cov2"
> [5]  "cov3" "cov4"
> [7]  "cov5" "cov6"
> [9]  "cov7" "respondernon:treatment1"
> [11] "responderresponder:treatment1"


Based on this, and info gathered from my basic LDA class, I interpreted the coefficients as such: Intercept - The treatment0 and cov1 thru cov7 at 0 effect in nonresponder responderresponder - the difference between the non-responders and responders at treatment0 and cov1 thru cov7 at 0 respondernon:treatment1 - the difference between nonresponders at treatment 0 and treatment 1 adjusted for cov1-cov7 responderresponder:treatment1 - the difference between responders at treatment 0 and treatment 1 adjusted for cov1-cov7

In this case, I would assume that the difference between responderresponder:treatment1 and respondernon:treatment1 would be the interaction of interest. However, I don't believe that this is the case, given that the results I run with the 9.5.2 results don't match these (a la 9.5.3) results.

modified 5 months ago • written 5 months ago by flippy230

Each time I have answered a question from you on this forum, you subsequently deleted the question so that my answer and my time was wasted. So I will keep away from your questions from now on.

Anyway, you might find it helpful to read the note on the contrasts.fit help page, which possibly answers your question in this case.