Hi Michael,

I would like to make sure that I am interpretating the p-values correctly when using the likelihood ratio test with a multifactor design.

We are interested in identifying significant genes that are differentially expressed between 3 disease states while controlling for sex. I understand that the likelihood ratio test is testing whether the reduced model (dispersions estimated on sex alone) or the full model (dispersions estimated on sex and condition) better fits the dataset.

So, if an adjusted p-value > 0.05 it indicates that those genes are better explained by the reduced model and are either affected by sex only or not at all. However if an adjusted p value < 0.05 , does it indicate the those genes are better explained / unique to disease state while controlling for sex OR better explained by disease state and sex combined?

(coldata = data.frame(row.names=colnames(countdata), disease, sex))
dds = DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~sex+disease)
dds

# Run the DESeq pipeline
ddsLRT = DESeq(dds, test="LRT", full=~sex+disease, reduced=~sex)
res=results(ddsLRT)
resSig=subset(res, padj<0.05)

"I understand that the likelihood ratio test is testing whether the reduced model (dispersions estimated on sex alone) or the full model (dispersions estimated on sex and condition) better fits the dataset."

Close, but not exactly. It tests whether the increase in the log likelihood from the additional coefficients would be expected if those coefficients were equal to zero. It doesn't mean the reduced model is a good model or a good fit.

The way to interpret the p-value is: if the adjusted p-value is small, then for the set of genes with those small adjusted p-values, the additional coefficient in full and not in reduced increased the log likelihood more than would be expected if their true value was zero.

See here for some more background:

https://en.wikipedia.org/wiki/Likelihood-ratio_test