diff --git a/lsi.Rmd b/lsi.Rmd index f7f4c9e..be75739 100644 --- a/lsi.Rmd +++ b/lsi.Rmd @@ -14,9 +14,8 @@ output: ## Brain imaging study -```{r fig.align="center", out.width = '80%', echo=FALSE} +```{r, fig.align="center", out.width = '80%', echo=FALSE} knitr::include_graphics("./figures/DTIColor.jpg") -library("magick") ``` - Diffusion Tensor Imaging (DTI) data @@ -42,6 +41,8 @@ The dataset `dti` contains library(tidyverse) library(locfdr) library(gganimate) +library(magick) +library(gridExtra) dti <- read_csv("https://raw.githubusercontent.com/statOmics/HDA2020/data/dti.csv", col_types = cols()) @@ -64,16 +65,11 @@ pZ <- dti %>% We will now plot the animated graph -```{r eval=FALSE} -pZ -``` - - __WARNING__: The animated graph will only be visible in the HTML output, not in PDF format. If you're reading the PDF version, check [online](https://statomics.github.io/HDDA21/lsi.html#111_Data_Exploration) for the animated graph. -```{r echo=FALSE, message=FALSE, eval=knitr::is_html_output()} +```{r, message=FALSE, eval=knitr::is_html_output()} animate(pZ, nframes = 103, end_pause = 3) ``` @@ -119,11 +115,7 @@ pPval <- dti %>% We will now plot the animated graph -```{r eval=FALSE} -pPval -``` - -```{r echo=FALSE, message=FALSE, eval=knitr::is_html_output()} +```{r, message=FALSE, eval=knitr::is_html_output()} animate(pPval, nframes = 103, end_pause = 3) ``` @@ -186,8 +178,7 @@ Hence, for any $0<\alpha<1$, The distribution of the z-statistic and the p-values under $H_0$ are illustrated below: ```{r} -library(gridExtra) - +set.seed(123) simData <- tibble( z.value = rnorm(20000) ) @@ -197,8 +188,8 @@ simData <- simData %>% mutate(p.value = 2*(1-pnorm(abs(z.value)))) p1 <- simData %>% ggplot(aes(x = z.value)) + geom_histogram( - aes(y=..density..), - color = "black") + + aes(y = after_stat(density)), + color = "black", bins = 50) + stat_function(fun = dnorm, args=list(mean=0, sd=1)) p2 <- simData %>% @@ -223,8 +214,8 @@ $\rightarrow$ Probability to have a false positive (FP) among all m simultatenou test $>>> \alpha= 0.05$ - Indeed for each non differential voxel we have a probability of 5% to return a FP. -- In a typical experiment the majority of the voxel are non differential. -- So an upperbound of the expected FP is $m \times \alpha$ or $15443 \times 0.05=`r round(15443*0.05,0)`$. +- In a typical experiment the majority of the voxel are non differential. +- So an upperbound of the expected FP is $m \times \alpha$ or $15443 \times 0.05=`r round(15443*0.05,0)`$. $\rightarrow$ Hence, we are bound to call many false positive voxels each time we run the experiment. @@ -290,15 +281,15 @@ FWER$=0.05$ and $m=15443$: $\alpha=0.0000033$. We will argue that this procedure is too stringent for large $m$. -### Bonferroni method +### Bonferroni method -The Bonferroni method is another method that is widely used to control the FWER: +The Bonferroni method is another method that is widely used to control the FWER: -- assess each test at +- assess each test at \[\alpha_\text{adj}=\frac{\alpha}{m}\] - The method does not assume independence of the test statistics. -- Again, the method is very conservative! +- Again, the method is very conservative! --- @@ -322,7 +313,7 @@ First we give a very general definition of an **adjusted $p$-value**. The corrected $p$-value should be reported. It accounts for the multiplicity problem and it can be compared directly to the nominal FWER level to make calls at the FWER level. -- adjusted p-values for Bonferroni method: +- adjusted p-values for Bonferroni method: \[p_\text{adj}=\text{min}\left(p \times m,1\right)\] --- @@ -331,7 +322,7 @@ First we give a very general definition of an **adjusted $p$-value**. ## Introduction -In large scale inference it would be more interesting to tolerate a few false positives as long as they do not dominate the toplist +In large scale inference it would be more interesting to tolerate a few false positives as long as they do not dominate the toplist We first introduce some notation: @@ -356,8 +347,8 @@ The table shows the results of $m$ hypothesis tests in a single experiment. --- -- Note that the table is not completely observable. -- Indeed, we can only observe the bottom row! +- Note that the table is not completely observable. +- Indeed, we can only observe the bottom row! - The table is introduced to better understand the concept of FWER and to introduce the concept of the false discovery rate (FDR). --- @@ -372,18 +363,17 @@ The FWER can now be reexpressed as \[ \text{FWER}=\text{P}\left[\text{reject at least one }H_{0i} \mid H_0\right] = \text{P}\left[FP>0\right] . \] - -- However, we know that the FWER is very conservative in large scale inference problems. -- Therefore it would be more interesting to tolerate a few false positives as long as they do not dominate the toplist +- However, we know that the FWER is very conservative in large scale inference problems. +- Therefore it would be more interesting to tolerate a few false positives as long as they do not dominate the toplist -The **False Discovery Proportion (FDP)** is the fraction of false positives that are returned, i.e. +The **False Discovery Proportion (FDP)** is the fraction of false positives that are returned, i.e. \[ FDP = \frac{FP}{R} \] -- However, this quantity cannot be observed because in practice we only know the number of voxels for which we rejected $H_0$, $R$. +- However, this quantity cannot be observed because in practice we only know the number of voxels for which we rejected $H_0$, $R$. - But, we do not know the number of false positives, $FP$. @@ -391,7 +381,7 @@ Therefore, Benjamini and Hochberg, 1995, defined The **False Discovery Rate (FDR \[ \text{FDR} = \text{E}\left[\frac{FP}{R}\right] =\text{E}\left[\text{FDP}\right] \] -the expected FDP, in their seminal paper Benjamini, Y. and Hochberg, Y. (1995). "Controlling the false discovery rate: a practical and powerful approach to multiple testing". Journal of the Royal Statistical Society Series B, 57 (1): 289–300. +the expected FDP, in their seminal paper Benjamini, Y. and Hochberg, Y. (1995). "Controlling the false discovery rate: a practical and powerful approach to multiple testing". Journal of the Royal Statistical Society Series B, 57 (1): 289–300. - An FDR of 1% means that on average we expect 1% false positive voxels in the list of voxels that are called significant. @@ -401,10 +391,10 @@ the expected FDP, in their seminal paper Benjamini, Y. and Hochberg, Y. (1995). Consider $m = 1000$ voxels -- Suppose that a researcher rejects all null hypotheses for which $p < 0.01$. +- Suppose that a researcher rejects all null hypotheses for which $p < 0.01$. -- If we use $p < 0.01$, we expect $0.01 \times m_0$ tests to return false positives. -- A conservative estimate of the number of false positives that we can expect can be obtained by considering that the null hypotheses are true for all features, $m_0 = m = 1000$. +- If we use $p < 0.01$, we expect $0.01 \times m_0$ tests to return false positives. +- A conservative estimate of the number of false positives that we can expect can be obtained by considering that the null hypotheses are true for all features, $m_0 = m = 1000$. - We then would expect $0.01 \times 1000 = 10$ false positives ($FP=10$). - Suppose that the researcher found 200 voxels with $p<0.01$ ($R=200$). @@ -419,7 +409,7 @@ Consider $m = 1000$ voxels 1. Let $p_{(1)}\leq \ldots \leq p_{(m)}$ denote the ordered $p$-values. -2. Find the largest integer $k$ so that +2. Find the largest integer $k$ so that $$ \frac{p_{(k)} \times m}{k} \leq \alpha $$ @@ -485,11 +475,11 @@ At the 5% FDR, `r sum(dti$padj < 0.05)` voxels are returned as significantly dif - Bonferroni: $\alpha_\text{adj}=`r format(0.05/nrow(dti),digits=2)` \rightarrow$ `r sum(dti$p.value<(0.05/nrow(dti)))` voxels are significant at the Bonferroni FWER -- BH-FDR: +- BH-FDR: 1. ordered $p$-values. -2. Find the largest integer $k$ so that +2. Find the largest integer $k$ so that $$ \frac{p_{(k)} \times m}{k} \leq \alpha $$ @@ -504,15 +494,15 @@ Otherwise none of the null hypotheses is rejected. ```{r echo=FALSE} alpha <- 0.05 -res <- dti %>% +res <- dti %>% select("z.value","p.value","padj") %>% arrange(p.value) res$padjNonMonoForm <- paste0(nrow(res)," x pval /",1:nrow(res)) res$padjNonMono <- res$p.value *nrow(res) /(1:nrow(res)) res$adjAlphaForm <- paste0(1:nrow(res)," x ",alpha,"/",nrow(res)) -res$adjAlpha <- alpha * (1:nrow(res))/nrow(res) -res$"pval < adjAlpha" <- res$p.value < res$adjAlpha -res$"padj < alpha" <- res$padj < alpha +res$adjAlpha <- alpha * (1:nrow(res))/nrow(res) +res$"pval < adjAlpha" <- res$p.value < res$adjAlpha +res$"padj < alpha" <- res$padj < alpha res[1:10,] %>% knitr::kable() res[11:20,] %>% knitr::kable() res[21:30,] %>% knitr::kable() @@ -564,7 +554,7 @@ animate(pFDR, nframes = 103, end_pause = 3) \text{FDR}=\text{E}\left[TP/R\right]=\text{E}\left[\text{FDP}\right] \leq \frac{m_0}{m} \alpha \leq \alpha . \] ---- +--- ### Extension @@ -579,7 +569,7 @@ The inequality \] shows that BH1995 is a conservative method, i.e. it controls the FDR at the safe side, i.e. when one is prepared to control the FDR at the nominal level $\alpha$, the BH95 will guarantee that the true FDR is not larger than the nominal level (when the assumptions hold). -- More interestingly is that $\frac{m_0}{m} \alpha$ is in between the true FDR and the nominal FDR. +- More interestingly is that $\frac{m_0}{m} \alpha$ is in between the true FDR and the nominal FDR. - Suppose that $m_0$ were known and that the BH95 method were applied at the nominal FDR level of $\alpha=m/m_0 \alpha^*$, in which $\alpha^*$ is the FDR level we want to control. Then the inequality gives \[ @@ -763,8 +753,8 @@ zSim %>% as_tibble %>% ggplot(aes(x = zSim)) + geom_histogram( - aes(y=..density..), - color = "black") + + aes(y = after_stat(density)), + color = "black", bins = 50) + stat_function(fun = dnorm, args = list( mean = 0,