Tests statistiques et multiplicité

Python works with packages. A large number of packages can be dowloaded all together in the Anaconda environment. Start downloading here (version >3) https://www.anaconda.com/download/ Anaconda comes with

  • Spyder, a GUI for python
  • Jupyter Notebook, how to create great documents

Enjoy!

Bee careful :

  • Python is sensible to spaces in the beginning of lines,
  • Python is sensible to number type : $3$ and $3.$ are totally different: integer and float types.
In [46]:
import warnings
warnings.filterwarnings('ignore')

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd

Les données

Le fichier Supra_data_s.txt contient les données relevées lors d’une étude prospective (SUPRA - Supervision of the PRogression of Atherosclerosis in HIV-1 infected patients) parmi des participants de la cohorte Aquitaine. L’étude portait sur l’évaluation des facteurs de risques associés aux maladies cardiovasculaires chez les patients infectés par le VIH, en utilisant le marqueur IMT - intima media thickness of the carotid arteries (colonne mesure). L’étude comportait le relevé d’autres marqueurs connus comme étant des facteurs de risques des maladies cardiovasculaires : sexe (SEXE = 1 si homme, 2 si femme), age (AGE en années), taille (en cm), poids (en kg), habitudes liées au sport (0=non sportif, 1=sportif), tabac et prise d’alcool (0=consommation nulle, 1=consommation occasionnelle, 2=consommation régulière).

In [24]:
data = np.loadtxt("Supra_data_s.txt",skiprows=1)
header = ["SEXE","AGE","taille","poids","tabac","SPORT","mesure","alcool"]
print(data[0:10,:])
[[  1.    33.   170.    70.     1.     0.     0.52   1.  ]
 [  2.    33.   177.    67.     2.     0.     0.42   1.  ]
 [  2.    53.   164.    63.     1.     0.     0.65   0.  ]
 [  2.    42.   169.    76.     1.     1.     0.48   1.  ]
 [  2.    53.   152.    54.     0.     0.     0.45   1.  ]
 [  2.    50.   162.    53.     2.     0.     0.49   1.  ]
 [  1.    22.   181.    72.     0.     1.     0.42   1.  ]
 [  2.    26.   162.    61.     1.     1.     0.45   1.  ]
 [  1.    50.   170.    80.     0.     0.     0.65   2.  ]
 [  2.    47.   178.    73.     0.     0.     0.52   1.  ]]

Creation du BMI

QUESTION Créer la variable BMI Indice de Masse Corporelle telle que $$ BMI=\frac{Poids}{Taille^2}, $$ et l'ajouter au dataset.

In [25]:
bmi = data[:,3:4]/data[:,2:3]**2
data = np.concatenate((data,bmi),axis=1)
header = header + ["bmi"]
print(data.shape)
(110, 9)

Tests unitaires

In [47]:
from scipy import stats

Test de student

Comparaison à une moyenne connue

QUESTION On suppose que l’IMT moyen dans la population générale est de $0.5$. L’échantillon de l’étude est-il représentatif de la population générale en terme d’IMT sous un seuil $\alpha=0.05$ ?

Pour ça on utilisera le test de Student en comparaison à une moyenne connue.

scipy.stats.ttest_1samp permet de faire ça.

In [27]:
imt_sample = data[:,6]
print(stats.ttest_1samp(imt_sample,0.5))
Ttest_1sampResult(statistic=3.6639873759501076, pvalue=0.00038515402414989104)

Pour un seuil $\alpha=0.05$, il nous est nécessaire de rejetter l'hypothèse d'appartenance de la population d'étude à la population générale.

Visualisation de données, la boite à moustache ou boxplot

QUESTION Réaliser la boîte à moustache de l’IMT en fonction de la variable sexe.

In [28]:
import seaborn as sns

bp = sns.boxplot(x=data[:,0], y=data[:,6]).set(xlabel='Sexe',ylabel='IMT')

On peut voir que les distributions semblent être les mêmes.

QUESTION Utiliser un test de comparaison de moyenne de Student à $2$ échantillons indépendants pour s'en assurer, ou pas.

Repérer les positions des hommes et des femmes :

In [29]:
men = [i for i,x in enumerate(data[:,0]) if x==1.]
women = [i for i,x in enumerate(data[:,0]) if x==2.]

Ce qui permet de sélectionner les données d'IMT pour les femmes d'un côté et les hommes de l'autre.

In [30]:
imt_men = data[men,6]
imt_women = data[women,6]
  • Première étape : Test de l'hypothèse de l'égalité des variances :

    On utilise le F-test et l'on compare à $\alpha=0.05$

In [31]:
F = np.var(imt_men)/np.var(imt_women)
p_value = stats.f.cdf(F, len(men)-1, len(women)-1)
print(p_value)
0.34570300779740415

On a une p-value supérieure à 0.5 donc on suppose que l'hypothèse d'égalité des variances est à conserver.

Nous pouvons donc passer à la deuxième étape.

  • Deuxième étape : Test de Student à $2$ échantillons :
In [32]:
print(stats.ttest_ind(imt_men,imt_women)) 
Ttest_indResult(statistic=-0.20707837588744438, pvalue=0.8363385826501982)

La p-value est supérieure à $\alpha=0.05$ donc il faut conserver l'hypotèse $\mathcal{H}_0$ : Les deux populations ont le même IMT moyen.

INTERPRETATION

QUESTION Tracer un Boxplot de l'IMT en fonction de la variable tabac

In [33]:
bp = sns.boxplot(x=data[:,4], y=data[:,6]).set(xlabel='Tabac',ylabel='IMT')
In [34]:
tabac_0 = data[[i for i,x in enumerate(data[:,4]) if x==0.],6]
tabac_1 = data[[i for i,x in enumerate(data[:,4]) if x==1.],6]
tabac_2 = data[[i for i,x in enumerate(data[:,4]) if x==2.],6]

F, p = stats.f_oneway(tabac_0,tabac_1,tabac_2)

print(p)
0.08181858704650491

La p-value est supérieure à $\alpha=0.05$ donc il faut conserver l'hypotèse $\mathcal{H}_0$ : Les trois populations ont le même IMT moyen.

INTERPRETATION

QUESTION Binariser la variable mesure en mesure_bin qui vaut $0$ si mesure est inférieure à $0.5$ et $1$ sinon.

In [35]:
mesure_bin = [0]*len(data[:,0])
mesure_1_pos = [i for i,x in enumerate(data[:,6]) if x>=0.5]
for i in mesure_1_pos:
    mesure_bin[i] = 1
print(mesure_bin)
[1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0]

QUESTION Creer une table de contingence qui permette de croiser la variable mesure_bin et la variable tabac

In [36]:
data_contin = np.array([data[:,4],mesure_bin]).T
pd_data = pd.DataFrame(data=data_contin,columns=["Tabagisme (X)","IMT (Y)"])
poids_tabagisme = pd.crosstab(index=pd_data["Tabagisme (X)"], 
                           columns=pd_data["IMT (Y)"])
poids_tabagisme
Out[36]:
IMT (Y) 0.0 1.0
Tabagisme (X)
0.0 41 31
1.0 5 13
2.0 5 15

QUESTION Faire un test du $\chi^2$ permettant de tester le lien entre le tabagisme et l'IMT binarisé.

In [37]:
chi2, pvalue, degrees, expected = stats.chi2_contingency(poids_tabagisme)
print(pvalue)
0.00903901150181748

La p-value est inférieure à $\alpha=0.05$ donc il faut rejetter l'hypotèse $\mathcal{H}_0$.

INTERPRETATION

Multiplicité des tests

Afin d'ajuster les tests, on chargera la librairie statsmodels.stats.multitest que l'on renomera smm

In [45]:
import statsmodels.stats.multitest as smm

La fonction alors utilisée est multipletests.

Où l'argument method peut être l'un des suivants :

  • b, bonferroni : one-step correction
  • s, sidak : one-step correction
  • hs, holm-sidak : step down method using Sidak adjustments
  • h, holm : step-down method using Bonferroni adjustments
  • sh, simes-hochberg : step-up method (independent)
  • hommel : closed method based on Simes tests (non-negative)
  • fdr_i, fdr_bh : Benjamini/Hochberg (non-negative)
  • fdr_n, fdr_by : Benjamini/Yekutieli (negative)
  • fdr_tsbh : two stage fdr correction (Benjamini/Hochberg)
  • fdr_tsbky : two stage fdr correction (Benjamini/Krieger/Yekutieli)
  • fdr_gbs : adaptive step-down fdr correction (Gavrilov, Benjamini, Sarkar)

    On ne va s'intéresser qu'à certaines de ces méthodes. On note par la suite :

  • $\alpha$ le niveau des tests,
  • $m$ le nombre de tests

Un test multiple

On considère maintenant les $m=8$ t-tests correspondant aux $8$ variables considérées précédemment, autres que le sexe.

QUESTION Créer deux datasets data_women et data_men qui contiennent les restrctions du dataset initial aux femmes et aux hommes respectivement.

In [58]:
men = [i for i,x in enumerate(data[:,0]) if x==1.]
women = [i for i,x in enumerate(data[:,0]) if x==2.]

data_men = data[men,1:9]
data_women = data[women,1:9]

QUESTION Créer le vecteur p qui contiennent les résultats des $m=8$ t-tests unitaires.

In [62]:
m = 8
p = [0]*m
for i in range(m):
   t, p[i] = stats.ttest_ind(data_men[:,i],data_women[:,i])
print(p)
[0.00962906732269608, 3.3451449978972865e-14, 8.333867926819101e-09, 0.22164597588141535, 0.3632450264859446, 0.8363385826501982, 0.006654207454445699, 1.9180208767383257e-09]

QUESTION Pour un seuil $\alpha=0.05$, indiquer quels tests sont significatifs.

In [63]:
alpha = 0.05
test_signif = [i for i,x in enumerate(p) if x<=alpha]

print(test_signif)
[0, 1, 2, 6, 7]

Une méthode contrôlant le FWER : Bonferroni

Les m tests sont effectués au niveau $\alpha^\star= \frac{\alpha}{m}$ pour contrôler le FWER au niveau $\alpha$.

La méthode est method=b ou bien method=bonferroi.

QUESTION Estimer les p-values corrigées dans le cadre de la correction de Bonferroni observer l'effet sur les pvalues ajustées et concllure.

In [67]:
alpha = 0.05
rej, pval_corr = smm.multipletests(p, alpha=alpha, method="b")[:2]

test_signif_bonferroni = [i for i,x in enumerate(pval_corr) if x<=alpha]
print(test_signif_bonferroni)
[1, 2, 7]

Une méthode contrôlant le FDR : Benjamini and Hochberg (1995)

La méthode s'applique de cette manière :

  • Ordonner les $m$ p-values par ordre croissant de sorte que $$ p_1\leq p_2\leq \cdots \leq p_m $$
  • Trouver $L$ tel que $$ L=\max\{j\in\{1,\cdots,m\}|p_j<\alpha\frac{j}{m}\} $$
  • Rejeter toutes les hypothèses $\mathcal{H}_{0,j}$ pour lesquelles $p_j<p_L$.

$p_L$ est appelé de seuil de rejet de BH (BH threshold).

L'argument dans la fonction est method=fdr_bh.

QUESTION Estimer les p-values corrigées dans le cadre de la correction de Benjamini-Hochberg observer l'effet sur les pvalues ajustées et concllure.

In [69]:
alpha = 0.05
rej, pval_corr_bh = smm.multipletests(p, alpha=alpha, method="fdr_bh")[:2]

test_signif_bonferroni = [i for i,x in enumerate(pval_corr_bh) if x<=alpha]
print(pval_corr_bh)
print(test_signif_bonferroni)
[1.54065077e-02 2.67611600e-13 2.22236478e-08 2.95527968e-01
 4.15137173e-01 8.36338583e-01 1.33084149e-02 7.67208351e-09]
[0, 1, 2, 6, 7]

QUESTION Comparer les deux methodes et conclure.