Standard error and effect size from Limma
2
1
Entering edit mode
Vani ▴ 20
@vani-8145
Last seen 5.7 years ago
United States

Is it possible to find the standard error and effect size of a dataset using limma's lmfit and toptable? Please advise.

effect size Limma toptable Standard Error • 4.4k views
5
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia

The effect sizes are contained in fit$coefficients. The standard errors can be obtained from SE <- sqrt(fit$s2.post) * fit$stdev.unscaled This gives a matrix containing standard errors for every coefficient and every gene. ADD COMMENT 0 Entering edit mode I am confused; shouldn't that be sqrt(fit$s2.post)? Similar to the line 255 of the toptable.R file in limma source package, the confidence intervals are calculated using

sqrt(eb$s2.post[top])*fit$stdev.unscaled[top,coef]*qt(alpha,df=eb$df.total[top]) ADD REPLY 0 Entering edit mode You are right. Now corrected. ADD REPLY 0 Entering edit mode I'm sorry, just tried this code, but in my fit object I don't have 's2.post'. How do I get this? ADD REPLY 0 Entering edit mode If you've run eBayes() on the fit object, then you will haves s2.post. ADD REPLY 0 Entering edit mode However, I'm following this protocol here: https://molepi.github.io/DNAmArray_workflow/06_EWAS.html#correct_for_bias_and_inflation to run a EWAS. So, we start with limma(), but correct for inflation using bacon(). And thus, this is my code. designp <- model.matrix(metadata(Mvalues)$formula, data = colData(Mvalues))
datap <- assays(Mvalues)$data fitp <- limma::lmFit(datap, designp) Given that tstatp <-fitp$coef/fitp$stdev.unscaled/fitp$sigma had gotten me the T-statistic, I assumed fitp$stdev.unscaled/fitp$sigma was equal to standard error. And I assumed that fitp$coef would give me the effect sizes. If I don't run eBayes() where do I get the standard error from? Or is that effectively not possible? Thanks! ADD REPLY 0 Entering edit mode If you want to ask questions about a non-standard workflow, you should start a new post with the appropriate tags. (In this case, it doesn't seem to be a Bioconductor package, so you might as well ask the authors directly.) The workflow in question is a bit bemusing as limma is run without EB shrinkage, which defeats the purpose - you might as well use lm.fit. Anyway, fitp$stdev.unscaled*fitp\$sigma is the standard error of the coefficient.

0
Entering edit mode

Ah, thank you for the answer, this was not immediately clear to me.

4
Entering edit mode
@steve-lianoglou-2771
Last seen 1 day ago
Denali

Yes.

Calling topTable(fit, ..., confint=TRUE) provides you with the effect size (logFC) and 95% confidence intervals (CI.L, CI.R), from which you can back calculate the standard error, if need be.

0
Entering edit mode

Cool thanks. Just a quick question: How would I calculate the standard error from the confidence intervals?

0
Entering edit mode

I'd say, (CI.R - CI.L)/3.92 (since the 95% confidence intervals should be plus/minus 1.96se)