Question

Standard error and effect size from Limma

1

Entering edit mode

Vani ▴ 20

@vani-8145

Last seen 8.9 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.5k views

ADD COMMENT • link updated 9.3 years ago by Gordon Smyth 51k • written 9.3 years ago by Vani ▴ 20

score 5 · Answer 1 · 2015-07-20

5

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 10 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 • link 9.3 years ago • updated 9.2 years ago Gordon Smyth 51k

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 • link 9.2 years ago January Weiner ▴ 240

0

Entering edit mode

You are right. Now corrected.

ADD REPLY • link 9.2 years ago Gordon Smyth 51k

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 • link 7.3 years ago s.w.vanderlaan ▴ 30

0

Entering edit mode

If you've run eBayes() on the fit object, then you will haves s2.post.

ADD REPLY • link 7.3 years ago Gordon Smyth 51k

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 • link 7.3 years ago s.w.vanderlaan ▴ 30

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 • link 7.3 years ago Aaron Lun ★ 28k

0

Entering edit mode

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

ADD REPLY • link 7.3 years ago s.w.vanderlaan ▴ 30

score 4 · Answer 2 · 2015-07-20

4

Entering edit mode

Steve Lianoglou ★ 13k

@steve-lianoglou-2771

Last seen 21 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 • link 9.3 years ago Steve Lianoglou ★ 13k

0

Entering edit mode

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

ADD REPLY • link 9.3 years ago Vani ▴ 20

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 • link 8.7 years ago Arnaud Amzallag ▴ 100