<p>The only way to get the alpha levels is to determine what they will be with <code>p.adjust()</code>, since they will depend on the distribution of your unadjusted p values. The general steps you should follow will be:</p>

<ol>
<li>Add a column of adjusted p-values to your dataframe (<code>mydata$padj = p.adjust(mydata, method="BH")</code>, which is the same as FDR and saves a character).</li>
<li>Use <code>which</code> and <code>max</code> to determine your two alpha threshold (e.g., <code>max(mydata$pvalue[mydata$padj &lt; 0.05])</code></li>
</ol>

<p>Then you can adjust your plots however you like (presumably with some horizontal lines at the various alphas). Whether you take the smallest non-significant value or the largest significant value is up to you, just describe what "dots on the line" represent.</p>

CC BY-SA 3.0
2017-06-13T10:20:14.420Z
Aren’t alpha thresholds usually exclusive? I.e. it should be `< 0.05`, not `<= 0.05`.
Konrad Rudolph
CC BY-SA 3.0
2017-06-13T10:23:31.717Z
If it matters, you're not using p-values correctly. See the [ASA position statement](http://www.amstat.org/asa/files/pdfs/P-ValueStatement.pdf).
gringer
CC BY-SA 3.0
2017-06-13T11:10:50.080Z
@KonradRudolph Yes, good catch, I've fixed that. I had >= and min() in my mind apparently.
Devon Ryan
CC BY-SA 3.0
2017-06-12T18:22:42.213Z
2017-06-13T11:14:15.507Z
2017-06-13T11:14:15.507Z
Devon Ryan
<p>When you're doing a p-value adjustment, the same unadjusted p-value in different genes can be given different adjusted p-values depending on other factors. That means that you can't directly draw a line associated with the FDR on a plot of unadjusted p-value.</p>

<p>One possibility would be to take a range of values that are close to the FDR threshold (e.g. the 20 values closest to threshold), and draw a p-value greyzone within that region:</p>

<pre><code>#!/usr/bin/Rscript
values &lt;- c(rnorm(10000),rnorm(100, mean=1.5));
val.mean &lt;- median(values);
val.diffs &lt;- abs(values - median(values));
val.reldiffs &lt;- (values - median(values));

val.pval &lt;- pnorm(val.diffs, mean = mean(val.diffs),
                   sd=sd(val.diffs), lower.tail=FALSE);
val.padj &lt;- p.adjust(val.pvals, method="BH");
fdr.threshold &lt;- 0.1;
close.bh &lt;- order(abs(val.padj - fdr.threshold))[1:20];

png("SE.663.png");
plot(val.reldiffs, -log10(val.pval),
     col=ifelse(1:10100 &lt;= 10000,"darkblue","darkgreen"));
abline(h=-log10(0.05), col="red");
text(0,-log10(0.05),"p=0.05", pos=1);
abline(h=range(-log10(val.pval[close.bh])), col="#00000040", lty="dashed");
rect(xleft=min(val.reldiffs)*2, xright=max(val.reldiffs)*2,
     ytop=max(-log10(val.pval[close.bh])),
     ybottom=min(-log10(val.pval[close.bh])), col="#00000020", border=NA);
text(0,min(-log10(val.pval[close.bh])),"FDR=0.1", pos=1);
invisible(dev.off());
</code></pre>

<p><a href="https://i.sstatic.net/l8Ldk.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/l8Ldk.png" alt="Volcano plot of p-values with FDR greyzone"></a></p>

CC BY-SA 3.0
2017-06-13T10:18:22.030Z
Why the `dummy <-` assignment?
Konrad Rudolph
CC BY-SA 3.0
2017-06-13T10:21:19.270Z
It's my default workaround to stop reporting `null device / 1` to standard out. I've changed it to `invisible(dev.off())`, which does the same thing (and is probably more descriptive about its purpose).
gringer
CC BY-SA 3.0
2017-06-13T10:22:06.380Z
Ah. The writeout never bothered me here. It’s not like it would be written to the output in an actual script/document.
Konrad Rudolph
CC BY-SA 3.0
2017-06-13T10:27:11.697Z
Sometimes I send stdout to files and parse it. I've just got into the habit of getting rid of things like this so that I don't have weird errors cropping up for unknown reasons. I also have an issue with the `scan` function, which I silence when I use it.
gringer
CC BY-SA 3.0
2017-06-12T20:51:34.303Z
2017-06-13T10:21:37.260Z
2017-06-13T10:21:37.260Z
gringer
<p>I have a set of differentially methylated/expressed/whatever entities with p-values attached (example below).</p>

<pre><code>entity_name    p-value    magnitude
entity1        0.04459    0.68
entity2        0.02283    0.99
...
entity_n       0.78       0.025
</code></pre>

<p>Typically, I apply the p.adjust function in R with the "fdr" (Benjamini-Hochberg) approach to leave me with p-values adjusted to control the FDR.</p>

<pre><code>adjusted &lt;- p.adjust(mydata,"fdr")
</code></pre>

<p>However, I am interested in showing a volcano plot with the unadjusted p-values, and two alpha levels: 0.05 and one that corresponds to the correction. What is the best way to get this alpha? Is it appropriate to set the "corrected alpha" to the lowest original p-value that doesn't pass FDR correction?</p>

<r>
How to correct alpha, and not p-values themselves, for visualization purposes
# How to correct alpha, and not p-values themselves, for visualization purposes

(Asked by: username_0 on 2017-06-12)

I have a set of differentially methylated/expressed/whatever entities with p-values attached (example below).
entity_name    p-value    magnitude
entity1        0.04459    0.68
entity2        0.02283    0.99
...
entity_n       0.78       0.025

Typically, I apply the p.adjust function in R with the "fdr" (Benjamini-Hochberg) approach to leave me with p-values adjusted to control the FDR.
adjusted <- p.adjust(mydata,"fdr")

However, I am interested in showing a volcano plot with the unadjusted p-values, and two alpha levels: 0.05 and one that corresponds to the correction. What is the best way to get this alpha? Is it appropriate to set the "corrected alpha" to the lowest original p-value that doesn't pass FDR correction?

**Tags:** r

---

## 2 answers

### Answer by username_1 on 2017-06-12. Score: 7 | Accepted answer

The only way to get the alpha levels is to determine what they will be with `p.adjust()`, since they will depend on the distribution of your unadjusted p values. The general steps you should follow will be:

1. Add a column of adjusted p-values to your dataframe (`mydata$padj = p.adjust(mydata, method="BH")`, which is the same as FDR and saves a character).
2. Use `which` and `max` to determine your two alpha threshold (e.g., `max(mydata$pvalue[mydata$padj < 0.05])`

Then you can adjust your plots however you like (presumably with some horizontal lines at the various alphas). Whether you take the smallest non-significant value or the largest significant value is up to you, just describe what "dots on the line" represent.

> **Comments:**
> **username_2**: Aren’t alpha thresholds usually exclusive? I.e. it should be `< 0.05`, not `<= 0.05`.  
> **username_3**: If it matters, you're not using p-values correctly. See the [ASA position statement](http://www.amstat.org/asa/files/pdfs/P-ValueStatement.pdf).  
> **username_1**: @KonradRudolph Yes, good catch, I've fixed that. I had >= and min() in my mind apparently.  

---

### Answer by username_3 on 2017-06-12. Score: 3

When you're doing a p-value adjustment, the same unadjusted p-value in different genes can be given different adjusted p-values depending on other factors. That means that you can't directly draw a line associated with the FDR on a plot of unadjusted p-value.
One possibility would be to take a range of values that are close to the FDR threshold (e.g. the 20 values closest to threshold), and draw a p-value greyzone within that region:
#!/usr/bin/Rscript
values <- c(rnorm(10000),rnorm(100, mean=1.5));
val.mean <- median(values);
val.diffs <- abs(values - median(values));
val.reldiffs <- (values - median(values));

val.pval <- pnorm(val.diffs, mean = mean(val.diffs),
                   sd=sd(val.diffs), lower.tail=FALSE);
val.padj <- p.adjust(val.pvals, method="BH");
fdr.threshold <- 0.1;
close.bh <- order(abs(val.padj - fdr.threshold))[1:20];

png("SE.663.png");
plot(val.reldiffs, -log10(val.pval),
     col=ifelse(1:10100 <= 10000,"darkblue","darkgreen"));
abline(h=-log10(0.05), col="red");
text(0,-log10(0.05),"p=0.05", pos=1);
abline(h=range(-log10(val.pval[close.bh])), col="#00000040", lty="dashed");
rect(xleft=min(val.reldiffs)*2, xright=max(val.reldiffs)*2,
     ytop=max(-log10(val.pval[close.bh])),
     ybottom=min(-log10(val.pval[close.bh])), col="#00000020", border=NA);
text(0,min(-log10(val.pval[close.bh])),"FDR=0.1", pos=1);
invisible(dev.off());

> **Comments:**
> **username_2**: Why the `dummy <-` assignment?  
> **username_3**: It's my default workaround to stop reporting `null device / 1` to standard out. I've changed it to `invisible(dev.off())`, which does the same thing (and is probably more descriptive about its purpose).  
> **username_2**: Ah. The writeout never bothered me here. It’s not like it would be written to the output in an actual script/document.  
> **username_3**: Sometimes I send stdout to files and parse it. I've just got into the habit of getting rid of things like this so that I don't have weird errors cropping up for unknown reasons. I also have an issue with the `scan` function, which I silence when I use it.  

---

