DU Bii - module 3: R and stats

---

Session 3: practical on statistics for omics data

Thursday 4th of March, 2021 Claire Vandiedonck

teachers: Claire Vandiedonck & Antoine Bridier-Nahmias; helpers: Anne Badel, Clémence Réda, Olivier Sand, Jacques van Helden

Content of this tutorial:

1. Test statisqtique de comparaison de moyennes: impact de la différence entre les moyennes et de la variance
2. Effets des tests multiples sur les taux de faux positifs

=> About this jupyter notebook

This a jupyter notebook in R, meaning that the commands you will enter or run in Code cells are directly understood by the server in the R language.
You could run the same commands in a Terminal or in RStudio.

In this tutorial, you will run one cell at a time.

Do not hesitate to try other commands by adding other cells.

0. Before going further

---

Check your your working directory, and sit it to another path if needed:

getwd()

'/shared/ifbstor1/home/cvandiedonck/dubii2021/RSession3'

#setwd('/shared/home/cvandiedonck/RSession3') #change with your login!!!
#getwd() #change is visible

I. Test statistique de comparaison de moyennes: impact de la différence entre les moyennes et de la variance

---

I.A. Deux populations avec une moyenne différente mais une même variance

Nous allons d'abord tirer deux échantillons de deux populations présentant des moyennes différentes et une même variance.

A titre d'exemple dans ce parctical nous prendrons la "taille des humains debout" comme trait d'intérêt. En France en 2001, les données de l'INSEE indiquaient une taille moyenne de 1.62 pour les femmes et de 1.74 m pour les hommes. Dans l'ensemble de la population, l'écacrt-type était d'environ 7 cm.

• Générez 2 échantillons sample1 et sample2 de taille 100 issus de la population de femmes et de la population d'hommes avec les paramètres indiqués en supposant que la densité de probabilité de la taille suit une loi normale.

Tip: utiliser la fonction rnorm()

sample1 <- rnorm(n=100, mean=162, sd=7)

sample2 <- rnorm(n=100, mean=174, sd=7)

• Regardez les valeurs résumées des distributions obersvées et faites un simple boxplot de ces distributions.

Tip: utiliser les fonctions summary() et boxplot

summary(sample1)
summary(sample2)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  147.2   158.2   163.3   162.5   166.9   184.4 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  157.0   167.8   172.8   172.4   177.3   193.4 

boxplot(sample1, sample2)

• Testez la différence des moyennes entre les deux populations d'où sont issus ces échantillons avec un test de Student bilatéral

Tip: utiliser la fonction t.test()

t.test(sample1, sample2, alternative = "two.sided")

	Welch Two Sample t-test

data:  sample1 and sample2
t = -10.223, df = 197.96, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.779396  -7.969723
sample estimates:
mean of x mean of y 
 162.4941  172.3687 

• Quelle est la valeur de la p-value?

Tip: utilisez la fonction str() pour voir comment extraire cette pvalue du test

str(t.test(sample1, sample2, alternative = "two.sided"))

List of 10
 $ statistic  : Named num -10.2
  ..- attr(*, "names")= chr "t"
 $ parameter  : Named num 198
  ..- attr(*, "names")= chr "df"
 $ p.value    : num 5.72e-20
 $ conf.int   : num [1:2] -11.78 -7.97
  ..- attr(*, "conf.level")= num 0.95
 $ estimate   : Named num [1:2] 162 172
  ..- attr(*, "names")= chr [1:2] "mean of x" "mean of y"
 $ null.value : Named num 0
  ..- attr(*, "names")= chr "difference in means"
 $ stderr     : num 0.966
 $ alternative: chr "two.sided"
 $ method     : chr "Welch Two Sample t-test"
 $ data.name  : chr "sample1 and sample2"
 - attr(*, "class")= chr "htest"

t.test(sample1, sample2, alternative = "two.sided")$p.value

5.72396607010878e-20

Illustrons à présent cette différence à l'aide du paquet dabestR https://github.com/ACCLAB/dabestr

• Installez le paquet dabestr s'il n'est pas déjà installé et chargez-le dans votre session R.

Tip: utilisez la fonction requireNamespace), install.packages(), library() et sessionInfo()

if (!requireNamespace("dabestr"))
  install.packages("dabestr")
library(dabestr)
sessionInfo()

Loading required namespace: dabestr

Loading required package: magrittr

Warning message:
“package ‘magrittr’ was built under R version 4.0.3”

R version 4.0.2 (2020-06-22)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.0.2/lib/libopenblasp-r0.3.10.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dabestr_0.3.0  magrittr_2.0.1

loaded via a namespace (and not attached):
 [1] munsell_0.5.0    tidyselect_1.1.0 uuid_0.1-4       colorspace_2.0-0
 [5] R6_2.5.0         rlang_0.4.10     stringr_1.4.0    dplyr_1.0.3     
 [9] tools_4.0.2      grid_4.0.2       gtable_0.3.0     DBI_1.1.1       
[13] htmltools_0.5.1  ellipsis_0.3.1   assertthat_0.2.1 digest_0.6.27   
[17] tibble_3.0.5     lifecycle_0.2.0  crayon_1.3.4     IRdisplay_0.7.0 
[21] ggplot2_3.3.3    purrr_0.3.4      repr_1.1.0       base64enc_0.1-3 
[25] vctrs_0.3.6      IRkernel_1.1.1   glue_1.4.2       evaluate_0.14   
[29] stringi_1.5.3    pbdZMQ_0.3-3.1   compiler_4.0.2   pillar_1.4.7    
[33] scales_1.1.1     generics_0.1.0   boot_1.3-25      jsonlite_1.7.2  
[37] pkgconfig_2.0.3  Cairo_1.5-12.2  

• Consultez la page github de dabestr voire sa vignette et faites une représentation graphique dite Gardner-Altman estimation plot des distributions qui affiche l'IC de la moyenne en prenant les hommes comme groupe de référence.

Pour cela vous devez d'abord générer un dataframe:

- Rassemblez les données des deux échantillons dans un dataframe avec deux colonnes `height` et `grp` pour le groupe.

Tip: utilisez la fonction data.frame() sur un vecteur concanténant les taille des deux échantillons avec la fonction c() et un vecteur répétant les valeurs 1 et 2 100 fois dans le bon ordre avec la fonction rep() et son argument each=

df <- data.frame(height=c(sample1, sample2), grp=rep(c(1,2), each=100))

- Affichez les premières, dernières lignes du dataframe et sa structure. Comptez également le nombre d'obersvations appartenant au groupe 1 ou 2.

Tip: utilisez les fonctions head() , tail(), str() et table().

head(df)

A data.frame: 6 × 2

	height	grp
	<dbl>	<dbl>
1	162.7547	1
2	160.3425	1
3	158.9158	1
4	164.2203	1
5	152.4980	1
6	163.6377	1

tail(df)

A data.frame: 6 × 2

	height	grp
	<dbl>	<dbl>
195	177.2447	2
196	167.3790	2
197	171.7305	2
198	172.6847	2
199	178.4805	2
200	175.7835	2

str(df)

'data.frame':	200 obs. of  2 variables:
 $ height: num  163 160 159 164 152 ...
 $ grp   : num  1 1 1 1 1 1 1 1 1 1 ...

table(df$grp)

  1   2 
100 100 

- Faites de Grander-Latman estimation plot:

Tip: utilisez les fonctions dabest() , unpaired_mean_diff(), et plot().

unpaired_mean_diff <- dabest(df, grp, height,
                             idx = c("1", "2"),
                             paired = FALSE) %>% 
                      mean_diff()
unpaired_mean_diff
plot(unpaired_mean_diff)

dabestr (Data Analysis with Bootstrap Estimation in R) v0.3.0
=============================================================

Good evening!
The current time is 12:23 PM on Tuesday March 09, 2021.

Dataset    :  df
X Variable :  grp
Y Variable :  height

Unpaired mean difference of 2 (n = 100) minus 1 (n = 100)
 9.87 [95CI  8; 11.8]


5000 bootstrap resamples.
All confidence intervals are bias-corrected and accelerated.

I.B. Impact de la taille des échantillons, de la différence des moyennes ou de la variance

Nous allons à présent étudier l'impact de la différence entre les moyennes des populations de départ ou celui de la variance sur la puissance de détecter une différence.

• Pour jouer sur ces différents paramètres, écrivez une fonction qui prend comme argument la taille des échantillons, la moyenne et la variance de chacune des populations, et qui renvoie la statistique t et la pvalue du test de Student.

compute_tv <- function(n1, m1, s1, n2, m2, s2){
    sample1 <- rnorm(n=n1, mean=m1, sd=s1)
    sample2 <- rnorm(n=n2, mean=m2, sd=s2)
    tval <- t.test(sample1, sample2, alternative = "two.sided")$statistic
    pval <- t.test(sample1, sample2, alternative = "two.sided")$p.value
    tv <- list("t"=tval, "p"=pval)
    return(tv)
}

• Faites le tests en prenant à présent les écarts types différents chez les hommes et les femmes qui sont respectivement de 7.1 cm et 6.5 cm dans la population Française en 2001.

compute_tv(100, 162, 6.5, 100, 174, 7.1)

$t

t: -13.987956846351

$p

4.06669423654191e-31

• Refaites le test en réduisant les effectifs à 30 pour cahque échantillon.

compute_tv(30, 162, 6.5, 30, 174, 7.1)

$t

t: -4.67094879499133

$p

2.05444475523577e-05

=> Comment a varié la significativité?

Les paramètres de la distribution de la taille des humains varient aussi entre les populations dans le monde.

• Refaites le test avec ces effectifs de 30, en gardant la même différence entre les moyennes mais en augmentant la variance à 25 pour chaque population.

compute_tv(30, 162, 25, 30, 174, 25)

$t

t: -1.41882851214672

$p

0.161303380927312

=> La différence, même importante entre les moyennes des populations, est-elle encore détectable?

• Refaites le test avec des populations présentant une petite différence entre les moyennes (3 cm avec par exemple 160 cm pour les femmes et 163 cm pour les hommes) et des variances faibles (3 cm)

compute_tv(30, 160, 3, 30, 163, 3)

$t

t: -2.24359804586757

$p

0.0289520612020721

=> Détectez-vous quand même une différence statistiquement significative?

Conclusion: La puissance du test augmente avec: - la taille des échantillons - la différence entre les moyennes des populations - une petite variance

II. Effet des tests multiples sur le taux de faux positifs

---

Dans les jeux de données omiques, on effectue autant de tests que de "features".

Nous n'allons pas générer dans ce practical de multiples features mais nous allons prendre notre trait de la taille et tester par simulation l'impact des tests multiples sur les faux positifs.

Nous repartons du dataframe généré ci-dessus.

• Nous allons simuler 10000 échantillons sous une hypothèse nulle où le groupe ne modifie pas les moyennes. Pour ce faire, générez 10000 permutations de Monte-Carlo en tirant au hasard les statuts pour les groupes. Et récupérez les valeurs de tests tvalues et pvalues .

Tip: Utilisez la fonction for() pour faire répéter votre code 1000 fois. Dans la boucle, générez un vecteur avec les valeurs perumutées pour le groupe et effectuez le test pour récupérer les valeurs p et de statsitique

tvalues <- c()
pvalues <- c()

for(i in 1:10000){
    perm_df <- df
    perm_df$grp <- sample(perm_df$grp, size=200)
    tval <- t.test(perm_df$height ~ perm_df$grp, alternative = "two.sided")$statistic
    tvalues <- c(tvalues, tval)
    pval <- t.test(perm_df$height ~ perm_df$grp, alternative = "two.sided")$p.value
    pvalues <- c(pvalues, pval)
}

• Affichez les valeurs résumées et un histogramme des statistiques t:

summary(tvalues)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-4.084798 -0.681368  0.007453  0.002084  0.672858  3.900867 

hist(tvalues)

• Combien de valeurs de statsitique sont significatives bien que conformes à H0?

Tip: Utilisez la fonction qnorm() pour retrouver la valeur seuil d'une loi normale au risque alpha = 5%

qnorm(0.025, 0, 1, lower.tail=FALSE)

1.95996398454005

length(which(abs(tvalues)>=1.96))

485

=> Ce résultat était-il attendu?

• Affichez les valeurs résumées et un histogramme des pvalues:

summary(pvalues)
hist(pvalues)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000643 0.2507573 0.4986034 0.4995971 0.7513176 0.9999545 

=> Quelle est la distribution des pvalues sous HO?

Conclusion: Sous HO: - on obtient environ une proportion alpha de valeurs significatives -> si 10000 tests ont été effectués, cela fait envriron 500 différences significatives! Il faut donc corriger pour les tests multiples pour éviter un taux de faux positifs trop élevé. - la distribution des pvalues est une distribution uniforme entre 0 et 1.

---

---

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.0.2/lib/libopenblasp-r0.3.10.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dabestr_0.3.0  magrittr_2.0.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         plyr_1.8.6         vipor_0.4.5        pillar_1.4.7      
 [5] compiler_4.0.2     RColorBrewer_1.1-2 forcats_0.5.0      base64enc_0.1-3   
 [9] tools_4.0.2        boot_1.3-25        digest_0.6.27      uuid_0.1-4        
[13] jsonlite_1.7.2     evaluate_0.14      lifecycle_0.2.0    tibble_3.0.5      
[17] gtable_0.3.0       pkgconfig_2.0.3    rlang_0.4.10       IRdisplay_0.7.0   
[21] DBI_1.1.1          IRkernel_1.1.1     beeswarm_0.2.3     repr_1.1.0        
[25] dplyr_1.0.3        stringr_1.4.0      generics_0.1.0     vctrs_0.3.6       
[29] cowplot_1.1.1      grid_4.0.2         tidyselect_1.1.0   glue_1.4.2        
[33] R6_2.5.0           ggbeeswarm_0.6.0   pbdZMQ_0.3-3.1     farver_2.0.3      
[37] tidyr_1.1.2        purrr_0.3.4        ggplot2_3.3.3      scales_1.1.1      
[41] ellipsis_0.3.1     htmltools_0.5.1    assertthat_0.2.1   colorspace_2.0-0  
[45] labeling_0.4.2     stringi_1.5.3      munsell_0.5.0      crayon_1.3.4      
[49] Cairo_1.5-12.2    

---

---

About jupyter notebooks:
- To add a new cell, click on the "+" icon in the toolbar above your notebook
- You can "click and drag" to move a cell up or down
- You choose the type of cell in the toolbar above your notebook:
- 'Code' to enter command lines to be executed
- 'Markdown' cells to add text, that can be formatted with some characters
- To execute a 'Code' cell, press SHIFT+ENTER or click on the "play" icon
- To display a 'Markdown' cell, press SHIFT+ENTER or click on the "play" icon
- To modify a 'Markdown'cell, double-click on it

To make nice html reports with markdown: html visualization tool 1 or html visualization tool 2, to draw nice tables, and the Ultimate guide.
Further reading on JupyterLab notebooks: Jupyter Lab documentation.