Is it possible to find the standard error and effect size of a dataset using limma's lmfit and toptable? Please advise.
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.
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
You are right. Now corrected.
I'm sorry, just tried this code, but in my `fit` object I don't have 's2.post'. How do I get this?
If you've run eBayes() on the fit object, then you will haves s2.post.
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?
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.
Ah, thank you for the answer, this was not immediately clear to me.
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.
topTable(fit, ..., confint=TRUE)
Cool thanks. Just a quick question: How would I calculate the standard error from the confidence intervals?
I'd say, `(CI.R - CI.L)/3.92` (since the 95% confidence intervals should be plus/minus 1.96se)
Login before adding your answer.
Use of this site constitutes acceptance of our User Agreement and Privacy