Intuitive statistics

“Pose ton esprit, Michel”
Catherine Pirenne, my grand-mother

Introduction

The aim of this course is:

1. To make you think more deeply about some fundamental concepts in statistics, including:
   • the meaning of p-value,
     
   • the issue of multiple testing,
     
   • the meaning of the false discovery rate (FDR), and
     
   • the winners’ curse
2. To introduce you to basic programming in R, including to compute p-values by simulation and permutation.

I have found this book especially to be an excellent introductions to R:

• “Hands-on programming with R” by Garrett Grolemund (https://rstudio-education.github.io/hopr/)

Other references which may be useful to you include:

• “R for data science” by Hadley Wickham & Garrett Grolemund (https://r4ds.had.co.nz)
  
• “R graphics cookbook” by Winston Chang (http://users.metu.edu.tr/ozancan/R%20Graphics%20Cookbook.pdf)

You can download this R notebook and other material from: https://www.gigauag.uliege.be/cms/c_4766333/en/gigauag-teaching

I am relatively new to R myself. So be clement please :-).

You will need to install the following R packages and libraries to use this Notebook

# packages and libraries to install
###################################
#install.packages("pryr")
#install.packages("ggplot2")
#install.package("RColorBrewer)
library("pryr")
library("ggplot2")
library("RColorBrewer")
brewer.pal(9,"Blues")

[1] "#F7FBFF" "#DEEBF7" "#C6DBEF" "#9ECAE1" "#6BAED6" "#4292C6" "#2171B5" "#08519C"
[9] "#08306B"

You may also want to check your working directory before starting

# Check / set your working directory
####################################
setwd("/Users/michelgeorges/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/R_INT_STAT")
getwd()

[1] "/Users/michelgeorges/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/R_INT_STAT"

About p-values, null (H0) and alternative (H1) hypotheses

PLAYING WITH COINS

A common modus operandi in science is to challenge a preconceived model (H0) with observations. If the likelihood of the observations is very small under H0, one may put the validity of H0 in question and possibly “reject” it.

To illustrate this, we will play at “heads or tails” (i.e. “pile ou face”) to test whether coins are balanced, i.e. have an equal chance to fall on the head vs tail face when tossed. Each of us will test a 1-€ coin by tossing it, say 50 times and record the number of times it falls on head. If the coin is balanced (=H0) we expect 25 heads and 25 tails. All different values are a “sign of bias”, i.e. an excess of heads (n_H > 25) or an excess of tails (n_T > 25). Yet we are not worried about the balance of a coin if we observe 26 heads and 24 tails (or vice versa), or even 25 heads and 23 tails (or vice versa). But should we be worried about the coin if we observe f.i. 35 heads and 15 tails?

To address this, we need to be able to compute the probability of the observations under H0 (= balanced coin). In this peculiar case, the probability of a given outcome (n_H and n_T) can easily be computed using the binomial distribution:

[(n_H+nT)!/n_H!n_T!] x 0.5^(n_H+n_T)

We will now use R to generate the corresponding probability distribution. For that we will generate a data frame in which the first column (NR) corresponds to the observed number of “heads”, and six other columns are prepared to receive (i) the corresponding probability (which we will call F), (ii) the deviation from expectation (DEV), (iii) the probability to observe an equally or more deviant series (P for p-value), (iv) a column for a logical (T), (v) …, (vi) … We will assume that - to test its balance - we toss the coin 50 times.

nr <- 50
# Generate data frame
toss <- data.frame(NR=0:nr, F=0, DEV=0, P=0, T=0, BIN=0, SIM=0)
toss

# Generate function computing combinations of n objects taken x at a time
comb = function(n, x) {
  factorial(n) / factorial(n-x) / factorial(x)
}
# Compute probabilities "f"
toss[,2] = 0.5^nr * comb(nr,toss[,1])
toss

# Plot probabilities
plot(toss[,1],toss[,2],col= "blue")

We see that this distribution resembles a normal distribution, “balanced” outcomes are more likely than “unbalanced” ones due to the “nCr” term in the equation.

To test whether our coin is tricked or not we actually utilize a different probability which we will refer to as “the p-value” found in all publications. The question we actually want to ask when we throw 35 heads and 15 tail, is “What was the probability to throw an equally unbalanced or even more unbalanced series?” Well, we could have thrown 50/0, 49/1, 48/2, … 35/15, but also 0/50, 1/49, 2/48, … 15/35, which are all outcomes that are equally or more “unbalanced”. Indeed, when we checked for the balance or not of the coin we didn’t decide whether the bias should favor heads or tails. We performed a so-called “two-tailed test”. Had we suspected in advance that the bias would favor the heads, we would only have considered 50/0, 49/1, 48/2, … 35/15 in a so-called “one-tailed” test. When looking at the probability distribution, and assuming it corresponds to a continuous distribution, our intend is to determine the part of the surface under the curve that lies outside of the boundaries between 34/16 and 35/15 and 16/34 and 15/35 in the two-tailed test, or on the left of the 34/16-35/15 boundary in the “one-tailed” test.

We will now use R to compute this probability for the different values of n_H using, first, our self-made little script and then check whether we obtain the same results as a dedicted R function. Note the use of “logical tests, subsetting and element-wise execution”.

# Compute the probability to have an equally or even more biased toss with our own litte script
################################################################################################
# Computing deviation from expectation
toss$DEV = abs(toss[,1]-nr/2)
toss

# Compute p-values
for (i in 0:nr){
  toss$T = toss$DEV >= toss[i+1,3]
  toss[i+1,4] = sum(toss$F*toss$T)
}
toss

# Plot p-values
plot(toss[,1],toss[,4],col="red", xlab = "NR HEADS", ylab = "P/F")
points(toss[,1],toss[,2],col="blue")

# Compare with R function
#########################
for (i in 0:nr){
  bin <- binom.test(x = toss[i+1,1],n = nr, p = 0.5, alternative = "two.sided")
  toss[i+1,6] = bin$p.value
}
toss

# Plot p-values
###############
lines(toss[,1],toss[,6],col="red")

# Plot -log(p)
##############
plot (toss[,1],-log10(toss[,6]),col="red",xlab = "NR HEADS", ylab = "-log(P)")
abline(h = -log10(0.05))

For 35 heads and 15 tails, the p-value is 0.0066. Maybe we have to worry about such a coin!

Assume that you forgot the equation of the binomial distribution, you could still have determined the corresponding p-value quite easily with your computer by “simulation”. Indeed, we can easily simulate series of 50 tosses under H0, i.e. the probability to throw a head is equal to that of a tail. We can then count how often – in the simulations - the number of heads departs as much or more from the expected (=25) than with our actual coin (f.i. 35), i.e. how often |Sim-25|≥|35-25|.
We will use R to generate the corresponding distribution by simulation and compare it with the theoretical one.

# Generate distribution of NR by simulation
###########################################
sample (x=0:1,size=nr, replace=TRUE)

 [1] 1 1 1 1 1 0 1 1 1 0 1 1 0 0 1 0 0 1 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 0 1 0 1 1 1 0 0 0 0
[44] 0 1 1 1 1 0 1

sum(sample (x=0:1,size=nr, replace=TRUE))

[1] 24

replicate(10,sum(sample (x=0:1,size=nr, replace=TRUE)))

 [1] 28 24 18 25 25 25 23 21 29 18

nr_repl = 100000
sim <- replicate(nr_repl,sum(sample (x=0:1,size=nr, replace=TRUE)))
qplot(sim,binwidth=1)

hist(sim,breaks=seq(-0.5,50.5,1), border = "red")

# Compare distribution of simulated and theoretical F values
############################################################
for (i in 0:nr){
  toss[i+1,6] = sum(sim == i)/nr_repl
}
toss

plot(toss[,1], toss[,2],col = "blue")
lines(toss[,1],toss[,6],col = "blue")

The p-value of an observed number of heads (or tails) under the null can very easily be computed as follows:

# Compute p-value using simple simulation
#########################################
test_nr = 35
sum(sim >= test_nr | sim <= nr-test_nr)/nr_repl

[1] 0.00652

GENETIC ASSOCIATION STUDIES

To further our discussions, I will invite you in one of my own research fields: genetic association studies. We are not all equal in the face of disease: some of us are more susceptible to asthma, others to cardiovascular issues, other to cancers, etc. Part of these differences are due to differences in “inherited predisposition”. This means that there are genetic variants in the population that increase (risk variants) or decrease (protective variants) disease risk of those that carry them. Imagine a genetic polymorphism with two variants (say allele A and allele a). Imagine that the allelic frequencies in the population are f_A and f_a = 1-f_A. Imagine that the a allele is a risk variant for Crohn’s disease. This means that if you carry a (i.e. your genotype is Aa or aa) your chance to develop Crohn’s disease is increased when compared to the AA individuals. As a consequence, if you sample a cohort of Crohn’s patients and measure the frequency of the a allele in this cohort of “cases”, it is expected to be higher (say f_a+δ) than the frequency of the a allele in a cohort of healthy “controls” (say ~f_a). It is interesting to identify genetic polymorphisms in our genome that contribute to differences in inherited disease predisposition in the population, as this is a first step towards identifying the genes whose function they perturb and which are potential drug targets. To test whether a polymorphism of interest is influencing disease predisposition, one will typically collect DNA samples from n_D cases, and n_H controls, genotype all individuals for the polymorphism of interest, and check whether the genotypes are associated with the disease, i.e. perform an association study. Stated more plainly, this simply means that one will compare the frequency of the A and a alleles in cases and controls and check whether there is any evidence that they differ significantly. Note that for a biallelic polymorphism, testing A is the same (albeit mirror image) as testing a: if one goes up by an amount δ, the other has to go down by the same amount.

Good … Let us do this. Imagine that we implement that with 100 cases and 100 controls, that the frequency of allele a is 0.10 in cases and 0.15 in controls. Is that interesting? Is the A/a polymorphism significantly associated with Crohn’s disease and a hence a genuine risk variant that may lead us to a real causative gene? In other words, is δ = 0.15-0.10 = 0.05 significantly different from zero?

Can we use an approach similar to the one used with the coins to evaluate this? I guess we can imagine something similar. Under H0, the frequency of a in cases and controls is the same: f_a. Thus, we can predict the probability distribution of n_D,a cases and n_H,a controls using the same equation as above (based on the binomial distribution). Assuming that what happens in cases and controls is independent, the probability of the combination of events is just the product of the individual probabilities. Each of these combinations generates a δ, and thus we have δ’s probability distribution under the null. Great! One problem though: with the coins we knew that under H0 f_H = f_T = 0.5, but what is f_a under the null? We actually don’t know for sure. But assuming that H0 is true, hence that f_a is the same in cases and controls, our best estimate of the true value of f_a is the frequency of a in cases and controls combined (largest sample we can assemble). Thus, we will generate an estimate f̂_a using both cases and controls. There we go … we have everything we need.

# Distribution of n_a in a sample of 100 individuals
####################################################
f_a = 0.125
nr_ind = 100
nr_gam = 2*nr_ind
ass_bin <- data.frame(NR_a = 0:nr_gam, FREQ=0, F_FREQ=0, DELTA = 0, F_DELTA = 0, P_DELTA = 0, P = 0)
ass_bin[,2] = ass_bin[,1]/nr_gam
ass_bin[,3] = dbinom(ass_bin[,1],nr_gam,f_a,log = FALSE)
ass_bin[,4] = ass_bin[,1]/nr_gam
ass_bin

plot(ass_bin[,2],ass_bin[,3],col="red", xlab = "FREQ", ylab = "F")

# Frequency distribution of delta between two samples of 100 individuals (f.i. cases and ctrls) under the null
################################################################################################################
for (i in 1:nr_gam+1){
  for (j in 1:nr_gam+1){
    ass_bin[abs(ass_bin[i,1]-ass_bin[j,1])+1,5] = ass_bin[abs(ass_bin[i,1]-ass_bin[j,1])+1,5] + (ass_bin[i,3] * ass_bin[j,3])
  }
}
ass_bin

plot(ass_bin[,4],ass_bin[,5],col="red", xlab = "DELTA", ylab = "F")

# Probability to have delta larger or equal than "x"
####################################################
for (i in 1:nr_gam){
  ass_bin[nr_gam-i,6] = ass_bin[nr_gam+1-i,6]+ass_bin[nr_gam-i,5]
}
ass_bin

plot(ass_bin[,4],ass_bin[,6],col="blue",xlab = "DELTA", ylab = "F/P", xlim = c(0,0.1))
points(ass_bin[,4],ass_bin[,5],col="red")

abline(h = 0.05)

In this specific case, we could fall back on the binomial distribution. But this will not always be the case or tractable. One of the things statisticians have done during the 20ies century is to imagine a number of “tricks” that will allow you in lots of circumstances to perform a simple manipulation of your data such that the outcome corresponds to a “statistic” with a known probability distribution. You have heard of Z (standard normal distribution with mean 0 and var/SD 1), chi-squared (sum of squared Z’s), t (Z/sqrt(n)) and F (ratio of variances) statistics. These are all “quantities” that have a known probability distribution under the null.

In the case of the association test, one can for instance generate a contingency table, and compute SUM(EXP-OBS)^2/EXP, which is known to have a chi-squared distribution with one degree of freedom.

# Contingency table
###################
SOE = 2*(1+25/175)
SOE 

[1] 2.285714

p_SOE = pchisq(SOE,1,lower.tail = FALSE)
p_SOE

[1] 0.13057

plot(ass_bin[,4],ass_bin[,6],col="blue",xlab = "DELTA", ylab = "P", xlim = c(0,0.1))
abline(h = 0.05)

points(0.05,p_SOE,col="green",pch = 19,cex = 2)

M<-as.table(rbind(c(30,170),c(20,180)))
M

    A   B
A  30 170
B  20 180

Xsq <- chisq.test(M,correct=FALSE)
Xsq$statistic

X-squared 
 2.285714 

Xsq$p.value

[1] 0.13057

An alternative is to generate a likelihood ratio test (LRT) which also has a chi-squared distribution.

# Likelihood Ratio Test
#######################
L0 = 0.125^(20+30) * 0.875^(170+180)
L1 = 0.15^30 * 0.85^170 * 0.10^20 * 0.90^180
log(L0)

[1] -150.7081

log(L1)

[1] -149.5584

lrt = -2*(log(L0)-log(L1))
lrt

[1] 2.299305

p_lrt = pchisq(lrt,1,lower.tail = FALSE)
p_lrt

[1] 0.1294319

plot(ass_bin[,4],ass_bin[,6],col="blue",xlab = "DELTA", ylab = "P", xlim = c(0,0.1))
abline(h = 0.05)

points(0.05,p_lrt,col="black",pch = 19,cex = 2)

A problem with these “standard” statistic tests is that to conform to the expected probability distribution, a number of conditions need to be met which may not always be the case. This can generate both false negatives, but more worrisomely false positives as well. So, don’t get excited too quickly when you obtain a low p-value ;-).

The game with the binomial distribution is not that easy, and you may not remember the standard statistic tricks too well. Fortunately, in the 21st century, computers really may help you out a lot. How about using in silico simulations as above, to evaluate whether δ = 0.05 is interesting or not? Why not generate genotypes for a case cohort of n_D individuals and a control cohort of n_H individuals using our best estimate of f̂_a and measure δ? If we repeat this a large number of times, we will in effect generate the distribution of δ under the null and directly evaluate how commonly values as large or larger the the δ seen with real data occur under the null.

# Using simulations
###################
f_a = 0.125
sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a))

  [1] 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 [44] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 [87] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[130] 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[173] 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

sum(sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a)))

[1] 23

sum(sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a)))/nr_gam

[1] 0.125

sim_f<-replicate(1000,sum(sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a)))/nr_gam)
mean(sim_f)

[1] 0.12452

sim_d<-replicate(10000,sum(sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a)))/nr_gam-sum(sample(c(0,1),nr_gam,replace = TRUE,prob = c(1-f_a,f_a)))/nr_gam)
mean(sim_d)

[1] -0.0002835

sum(abs(sim_d)>= 0.05)/10000

[1] 0.1395

plot(ass_bin[,4],ass_bin[,6],col="blue",xlab = "DELTA", ylab = "P", xlim = c(0,0.1))
abline(h = 0.05)

points(0.05,sum(abs(sim_d)>= 0.05)/10000,col="red",pch = 19,cex = 2)

EXERCICES

1 Compute the p-value of 38 euros tosses yielding 12 heads under the null.
2. Compute the p-value of finding the a allele at frequency of 0.20 in 150 cases, and at frequency of 0.24 in 125 cases.

#####################
nr_tosses = 38
nr_heads = 12
nr_repl = 10000
sim<-replicate(nr_repl,sum(sample (x=0:1,size=nr_tosses, replace=TRUE)))
sim
(abs(sim-nr_tosses/2) >= abs(nr_heads-nr_tosses/2))
sum(abs(sim-nr_tosses/2) >= abs(nr_heads-nr_tosses/2))
sum(abs(sim-nr_tosses/2) >= abs(nr_heads-nr_tosses/2))/nr_repl
#######################
nr_cases = 150
nr_ctrls = 125
a_cases = 0.20
a_ctrls = 0.24
nr_repl = 10000
f_a = (a_cases*nr_cases+a_ctrls*nr_ctrls)/(nr_cases+nr_ctrls)
sim_d<-replicate(nr_repl,sum(sample(c(0,1),2*nr_cases,replace = TRUE,prob = c(1-f_a,f_a) ) )/2/nr_cases - sum(sample(c(0,1),2*nr_ctrls,replace = TRUE,prob = c(1-f_a,f_a) ) )/2/nr_ctrls)
mean(sim_d)
sum(abs(sim_d)>= abs(a_cases-a_ctrls))/nr_repl
######################

We actually have an alternative to the simulations that may even be better. The question at hand is whether – in the real data – the genotypes are somehow correlated with disease status (i.e. being case vs control). What if we were to “shuffle” the data, meaning that we would assign the available genotypes randomly to any of the available phenotypes. For sure, any association between genotype and phenotype could only be fortuitous under such scenario. Hence repeating such permutations many times and each time measuring δ (as for the simulations) will yield its distribution under the null. With permutations, the data resemble the real data even better than with the simulations. Specifically, for this example, the genotypes in the permutations match the real data even better than the genotypes in simulations.

We will use real data to illustrate this approach. We will first load a “preprocessed” file with genotypes at ~ 420 SNPs for a case-control cohort for IBD totaling 4,000 individuals.

# Load real dataset
###################
GWAS <- read.csv("~/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/ibd_genotypes_76b541e64057bdfdaf86_/GWAS.data", header=FALSE, sep=";")
#View(GWAS)

The rows in the file correspond to individuals. The columns correspond to “id”, “sex” (1 vs 2), “affection status (1 vs 2), and then genotypes at a series of SNPs (0, 1 or 2 copies of allele a). Let’s us have a look at some of the features of the dataset.

# Have a look at the dataset
############################
nr_var <- ncol(GWAS)-3
nr_var

[1] 420

nr_ind_with_sex <- sum(GWAS$V2 == 1 | GWAS$V2 == 2)
nr_ind_with_sex

[1] 4000

nr_males <- sum(GWAS$V2 == 1)
nr_males

[1] 1893

nr_ind_with_phen <- sum(GWAS$V3 == 1 | GWAS$V3 == 2)
nr_ind_with_phen

[1] 4000

nr_healthy_ind <- sum(GWAS$V3 == 1)
nr_healthy_ind

[1] 2496

nr_ind_with_V4_genotype <- sum(GWAS$V4 == 0 |GWAS$V4 == 1 | GWAS$V4 == 2)
nr_ind_with_V4_genotype

[1] 4000

freq_all_alt_V4 <- sum(GWAS$V4)/(2*sum(GWAS$V4 == 0 | GWAS$V4 == 1 | GWAS$V4 == 2))
freq_all_alt_V4

[1] 0.688625

nr_cases_with_F4_genotype <- sum((GWAS$V4 == 0 | GWAS$V4 == 1 | GWAS$V4 == 2) & GWAS$V3 == 2) 
nr_cases_with_F4_genotype

[1] 1504

freq_cases_alt_V4 <- sum(GWAS$V4[GWAS$V3 == 2])/(2*sum((GWAS$V4 == 0 | GWAS$V4 == 1 | GWAS$V4 == 2) & GWAS$V3 == 2))
freq_cases_alt_V4

[1] 0.6884973

nr_ctrls_with_F4_genotype <- sum((GWAS$V4 == 0 | GWAS$V4 == 1 | GWAS$V4 == 2) & GWAS$V3 == 1) 
nr_ctrls_with_F4_genotype

[1] 2496

freq_ctrls_alt_V4 <- sum(GWAS$V4[GWAS$V3 == 1])/(2*sum((GWAS$V4 == 0 |GWAS$V4 == 1 | GWAS$V4 == 2) & GWAS$V3 == 1))
freq_ctrls_alt_V4

[1] 0.6887019

i=5
sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2))

[1] NA

sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2), na.rm = TRUE)

[1] 3996

sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE)

[1] 2492

sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 2, na.rm = TRUE)

[1] 1504

sum(GWAS[,i+3][GWAS[,3] == 1], na.rm = TRUE)/(2*sum((GWAS[,i+3] == 0 |GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE))

[1] 0.9390048

For reasons which you will understand later I will save the original “affection status” as a distinct vector.

# Save original version of "affection status" as a vector
#########################################################
AFF_STATUS <- GWAS$V3
#AFF_STATUS

We will first compute the association (with disease) p-values for the 420 SNPs using both the contingency table approach and the LRT.

# Preparing a data frame
########################
RES_GWAS <- data.frame(SNP = 1:nr_var, CHR = 0, POS=0, NR_GEN_CTR=0, ALT_FR_CTR=0, NR_GEN_CAS=0, ALT_FR_CAS=0, ALT_FR_ALL=0, 
                       CHISQ=0, LP_CHISQ_1=0, LP_CHISQ_2=0, LRT=0, LP_LRT=0, LP_PERM=0, LP_MULTEST=0)
RES_GWAS

# Generating a function computing the p-value of a contingency table using a chi squared test
#############################################################################################
mychi <- function(i){
  ct<-matrix(c(2*RES_GWAS[i,4]*RES_GWAS[i,5],2*RES_GWAS[i,4]*(1-RES_GWAS[i,5]),2*RES_GWAS[i,6]*RES_GWAS[i,7],2*RES_GWAS[i,6]*(1-RES_GWAS[i,7])), nrow=2, byrow=TRUE)
  Xsq <- chisq.test(ct,correct=FALSE)
  Xsq$p.value
}
# Generating and storing summary data and p-values in data frame
################################################################
for (i in 1:nr_var){
  RES_GWAS[i,4] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE)
  RES_GWAS[i,5] = round(sum(GWAS[,i+3][GWAS[,3] == 1], na.rm = TRUE)/(2*RES_GWAS[i,4]),digits=3)
  RES_GWAS[i,6] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 2, na.rm = TRUE)
  RES_GWAS[i,7] = round(sum(GWAS[,i+3][GWAS[,3] == 2], na.rm = TRUE)/(2*RES_GWAS[i,6]),digits=3)
  RES_GWAS[i,8] = round((RES_GWAS[i,4]*RES_GWAS[i,5]+RES_GWAS[i,6]*RES_GWAS[i,7])/(RES_GWAS[i,4]+RES_GWAS[i,6]),digits=3)
  RES_GWAS[i,9] = 2*( (2*RES_GWAS[i,4]*RES_GWAS[i,5] - 2*RES_GWAS[i,4]*RES_GWAS[i,8]) ^2) / (2*RES_GWAS[i,4]*RES_GWAS[i,8])
                + ( (2*RES_GWAS[i,4]*(1-RES_GWAS[i,5]) - 2*RES_GWAS[i,4]*(1-RES_GWAS[i,8])) ^2) / (2*RES_GWAS[i,4]*(1-RES_GWAS[i,8]))
  RES_GWAS[i,10] = -log10(pchisq(RES_GWAS[i,9],1,lower.tail = FALSE))
  RES_GWAS[i,11] = -log10(mychi(i))
  LH1 = (2*RES_GWAS[i,4]*RES_GWAS[i,5])*log(RES_GWAS[i,5]) + (2*RES_GWAS[i,4]*(1-RES_GWAS[i,5]))*log(1-RES_GWAS[i,5]) + (2*RES_GWAS[i,6]*RES_GWAS[i,7])*log(RES_GWAS[i,7]) + (2*RES_GWAS[i,6]*(1-RES_GWAS[i,7]))*log(1-RES_GWAS[i,7])
  LH0 = (2*RES_GWAS[i,4]*RES_GWAS[i,5])*log(RES_GWAS[i,8]) + (2*RES_GWAS[i,4]*(1-RES_GWAS[i,5]))*log(1-RES_GWAS[i,8]) + (2*RES_GWAS[i,6]*RES_GWAS[i,7])*log(RES_GWAS[i,8]) + (2*RES_GWAS[i,6]*(1-RES_GWAS[i,7]))*log(1-RES_GWAS[i,8])
  RES_GWAS[i,12] = 2*(LH1-LH0)
  RES_GWAS[i,13] = -log10(pchisq(RES_GWAS[i,12],1,lower.tail = FALSE))
}
RES_GWAS

# Generate mini Manhattan plot
##############################
plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")
abline(h = -log10(0.05), lty = 4)

points(RES_GWAS[,1],RES_GWAS[,13],col="blue")

# Compare p-values CT vs LRT
############################
plot(RES_GWAS[,11],RES_GWAS[,13],col="red", xlab = "-log(p) CT", ylab = "-log(p) LRT")

As you can see, both approaches yield virtually identical results.

We will now apply the permutation test. Thus, we will re-assign affection status randomly to the 4,000 individuals yet maintaining the same number of cases and controls. By doing this we obviously disconnect the genotypes from the phenotype as discussed above. Any association could now only be fortuitous (although – see below). We will compute a “statistic” that measures the difference in f_a between cases and controls. I choose just to compute delta. We repeat this n_perm times, hence providing us with the distribution of delta under the null. The p-value of the real delta’s under the null then corresponds to the number of permutations yielding delta’s that are larger or equal to the real ones.

The following code may take a long time to run if you have lots of SNPs and/or many individuals. If you prefer not to run it now, you can go to the next “R chunk” to import precomputed results.

# Generating a data frame (with same structure as RES_GAWAS) to store the results of the permutations (SH-uffles)
##################################################################################################################
SH_GWAS <- data.frame(SNP = 1:nr_var, CHR = 0, POS=0, NR_GEN_CTR=0, ALT_FR_CTR=0, NR_GEN_CAS=0, ALT_FR_CAS=0, ALT_FR_ALL=0, 
                       CHISQ=0, LP_CHISQ=0, LP_CHISQ_2=0, LRT=0, LP_LRT=0, LP_PERM=0, LP_MULTEST=0)
# Generate a function that will shuffle the affection status in GWAS,
# generate and store the relevant summary data and p-values
#####################################################################
shuffle <- function(i){
  GWAS[,3] <- sample(AFF_STATUS,size = 4000,replace = FALSE)
  SH_GWAS[i,4] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE)
  SH_GWAS[i,5] = round(sum(GWAS[,i+3][GWAS[,3] == 1], na.rm = TRUE)/(2*SH_GWAS[i,4]),digits=3)
  SH_GWAS[i,6] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 2, na.rm = TRUE)
  SH_GWAS[i,7] = round(sum(GWAS[,i+3][GWAS[,3] == 2], na.rm = TRUE)/(2*SH_GWAS[i,6]),digits=3)
  DELTA = abs(SH_GWAS[i,5]-SH_GWAS[i,7])
  DELTA
}
# Compute distribution of delta under the null and resulting p-value of real delta for each SNP separately
##########################################################################################################
nr_perm=10000
for (i in 1:nr_var) {
  perm <- replicate(nr_perm,shuffle(i))
  nr_pos <- sum(perm >= abs(RES_GWAS[i,5]-RES_GWAS[i,7]))
  if (nr_pos == 0){
    RES_GWAS[i,14] = -log10(1/nr_perm)
  } else {
    RES_GWAS[i,14] = -log10(nr_pos/nr_perm)
  }
}
RES_GWAS

# Compare the p-values obtained with the different approaches
#############################################################
plot(RES_GWAS[,11],RES_GWAS[,14],col="red")
points(RES_GWAS[,13],RES_GWAS[,14],col="green")

points(RES_GWAS[,11],RES_GWAS[,13],col="black")
# Compare the Manhattan plots obtained with the different approaches
####################################################################
plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")

abline(h = -log10(0.05), lty = 4)
points(RES_GWAS[,1],RES_GWAS[,13],col="green")

points(RES_GWAS[,1],RES_GWAS[,14],col="black")
# Write RES_GWAS to file
########################
write.csv(RES_GWAS,"GWAS_RES",row.names = FALSE)
# Read RES_GWAS
###############
#GWAS_RES <- read.csv("~/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/R_INT_STAT/GWAS_RES")
#View(GWAS_RES)

Run this R chunk to import precomputed results and import them in RES_GWAS

# Read RES_GWAS
###############
GWAS_RES <- read.csv("~/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/R_INT_STAT/GWAS_RES")
#View(GWAS_RES)
RES_GWAS[,14] <- GWAS_RES[,14]
#RES_GWAS
# Compare the p-values obtained with the different approaches
#############################################################
plot(RES_GWAS[,11],RES_GWAS[,14],col="red")
points(RES_GWAS[,13],RES_GWAS[,14],col="green")

points(RES_GWAS[,11],RES_GWAS[,13],col="black")

# Compare the Manhattan plots obtained with the different approaches
####################################################################
plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")
abline(h = -log10(0.05), lty = 4)

points(RES_GWAS[,1],RES_GWAS[,13],col="green")
points(RES_GWAS[,1],RES_GWAS[,14],col="black")

In this case, the p-values obtained by permutation are very similar to those obtained with the other approaches [at least for p-values > 1/(number of permutations)] . But this is not always the case. Assuming that the permutation tests was properly executed, it is the most trustworthy. As you can see, a disadvantage of the permutation test is that it may be time consuming with large datasets (as is the case for GWAS).

CORRECTING FOR FIXED EFFECTS IN PERMUTATION TESTS

Finding a significant association between a SNP and a disease suggests that the SNP may perturb the function of a nearby gene (i.e. in cis) hence increasing disease risk. Indeed, uncovering such scenarios is the main reason why these experiments are conducted in the first place. However, one has to be cautious in interpreting the results, as more trivial explanation may account for the observed association. One of these is so-called population stratification. As an example, the case-control cohort may combine samples originating from different countries. The frequency of the a allele may differ between countries as a results of random drift. If the proportion of cases vs controls differs between countries, these two factors combined may cause a difference in the frequency of the a allele between “all” cases and “all” controls (i.e. across countries) that is not reflecting a genuine “causal” association with disease, but rather population stratification. There are several ways to correct for stratification. As we are discussing the permutation test here, we will only describe a way to ajust the permutation test to account for possible “confounding” effects. What one can do in the example cited above (i.e. an association study conducted with samples originating from different countries having different case vs control ratios) is to perform the permutations within country (i.e. within “strata”). By doing so, the permutations remain “affected” by the same stratification problem as the real data. Thus the delta’s obtained from the permuted data will tend to be inflated as well (as a result of stratification). The delta seen with real data will therefore not be seen as particularly exceptional in regard to the distribution of delta’s obtained across the permutations. Hence, the corresponding p-value will not be deemed particularly small nor interesting.

EXERCICE

As an exercice, you may want to write an R script that permutes for instance within sex in the GWAS dataset.

Correcting for multiple testing

PLAYING WITH COINS AGAIN

Let us go back to our one Euro coins. We will each toss our Euro coin 50 times, record the number of heads, and compute the corresponding p-value, i.e. assuming that our Euro is balanced. I assume (there should be around 20 students in the room) that we will observe that - for at least one of us - the p-value turns out to be <= 0.05. Are there thus “unbalanced”" euro coins in circulation? Or not?

Let us ask the computer. We will simulate what we just did in the room for 10, 20, 30, … 100 students and see how many coins are - on average - declared “unbalanced” (p <= 0.05) per “class”. Given that we simulated balanced coins in the computer (and assuming we did this correctly) these apparently unbalanced coins should all be false positives.

# Building a function that mimics series of "heads and tail"
############################################################
nr_tosses = 300
series <- sample(0:1,nr_tosses, replace = TRUE)
series

  [1] 1 1 0 1 1 1 1 1 0 1 0 0 1 0 1 0 0 1 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0
 [37] 0 1 0 1 1 0 0 1 1 0 1 0 1 1 1 0 1 1 1 0 0 1 1 1 0 0 0 0 1 0 0 1 1 1 1 1
 [73] 0 1 0 0 1 0 0 1 0 1 0 0 1 1 1 0 0 0 1 1 0 0 1 0 1 0 1 1 0 1 0 0 0 0 1 1
[109] 1 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 0 1 0 1 1 0 1 1 1 1 0 0 0 0 1 0 1 0 1 1
[145] 1 0 1 0 0 0 0 0 0 1 0 0 0 1 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0
[181] 0 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 1 0 1 0 0 0 1 0 1 0 1 1 0 1 1 1 0 0 1 0
[217] 1 0 1 0 1 0 1 1 0 0 1 1 1 1 0 0 0 1 0 1 1 0 1 1 0 0 1 0 1 1 0 1 1 0 0 1
[253] 1 0 1 1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1
[289] 1 1 0 0 0 1 1 1 0 0 0 0

nr_heads <- sum(series)
nr_heads

[1] 138

bin_test <- binom.test(x = nr_heads, nr_tosses, p = 0.5, alternative = "two.sided")
bin_test <- binom.test(x = sum(sample(0:1,nr_tosses, replace = TRUE)), nr_tosses, p = 0.5, alternative = "two.sided")
bin_p <- bin_test$p.value
bin_p

[1] 1

coin_test <- function(n){
  bin_test <- binom.test(x = sum(sample(0:1,n, replace = TRUE)),n , p = 0.5, alternative = "two.sided")
  bin_p <- bin_test$p.value
  bin_p
}
coin_test(nr_tosses)

[1] 0.2726299

# Applying the test to classes ranging from 120 to 100 students, with 100 repetitions
######################################################################################
NR_FP<-matrix(1:1000, nrow=100, byrow=TRUE)
#NR_FP <-data.frame(ST10=1:1000,ST20=0,ST30=0,ST40=0,ST50=0,ST60=0,ST70=0,ST80=0,ST90=0,ST100=0)
for (i in 1:100){
  for (j in 1:10){
    stud_ser <- replicate(j*10,coin_test(nr_tosses))
    NR_FP[i,j] <- sum(stud_ser <= 0.05)
  }
}
head(NR_FP,3)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    2    1    2    1    0    4    2    3    3     3
[2,]    0    0    1    1    4    1    0    2    2     7
[3,]    1    0    1    1    4    3    2    4    2     7

boxplot(NR_FP)
for (i in 1:10){
  mean_fp <- mean(NR_FP[,i])  
}
mean_fp

[1] 4.42

Thus, even when all coins are balanced, a certain proportion of students will throw a series of heads and tails with p <= 0.05 under the null. In the example, this proportion is of the order of 0.83 students on average per class of 20, and 4 students on average per class of 100. You will notice that this value is slighltly below 5%.

In fact, and by definition, we expect 5% of the realized tests to be significant at the 5% level, i.e. with p <= 0.05. Using this threshold, we expect to have one false-postive every 20 tests, i.e. the Type I error rate.

The reason why we observed a slightly lower Type I error rate in the example reflects some idiosyncrasies of the binomial distribution. If we increase the size of the series (i.e. the number of tosses), we indeed obtain very close to the expected Type Error rate of I with “nominal” p-value threshold of 0.05.

nr_tosses=15000
sum(sample (x=0:1,size=nr_tosses, replace=TRUE))

[1] 7560

replicate(20,2*pbinom(nr_tosses/2-abs(sum(sample (x=0:1,size=nr_tosses, replace=TRUE))-nr_tosses/2),size=nr_tosses,prob=0.5))

 [1] 0.62999710 0.60697883 0.96743573 0.65338035 0.25640378 0.92843492
 [7] 0.01749831 0.45747560 0.43794384 0.36477163 0.72551946 0.24962436
[13] 0.70116258 0.11506006 0.59561237 0.64164430 0.59561237 0.36477163
[19] 0.12679585 0.49796819

psim2<-replicate(20000,2*pbinom(nr_tosses/2-abs(sum(sample (x=0:1,size=nr_tosses, replace=TRUE))-nr_tosses/2),size=nr_tosses,prob=0.5))
hist(psim2,breaks=seq(-0.1,1.1,0.05), border = "red")

The graph shows that the distribution of p-values under the null is “uniform”, i.e. 5% of the tests have p-value 0 < p < 0.05, 5% have p-value 0.05 < p < 0.10, 5% have p-value 0.10 < p < 0.15, etc … And this is exactly what one would expected given the very way the p-values were defined (think about this …).

A very common way to examine the distribution of p-values for a number of tests is by means of “QQ-plots” (or Quantile-Quantile Plots). The following code generates a QQ plot for the same set of simulations.

# Function generating a QQ plot of p-values "borrowed" from the internat
########################################################################
ggd.qqplot = function(pvector, main=NULL, ...) {
  o = -log10(sort(pvector,decreasing=F))
  e = -log10( 1:length(o)/length(o) )
  plot(e,o,pch=19,cex=1, main=main, ...,
       xlab=expression(Expected~~-log[10](italic(p))),
       ylab=expression(Observed~~-log[10](italic(p))),
       xlim=c(0,max(e)), ylim=c(0,max(o)))
  lines(e,e,col="red")
}
ggd.qqplot(psim2)

Clearly, we have to do something about this if we want to avoid generating false positives. If our experiment entails the realization of multiple tests (f.i. one experiment could be a class of 20 students realizing together 20 tests), we need to adjust the threshold value such that false positives are only expected once every 20 experiments. Such a threshold is referred to as an experiment-wide threshold. A test that yields a p-value < to this threshold is said to be experiment-wide significant. A test that yields a p-value < to the uncorrected (or nominal) significance threshold of 0.05 is said to be nominally significant. People sometimes also use suggestive thresholds. The suggestive threshold is expected to be “exceeded” on average once per experiment (under the null). I will not going in a detailed description of the suggestive threshold.

THE FACT THAT MANY SCIENTISTS DO NOT REALIZE THIS SUFFICIENTLY IS ONE OF THE MAIN CAUSES WHY TOO MANY RESULTS REPORTED IN THE LITERATURE ARE FALSE POSITIVES. JUST REALIZE 20 BOGUS TESTS (EXPRIMENTS) EVERY YEAR, APPLY A 5% SIGNIFICANCE THRESHOLD TO ALL YOUR TESTS, AND YOU ARE BOUND TO REPORT ONE FP IN YOUR YEARLY PAPER. I GUESS YOU MUST START TO SEE/FEEL ONE OF THE DIFFICULTIES: WHERE DO EXPERIMENTS BEGIN AND END? DO WE NEED A LIFE-TIME SIGNIFICANCE THRESHOLD?

How do we properly adjust the significance threshold when performing multiple tests?

Assume that we perform n tests, the adjusted threshold can be derived using Sidak’s correction:

0.05 = 1-(1-p_c)^n

hence corresponding to

p_c = 1 - 10^(log10(0.95)/n)

This threshold is more often approximated using the simpler Bonferroni correction:

p_c = 0.05/n

t = 0.05
SID_BON <- data.frame(n=1:1000,SIDAK=0,BONF=0,LSIDAK=0,LBONF=0)
SID_BON[,2] = round(1-10^(log10(1-t)/SID_BON[,1]),digits=10)
SID_BON[,3] = round(t/SID_BON[,1],digits=10)
SID_BON[,4] = -log10(1-10^(log10(1-t)/SID_BON[,1]))
SID_BON[,5] = -log10(t/SID_BON[,1])
SID_BON

plot(SID_BON[,4],SID_BON[,5],col="red", xlab = "L_SIDAK", ylab = "L_BONF")

plot(SID_BON[,1],SID_BON[,4], col="red", xlab = "N_tests", ylab = "L_BONF threshold")

The graphs show that (i) the Bonferroni correction is an excellent approximation, and (ii) the effect of the number of tests (X axis) on the p-value (in -log units) that needs to be applied to control the experiment-wide type I error at 0.05.

EXERCICE

Try to generate a small script that shows that when you apply Sidak’s or Bonferroni’s correction you only obtain false positives (FP) every 20 experiments conducted. Some experiments will yield 0 FP’s, some 1, some 2, some 3, … The threshold is such that the proportion of experiments with 0 FP corresponds to 95% of the experiments.

RETURNING TO GENETIC ASSOCIATION STUDIES

When performing “Genome-Wide Association Studies” or GWAS, one performs (as in all -omic studies) many, many tests (hundreds of thousands to millions). In this case, every SNP is tested for its association with the disease. Applying a nominal threshold of 0.05 is therefore bound to give a FP every 20 SNPs tested. We therefore have to correct the threshold accordingly to control the rate of FPs.

In the previous example we tested 420 SNPs. A simple way is therefore to correct the 0.05 threshold for the realization of 420 tests using either Sidak’s or Bonferroni’s equation.

plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")
abline(h = -log10(0.05))

abline(h = -log10(0.05/(nr_var)), lty = 2)

The plain horizontal line corresponds to the nominal (uncorrected) threshold of 0.05, while the dashed horizontal line corresponds the threshold, Bonferroni-corrected for the realization of 420 tests, that should control the experiment-wide false positive rate to no more than 0.05.

This may, however, be overkill. Indeed, our genome is characterized by a phenomenon called “linkage disequilibrium” or LD. This means that the genotypes of closely linked SNPs are “correlated”, i.e. they can be very similar to each other. Individuals that are heterozygous for the first SNP (say with alleles A and a) will tend to be heterozygous for the second SNP (say with alleles B and b). Individuals that are homozygous AA may tend to be more often bb than expected by chance. In other words, the genotypes at the “A” and “B” SNPs are not independent. At the extreme (“perfect LD”) the genotypes of two adjacent SNPs may be perfectly redundant in the studied cohort(s): if you know one, you can perfectly predict the other. Two SNPs in perfect LD clearly do not correspond to two different tests; they are one and the same. If x such pairs existed in our dataset, the actual number of realized tests would by 420-x, hence yielding a less stringent experiment-wide threshold (often referred to as genome-wide in the case of GWAS).

For two SNPs that are not perfectly redundant, yet still in LD, the outcome of the association test will tend to be correlated. They can not be seen as two independent tests either. Together, they account maybe for 1.5 but not for 2 tests. The question thus becomes: to how many independent tests does our little GWAS corresponds? If we knew that we could apply to appropriate Sidak/Bonferroni correction to determine the correct experiment-wide threshold.

There are a number of ways to achieve this. Permutation testing provides a very interesting way to correct for the effective number of tests performed, as we will now see.

We previously (see above) used a permutation test to evaluate the significance of the association between one SNP and the disease. The affection status was shuffled across the individuals hence disconnecting any possible relationship between genotype and phenotype (= disease). We are going to do the same thing here, but rather than testing one SNP only, we are going to test all the SNPs with the same phenotype permutation. Thus we are going to perform the same multiple tests as we did with the real data. We want to ensure that the correlation structure between tests that existed with the real data (due to LD between the SNPs), is faithfully maintained in the permutations. Thus, we shuffle the phenotypes between the individuals, but the “vectors” of individual genotypes remain unaltered. Thus, if the Aa genotype at SNP “A” was preferentially associated with the Bb genotype at SNP “B” in the real data, this will remain so in the permuted data.

We still want to compare the allelic frequencies between cases and controls in the permutations. But this time, we have to pay attention that our measure be comparable across SNPs. It shouldn’t for instance be dependent on f_a, which delta is. A good idea in this case is to find inspiration in classic statistics, i.e. use the transformations similar to those yielding the chi-squared, t and F test, or even the ensuing p-values. This is what I will be doing here, using the “CT” test.

Thus, for each permutation we will keep the best p-value (or -log(p-value)). If we perform n_perm such permutations, we will therefore end up with a list of n_perm best p-values. We can rank or sort these best p-values. The 5% experiment-wide significance threshold then corresponds to the 5% percentile of this list, i.e. a value that was only “matched” or “exceeded” (meaning an equal or lower best p-value found) for 5% of permutations. [We can do the same thing using -log(p) values, and determine the 95% percentile]. Thus, this is an elegant way to determine the experiment-wide significance threshold that appropriately corrects for the number of tests that were performed and their correlations.

We can use the same approach to determine the corrected p-value for each SNP. One just has to determine in what fraction of the permutations an equally “good” or “better” p-value (or -log(p)) value was obtained. It is noteworthy that accurate estimation of small (experiment-wide) p-values by permutations requires running a commensurate number of them.

Let us illustrate this with real data.

# Previous mini Manhattan plot with nominalk and Sidak/Bonferroni corrected significance thresholds
###################################################################################################
plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")
abline(h = -log10(0.05))

abline(h = -log10(0.05/(ncol(GWAS)-3)), lty = 2)

# Generating a data frame (with same structure as RES_GAWAS) to store the results of the permutations (SH-uffles)
################################################################################################################## 
SH_GWAS <- data.frame(SNP = 1:(ncol(GWAS)-3), CHR = 0, POS=0, NR_GEN_CTR=0, ALT_FR_CTR=0, NR_GEN_CAS=0, ALT_FR_CAS=0, ALT_FR_ALL=0, 
                      CHISQ=0, LP_CHISQ=0, LP_CHISQ_2=0, LRT=0, LP_LRT=0, LP_PERM=0)
SH_GWAS

# Define function to compute p-value based on CT test
#####################################################
mychi_sh <- function(i){
  ct<-matrix(c(2*SH_GWAS[i,4]*SH_GWAS[i,5],2*SH_GWAS[i,4]*(1-SH_GWAS[i,5]),2*SH_GWAS[i,6]*SH_GWAS[i,7],2*SH_GWAS[i,6]*(1-SH_GWAS[i,7])), nrow=2, byrow=TRUE)
  Xsq <- chisq.test(ct,correct=FALSE)
  Xsq$p.value
}
# Compute relevant stats for all SNPs (We will use the CT test in this example)
###############################################################################
for (i in 1:nr_var){
  SH_GWAS[i,4] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE)
  SH_GWAS[i,5] = round(sum(GWAS[,i+3][GWAS[,3] == 1], na.rm = TRUE)/(2*SH_GWAS[i,4]),digits=3)
  SH_GWAS[i,6] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 2, na.rm = TRUE)
  SH_GWAS[i,7] = round(sum(GWAS[,i+3][GWAS[,3] == 2], na.rm = TRUE)/(2*SH_GWAS[i,6]),digits=3)
  SH_GWAS[i,11] = -log10(mychi_sh(i))
}
SH_GWAS

# Define function to generate shuffled data, compute the corresponding stats FOR ALL SNPs (vs previous function), 
#and extract the top -log(p) value acroos all SNPs
#################################################################################################################
shuffle_2 <- function(){
  GWAS[,3] <- sample(AFF_STATUS,size = 4000,replace = FALSE)
  for (i in 1:nr_var){
    SH_GWAS[i,4] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 1, na.rm = TRUE)
    SH_GWAS[i,5] = round(sum(GWAS[,i+3][GWAS[,3] == 1], na.rm = TRUE)/(2*SH_GWAS[i,4]),digits=3)
    SH_GWAS[i,6] = sum((GWAS[,i+3] == 0 | GWAS[,i+3] == 1 | GWAS[,i+3] == 2) & GWAS[,3] == 2, na.rm = TRUE)
    SH_GWAS[i,7] = round(sum(GWAS[,i+3][GWAS[,3] == 2], na.rm = TRUE)/(2*SH_GWAS[i,6]),digits=3)
    SH_GWAS[i,11] = -log10(mychi_sh(i))
    assign("SH_GWAS",SH_GWAS,envir = globalenv())
  }
  MAX <- max(SH_GWAS[,11])
  MAX
  #list(ran.in = environment(),parent = parent.env(environment()), objects = ls.str(environment()))
}
# Note that the outcome of the loop needs to be assigned to GWAS in globalenv
# Perform permutations and generate list of best p-values
#########################################################
nr_perm=100
mult_test <- replicate (nr_perm,shuffle_2())
mult_test

  [1] 7.781412 1.789614 2.404158 2.939171 1.856394 2.309214 3.151856 1.973791
  [9] 2.534976 3.410673 2.990641 1.534577 1.906961 2.242938 2.177275 1.593806
 [17] 2.702784 4.238273 1.719955 2.249026 2.430571 4.190491 2.007094 2.539132
 [25] 2.705255 1.960397 3.312652 2.640095 2.742447 2.007550 2.213136 1.851350
 [33] 2.730949 1.991178 2.960982 3.926220 1.905751 3.039103 1.893748 1.840908
 [41] 2.335767 2.573776 2.817796 2.832877 2.204520 2.119171 2.061094 2.122357
 [49] 1.626046 2.750356 1.582762 1.717207 3.577824 2.816531 1.980708 3.279910
 [57] 2.977770 4.502004 2.888258 2.532667 3.150311 2.328889 2.594226 2.721036
 [65] 1.425274 2.172181 2.072350 1.954385 2.453251 3.113615 2.594229 1.665414
 [73] 2.885086 1.954591 2.147106 1.656673 2.601734 1.443936 1.748885 1.819391
 [81] 2.223627 2.334531 1.945039 3.982961 2.442108 2.306314 2.256293 1.817003
 [89] 2.042623 2.055231 3.390739 2.165115 1.680421 1.884140 1.872650 2.257531
 [97] 2.284148 2.188978 2.810305 2.835514

# Extract experiment-wide threshold
###################################
mt_thresh <- quantile(mult_test,0.95)
mt_thresh

     95% 
3.929058 

# Add new threshold to mini Manhattan plot
##########################################
plot(RES_GWAS[,1],RES_GWAS[,11],col="red", xlab = "ORDERED SNPs", ylab = "-log(p)")
abline(h = -log10(0.05))

abline(h = -log10(0.05/nr_var), lty = 2)
abline(h = mt_thresh, lty = 4)

# Determine corrected p-values for indivdiual SNPs
##################################################
x=5.495316e+00
mult_test

  [1] 7.781412 1.789614 2.404158 2.939171 1.856394 2.309214 3.151856 1.973791
  [9] 2.534976 3.410673 2.990641 1.534577 1.906961 2.242938 2.177275 1.593806
 [17] 2.702784 4.238273 1.719955 2.249026 2.430571 4.190491 2.007094 2.539132
 [25] 2.705255 1.960397 3.312652 2.640095 2.742447 2.007550 2.213136 1.851350
 [33] 2.730949 1.991178 2.960982 3.926220 1.905751 3.039103 1.893748 1.840908
 [41] 2.335767 2.573776 2.817796 2.832877 2.204520 2.119171 2.061094 2.122357
 [49] 1.626046 2.750356 1.582762 1.717207 3.577824 2.816531 1.980708 3.279910
 [57] 2.977770 4.502004 2.888258 2.532667 3.150311 2.328889 2.594226 2.721036
 [65] 1.425274 2.172181 2.072350 1.954385 2.453251 3.113615 2.594229 1.665414
 [73] 2.885086 1.954591 2.147106 1.656673 2.601734 1.443936 1.748885 1.819391
 [81] 2.223627 2.334531 1.945039 3.982961 2.442108 2.306314 2.256293 1.817003
 [89] 2.042623 2.055231 3.390739 2.165115 1.680421 1.884140 1.872650 2.257531
 [97] 2.284148 2.188978 2.810305 2.835514

(mult_test >= x)

  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [97] FALSE FALSE FALSE FALSE

sum(mult_test >= x)

[1] 1

sum(mult_test >= x)/nr_perm

[1] 0.01

-log10(sum(mult_test >= x)/nr_perm)

[1] 2

for (i in 1:nr_var){
  if (sum(mult_test >= RES_GWAS[i,11]) == 0){
    RES_GWAS[i,15] = -log10(1/nr_perm)
  } else {
    RES_GWAS[i,15] = -log10(sum(mult_test >= RES_GWAS[i,11])/nr_perm)
  }
}
RES_GWAS

# Add corrected p-values to Manhattan plot
###########################################
points(RES_GWAS[,1],RES_GWAS[,15],col="black")

Exploiting the distribution of observed p-values to estimate the False Discovery Rate or FDR

In the previous sections, we came to realize that multiple testing is an important issue that needs to be addressed by adjusting the significance thresholds accordingly. To do this, we determined the (effective) number of test performed (one way or the other) and then corrected the thresholds (one way or the other).

A quite radically different and interesting approach consists in examining the distribution of nominal p-values, see how it compares with the uniform distribution expected if all performed tests are real H0 test (see above), and learn from it.

Let us go back to our Euro coins, and first re-generate the uniform distribution of p-values expected if all coins are balanced. As above, I will test 20,000 coins by tossing each of them 15,000 times.

# Generating the distribution of p-values for 20,000 balanced coins
###################################################################
nr_dr=15000
sum(sample (x=0:1,size=nr_dr, replace=TRUE))

[1] 7643

replicate(20,2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE))-nr_dr/2),size=nr_dr,prob=0.5))

 [1] 0.85103876 0.33942574 0.48766961 0.29975865 0.61844117 0.76257415
 [7] 0.20566566 0.86385944 0.21157597 0.27750486 0.47748708 0.04041883
[13] 0.75015721 0.26331002 0.33123434 0.19987625 0.96743573 0.05100267
[19] 0.55114826 0.33123434

psim2<-replicate(20000,2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE))-nr_dr/2),size=nr_dr,prob=0.5))
hist(psim2,breaks=seq(-0.1,1.1,0.05), border = "red", xlab = "Binned p-values")

We will now repeat the experiment with 20,000 unbalanced coins. The probability to obtain heads is slighly increased to 0.51 (instead of 0.5 for balanced coins).

# Generating the distribution of p-values for 20,000 unbalanced coins
#####################################################################
nr_dr=15000
sum(sample (x=0:1,size=nr_dr, replace=TRUE, prob = c(0.49,0.51)))

[1] 7555

replicate(20,2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE, prob = c(0.49,0.51)))-nr_dr/2),size=nr_dr,prob=0.5))

 [1] 2.091127e-06 2.371302e-02 1.409665e-03 1.227846e-01 1.122022e-03
 [6] 2.847943e-01 4.326498e-04 2.794332e-04 1.463076e-02 3.884984e-02
[11] 3.747347e-03 2.450887e-03 1.041960e-01 2.176697e-02 2.371302e-02
[16] 2.272190e-02 1.832235e-01 8.200710e-02 3.374962e-03 4.606238e-03

psim3<-replicate(20000,2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE, prob = c(0.49,0.51)))-nr_dr/2),size=nr_dr,prob=0.5))
hist(psim3,breaks=seq(-0.1,1.1,0.05), border = "red", xlab = "Binned p-values")

The distribution is severely shifted towards lower p-values suggesting either that our statistical test has a problem, or (as is the case here) that many of our tests (in this case all) are true H1 hypotheses.

Let us know simulate a situation that resembles what we often face in the real world, i.e. that the tests performed include both true H0 and true H1. The proportion of true H0 is often referred to as “pi0”. In the following example “pi0” is simulated to be approximately 50%.

# Generating the distribution of p-values for a mixture of balanced and unbalanced coins
########################################################################################
mixt <- function(){
  if (runif(1,0,1) > 0.5){
    2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE))-nr_dr/2),size=nr_dr,prob=0.5)
  } else{
    2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE, prob = c(0.49,0.51)))-nr_dr/2),size=nr_dr,prob=0.5)
  }
}
nr_dr=15000
nr_exp=20000
psim4 <- replicate(nr_exp,mixt())
hist(psim4,breaks=seq(-0.1,1.1,0.05), border = "red", xlab = "Binned p-values")

The distribution remains heavily shifted towards low p-values due to the high proportion of H1 hypotheses amongst the tests. We normally expect around 1,000 observations per 5% bin under the null. This is approximately the case for bins with p-values > 0.3 in our example. But bins with p-values < 0.3 tend to be over-represented. The first bin, for instance,contains > 7,000 observations, while only 1,000 are expected under the null. Thus, we can reasonably conclude that in this bin approximately 6 out of 7 tests probably correspond to true H1 hypotheses. Or stated otherwise, if we select the tests in the first bin as potentially interesting (we tend to reject the null hypothesis), given the observed distribution of p-values, the proportion of false positives in this bin will be approximately 1 in 7. The right way to say this is that, in the first bin, the False Discovery Rate (FDR) has to be approximately 1/7.

Let us compute this a bit more accurately:

# Computing the FDR for the group of tests in the 0-0.05 bin
############################################################
sum(psim4 <= 0.05)

[1] 7294

FDR = (nr_exp*0.05)/sum(psim4 <= 0.05)
FDR

[1] 0.137099

Another common way to look at this distributions is by means of the QQ plots that were already discussed before.

# QQ plot for mixture of balanced and unbalanced coins
######################################################
ggd.qqplot = function(pvector, main=NULL, ...) {
  o = -log10(sort(pvector,decreasing=F))
  e = -log10( 1:length(o)/length(o) )
  plot(e,o,pch=19,cex=1, main=main, ...,
       xlab=expression(Expected~~-log[10](italic(p))),
       ylab=expression(Observed~~-log[10](italic(p))),
       xlim=c(0,max(e)), ylim=c(0,max(o)))
  lines(e,e,col="red")
}
ggd.qqplot(psim4)

In the previous example, we computed the FDR for the group of tests falling in the lowest bin. We can also compute the FDR for individual tests. These values correspond to so-called q-values. This is in fact a bit misleading. As you will soon understand, the q-value computed for an individual test is actually the joint FDR for all the test with p-value equal or smaller than the examined one. Let us compute the q-values for the “psim4” (mixed balanced and unbalanced coins) dataset.

# Computing anb examining q-values
###################################
#psim4
rank_sim4<-rank(psim4,na.last = TRUE, ties.method = "random")
#rank_sim4
#sort_sim4<-sort(psim4,decreasing=FALSE)
#sort_sim4
fdr<-(nr_exp*psim4)/rank_sim4
lpsim4<-log10(1/psim4)
lfdr<-log10(1/fdr)
#plot(lpsim4,lfdr)
#bonf<-psim4*nr_exp
#lbonf<-log10(1/bonf)
#plot(lbonf,lpsim4)
sid<-1-(1-psim4)^nr_exp
lsid<-log10(1/sid)
#plot(lsid,lpsim4)
plot(lsid,lfdr)

[+ interpretation]

As a matter of fact we can do better. From the distribution of p-values, we can have an idea of “pi0”, i.e. the proportion of H0 hypotheses amongst the realized tests. We will assume for that, that the “high p-value bins”" primarily encompass H0 hypotheses. Indeed, typically the height of these bins “flattens out”. Thus the average proportion of observations in these bins multiplied by the total number of bins (say 20) gives us an estimate of the proportion of H0 hypotheses in the data. This estimate can then be used to compute improved q-values. The idea behind this approach is described in Storey et al. (2003)[https://www.ncbi.nlm.nih.gov/pubmed/12883005].

Let us implement it to our psim4 data, using Storey’s software, which we will first install.

# Installing Storey's package
#############################
#install.packages("devtools")
library("devtools")
#install_github("jdstorey/qvalue")
library(qvalue)

While playing with this I realized that psim4 contains values > 1. We will first correct these values and set them to 1.

# Eliminating p-values > 1 from psim4
#####################################
max(psim4)

[1] 1.006515

psim4[psim4>1]<-1
max(psim4)

[1] 1

We will then recompute FDR’s using “our” method.

# Computing FDR with "our" method
##################################
rank_sim4<-rank(psim4,na.last = TRUE, ties.method = "random")
#rank_sim4
fdr<-(nr_exp*psim4)/rank_sim4
lfdr<-log10(1/fdr)

Let us now use Storey’s program.

# Computing FDR with Storey's method
#####################################
qobj <- qvalue(p = psim4)
#qobj
lq<-log10(1/qobj$qvalues)
lp<-log10(1/qobj$pvalues)

We can now compare the two sets of values.

plot(lfdr,lq)
lines(lfdr,lfdr,col = "red")

We can also use Storey’s program to estimate pi0 values.

# Generating the distribution of p-values for a mixture of balanced and unbalanced coins
########################################################################################
mixt <- function(pi_zero){
  if (runif(1,0,1) < pi_zero){
    2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE))-nr_dr/2),size=nr_dr,prob=0.5)
  } else{
    2*pbinom(nr_dr/2-abs(sum(sample (x=0:1,size=nr_dr, replace=TRUE, prob = c(0.49,0.51)))-nr_dr/2),size=nr_dr,prob=0.5)
  }
}
nr_dr=15000
nr_exp=10000
psim_series <- data.frame(V0=1:nr_exp,V25=0,V50=0,V75=0,V100=0)
#psim_series
for (i in 1:5){
  psim_series[,i] <- replicate(nr_exp,mixt(pi_zero = 0.10*(i-1)))
}
#psim_series
pi0 = data.frame(V1=0:1,V2=NA,V3=NA,V4=NA,V5=NA)
for (i in 1:5){
  pi0[1,i] = 0.10*(i-1)
  psim_series[,i][psim_series[,i]>1]<-1
  qobj <- qvalue(p = psim_series[,i])
  pi0[2,i]<-qobj$pi0
}
pi0

In these simulations, Storey’s program seem to slightly but systematically overestimate pi0 … :-)

Examples of the design and use of empirical statistical test in genome studies

1. Rosewick, N. et al. Cis-perturbation of cancer drivers by the HTLV-1/BLV proviruses is an early determinant of leukemogenesis. Nat Commun 8:15264 (2017) [https://www.nature.com/articles/ncomms15264]

This paper makes extensive use of simulations to evaluate the statistical significance of genomic features associated with BLV proviral insertion sights. The BLV and HTLV-I retroviruses cause leukelia in cattle/sheep and human, respectively, and this after a long incubation period. It was previously assumed that the proviral integration site was random, and that tumor development was due to (i) the transformation potential of virus-encoded proteins including Tax in the early phases of the disease, and (ii) somatic mutations accruing in cancer driver genes across the genome in the later stages of the disease. This paper shows that the integration of BLV proviruses in sheep is absolutely not random. Proviruses tend to preferentially “land” in hotspots of integration. Moreover, these hotspots are enriched in known cancer driver genes, strongly suggesting that the integration site contributes to early cell transofmration. Simulation was used in several experiments to evaluate the statistical significance of the observed proviral clustering.

Identifying integration hotspots: “Hotspots were then defined as follows: for each IS we counted the number of IS located in a 50 kb window centred on that IS. We then randomly picked 66,557 IS in the ovine genome, performed the same counting procedure and retained the simulated IS window harbouring the highest number of IS. We performed N=1,000,000 iterations of this procedure to generate a distribution reflecting random integration. We then assigned a P-value to each IS. Hotspots were defined by grouping consecutive IS windows that showed a P-value <= 0.05, the position of the two extreme IS 5 kb determining the boundaries of the hotspot.”

2. Rossin, E.J. et al. Proteins encoded in genomic regions associated with immune-mediated disease physically interact and suggest underlying biology. PLoS Genetics 7: e1001273 (2011) [https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1001273]

“Therefore, we apply a permutation method that is robust to non-specific binding and differences in publication density. We perform a within-degree node-label permutation that is carried out as follows: a random network is built that has nearly the exact same structure as the original InWeb network, only the node labels (i.e. the protein names) are randomly re-assigned to nodes of equal binding degree; this method assumes a null distribution of connectivity that is entirely a function of the binding degree of individual proteins. Random networks will have the same size, number of edges and per-protein binding degree as InWeb; we build 50,000 different random networks. With this method, we are able to test the non-randomness of our network conditional on the exact binding degree distribution of our disease proteins.”

The winners’ curse

In the previous section we showed how not properly accounting for multiple testing is a common source of false positive, i.e. non reproducible results. Another possible cause of the apparent non-reproducibility of results is when experiments are “underpowered”. This is well illustrated in GWAS studies. Several groups around the world initially performed “their own” GWAS, i.e. they assempled a case-control cohort of - say - 1,000 cases and 1,000 controls - and ran their own analyses. This would typically yield one or two genome-wide significant signals. All of these discoveries were accompanied by a confirmation in an independent “confirmation” cohort, hence strongly supporting their genuine nature [Note that the significance thresholds used in the confirmation step are “leanier” as much fewer tests are conducted in the confirmation than in the discovery phase]. Yet, the signals that were discovdered in these initial GWAS nearly always differed between studies. Why is that?
The reason is that the SNP effects were smaller than anticipated. As a result, the initial GWAS studies were initially “underpowered” to detect them. Assuming that you only have 10% power to detect the signals of interest given your sample size, you will by definition only find 10% of the existings signal in your study. Your colleague, who is in the same situation, may by chance detect another of the existing signals in hers/his study, etc. What lies behind the “by chance” statement in this context? Well to detect a signal in an underpowered study, the corresponding effect has to be “accidentally” slightly inflated in that study. Otherwise it would not exceed the significance threshold. This is the essence of what is commonly referred to as the winners’ curse. That the effects are indeed inflated in the discovery phase is confirmed when re-estimating the effect in the confirmation cohort: they are nearly always inferior to what they were in the discovery phase. All this made people quickly realize that salvation lays in collaboration. A new layer of effects that would be impossible to detect in individual studies indeed emerged when performing meta-analayses by merging the indiviual datasets. This is indeed what happened for nearly all studied diseases, yielding an amazing “crop” of GWAS hits for nearly all of them [see EBI’s GWAS catalog at https://www.ebi.ac.uk/gwas/].

Let us explore this a bit more quantitatively. We are going to simulate the identification, by association study, of a genetic variant affecting a quantitative trait. For simplicity we are going to assume that we work with a haploid organism. Thus the genotype of individuals is either A (“0” in the simulations) or a (“1” in the simulations). We will set the frequency of the a allele at 0.25 (the script below allows to vary f_a as well as to test a range of values). We will assume that the trait of interest is normally distributed with mean 0 and standard deviation 1 within genotype (i.e. “residual” variance). We will assume that the effect of the a allele is to increase the phenotype by delta. In our simulations, delta will vary from say 0.25 to 1 residual standard deviation. We will vary the sample size from, say 200 to >1,400. We will test the significance of the effect in our simulated data set using a t-test. We will select a significance threshold (f.i. p < 0.000001) to decide whether delta is significant or not. We will examine the impact of the studied parameters on (i) the power to detect the signal at the chosen significance threshold, (ii) the impact of power on the estimate of delta.

The following code may take a long time to run if run large numbers of simulations. If you prefer not to run it now, you can go to the next “R chunk” to import precomputed results.

# Set all_comb at zero prior to starting
########################################
all_comb = NULL
# Select parameter (allele freqiency, size of the effect, sample size) space which you want to explore
######################################################################################################
freq<-c(0.25)
effect<-c(0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,1)
ssize<-c(200,400, 600, 800, 1000, 1200, 1400)
# Generate all_comb matrix using "expand" such that all possible paremeter combinations correspond to one row
#############################################################################################################
all_comb<-expand.grid(freq,effect,ssize)
all_comb[,4:9]=0 # I just put some columns in which I will do additions at zero. May not be necessary 
names(all_comb)=c("AF","EFF","SS","AD","SUCC","DS","ADP","SUCCP","DSP")
# Select significance threshold to use
######################################
p_tresh=0.000001
# Generate function to simulate a population with given values to be found in all_comb and perform its analysis
###############################################################################################################
run_sim<-function(j){
  sim=NULL
  # Generate a data.frame in which to store individual information
  sim<-data.frame(GEN=1:all_comb[j,3],RES=NA,PHEN_0=NA,PHEN_1=NA)
  sim$GEN=sample(x=0:1,size=all_comb[j,3],replace=TRUE, prob = c(1-all_comb[j,1],all_comb[j,1]))
  sim$RES=rnorm(all_comb[j,3],0,1)
  sim$PHEN_0 = sim$RES
  sim$PHEN_1 = sim$RES+all_comb[j,2]
  sim$PHEN_0[sim$GEN==1] = NA
  sim$PHEN_1[sim$GEN==0] = NA
  # Calculate realized delta
  delta = mean(sim$PHEN_1,na.rm = TRUE) - mean(sim$PHEN_0,na.rm = TRUE)
  # Store realized data in appropriate place in all_comb
  all_comb[j,4]=all_comb[j,4]+delta
  # test significance of delta using t-test    
  t_res<-t.test(sim$PHEN_0,sim$PHEN_1, alternative = "two.sided")
  # Record whether significant and if it is, store results in appropriate place in all_comb
  if (t_res$p.value < p_tresh){
    all_comb[j,5]=all_comb[j,5]+1
    all_comb[j,6]=all_comb[j,6]+delta
  } 
  # Reset SIM
  sim=NULL
  # Transfer results to globalenv
  assign("all_comb",all_comb,envir = globalenv())
}
# Define numner of simulations to perform and run these one row at a time
#########################################################################
nr_sims = 20000
for (j in 1:nrow(all_comb)){
  replicate(nr_sims,run_sim(j))
}

Run this R chunk to import precomputed results

# Load precomputed all_comb fine
################################
all_comb <- read.csv("~/Documents/Teaching/GIGA_Doctoral_School/Intuitive_statistics/R_INT_STAT/WIN_CURSE")
#View(all_comb)
# Generate plots
################
plot(all_comb[1:9,2],all_comb[1:9,2],type = "l", xlab = "Real effect",ylab = "Estimated effect / Power", col = "red")
lines(all_comb[1:9,2],all_comb[1:9,9], col = "#08306B")

lines(all_comb[10:18,2],all_comb[10:18,9], col = "#08519C")
lines(all_comb[19:27,2],all_comb[19:27,9], col = "#2171B5")

lines(all_comb[28:36,2],all_comb[28:36,9], col = "#4292C6")
lines(all_comb[37:45,2],all_comb[37:45,9], col = "#6BAED6")

lines(all_comb[46:54,2],all_comb[46:54,9], col = "#9ECAE1")
lines(all_comb[55:63,2],all_comb[55:63,9], col = "#C6DBEF")

lines(all_comb[1:9,2],all_comb[1:9,8], col = "#08306B", lty = 2)
lines(all_comb[10:18,2],all_comb[10:18,8], col = "#08519C", lty = 2)

lines(all_comb[19:27,2],all_comb[19:27,8], col = "#2171B5", lty = 2)
lines(all_comb[28:36,2],all_comb[28:36,8], col = "#4292C6", lty = 2)

lines(all_comb[37:45,2],all_comb[37:45,8], col = "#6BAED6", lty = 2)
lines(all_comb[46:54,2],all_comb[46:54,8], col = "#9ECAE1", lty = 2)

lines(all_comb[55:63,2],all_comb[55:63,8], col = "#C6DBEF", lty = 2)

The dashed lines correspond to the power to detect the simulated effects (X axis) with population of increasing size (dark to light blue). The red diagonal corresponds to the expectation assuming that the effects are perfectly estimated. The regular lines correspond to the estimates effects in the populations of increasing size (dark to light blue). What we see is that when the power departs from 100%, the effects tend to be overestimated. Quite interestingly, at suboptimal power the estimated effects appear to be determined by the sample size more than by the actual effect. Thus, under these conditions, the estimates tell us more about our experimental design than about the biology we are trying to study.

How can we address this “curse”? Especially in -omic studies, it is a good idea to distinguish a discovery phase from a confirmation phase. Effects are first discovered in a “discovery” sample, applying stringent significance thresholds that properly account for multiple testing. In a second stage, the “discovered” effects are then confirmed in an independent “confirmation” sample, using a more lenient threshold that only corrects for the number of discoveries. If significant in the confirmation sample, the effect is deemed real, while the conformation cohort provides a less biased estimate of the size of the true effect.

Take home messages

1. The p-value of a test is the probability to have an equally or more unlikely outcome of a test under the null hypothesis (H0).

2. p-values can be computed with standard statistic tests or, increasingly, using self-designed empirical tests using in silico simulation and/or permutation.

3. When performing multiple tests, the threshold p-value to declare significance has to be adjusted in order to control the number of experiment-wide false positives. NOT DOING THIS APPROPRIATELY IS A MAJOR SOURCE OF FALSE POSITIVES IN THE LITTERATURE.

4. Rather than adjusting the threshold value for the effective number of realized tests, one may exploit the obtained distribution of p-values to compute FDR or q-values. This is especially effective when the proportion of true H1 hypotheses is high amongst the realized tests, hence causing a major shift of the p-value distribution to low values.

5. When the statistical power to detect a true H1 is limited (for instance because the sample size is too small), the effect will often be missed (false negatives, type II error). THIS IS A VALID CAUSE OF NON-REPRODUCIBILITY OF PUBLISHED RESULTS. To detect the effect one has to be a bit lucky, i.e. the effect needs to be - by chance - overestimated in your data. Thus, when one detects an effect that one had a small chance to detect, it will be overestimated in the “discovery” sample (= Winners’ curse). As matter of fact, under conditions of low statistical power, the estimates of the effects are determined more by the sample size than by the actual effect size.

Acknowledgements

I thank the 2019 students for a very pleasant and interactive three days (21-23/10). I would also like to thank one of the students (2019, Renaud Nivelle) for sharing his excellent R scripts, and helping us all to progress.