Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:20:14.420Z
Answers.Text: Aren’t alpha thresholds usually exclusive? I.e. it should be `< 0.05`, not `<= 0.05`.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:23:31.717Z
Answers.Text: 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).
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T11:10:50.080Z
Answers.Text: @KonradRudolph Yes, good catch, I've fixed that. I had >= and min() in my mind apparently.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-12T18:22:42.213Z
Answers.LastActivityDate: 2017-06-13T11:14:15.507Z
Answers.LastEditDate: 2017-06-13T11:14:15.507Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:18:22.030Z
Answers.Text: Why the `dummy <-` assignment?
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:21:19.270Z
Answers.Text: 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).
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:22:06.380Z
Answers.Text: Ah. The writeout never bothered me here. It’s not like it would be written to the output in an actual script/document.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:27:11.697Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-12T20:51:34.303Z
Answers.LastActivityDate: 2017-06-13T10:21:37.260Z
Answers.LastEditDate: 2017-06-13T10:21:37.260Z
Answers.UserDisplayName: gringer
Body: <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>
  
Tags: <r>
Title: How to correct alpha, and not p-values themselves, for visualization purposes
ThreadText: # 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**: @username_2 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());
  
  ```
  
  [![Volcano plot of p-values with FDR greyzone](https://i.sstatic.net/l8Ldk.png)](https://i.sstatic.net/l8Ldk.png)
  
  > **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.  
  
  ---
  
  