Search
Question: Standard error and effect size from Limma
1
gravatar for Vani
2.6 years ago by
Vani20
United States
Vani20 wrote:

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

ADD COMMENTlink modified 2.6 years ago by Gordon Smyth32k • written 2.6 years ago by Vani20
5
gravatar for Gordon Smyth
2.6 years ago by
Gordon Smyth32k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth32k wrote:

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 COMMENTlink modified 2.5 years ago • written 2.6 years ago by Gordon Smyth32k

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 REPLYlink modified 2.5 years ago • written 2.5 years ago by January Weiner230

You are right. Now corrected.

ADD REPLYlink modified 2.5 years ago • written 2.5 years ago by Gordon Smyth32k

I'm sorry, just tried this code, but in my `fit` object I don't have 's2.post'. How do I get this?

ADD REPLYlink written 6 months ago by s.w.vanderlaan20

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

ADD REPLYlink written 6 months ago by Gordon Smyth32k

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 REPLYlink modified 6 months ago • written 6 months ago by s.w.vanderlaan20

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 REPLYlink modified 6 months ago • written 6 months ago by Aaron Lun18k

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

 

ADD REPLYlink written 6 months ago by s.w.vanderlaan20
4
gravatar for Steve Lianoglou
2.6 years ago by
Genentech
Steve Lianoglou12k wrote:

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 COMMENTlink written 2.6 years ago by Steve Lianoglou12k

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

ADD REPLYlink written 2.6 years ago by Vani20

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

ADD REPLYlink written 23 months ago by Arnaud Amzallag100
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 333 users visited in the last hour