Standard error and effect size from Limma
2
1
Entering edit mode
Vani ▴ 20
@vani-8145
Last seen 9.0 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 • 8.6k views
ADD COMMENT •
5
Entering edit mode
@gordon-smyth
Last seen 13 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.

ADD REPLY •
0
Entering edit mode

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

 

ADD REPLY •
4
Entering edit mode
@steve-lianoglou-2771
Last seen 22 months ago
United States

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.
ADD COMMENT •
0
Entering edit mode

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

ADD REPLY •
0
Entering edit mode

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

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 565 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6