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.
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).
data = np.loadtxt("Supra_data_s.txt",skiprows=1)
header = ["SEXE","AGE","taille","poids","tabac","SPORT","mesure","alcool"]
print(data[0:10,:])
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.
bmi = data[:,3:4]/data[:,2:3]**2
data = np.concatenate((data,bmi),axis=1)
header = header + ["bmi"]
print(data.shape)
Tests unitaires¶
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.
imt_sample = data[:,6]
print(stats.ttest_1samp(imt_sample,0.5))
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.
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 :
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.
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$
F = np.var(imt_men)/np.var(imt_women)
p_value = stats.f.cdf(F, len(men)-1, len(women)-1)
print(p_value)
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 :
print(stats.ttest_ind(imt_men,imt_women))
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
bp = sns.boxplot(x=data[:,4], y=data[:,6]).set(xlabel='Tabac',ylabel='IMT')
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)
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.
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)
QUESTION Creer une table de contingence qui permette de croiser la variable mesure_bin et la variable tabac
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
QUESTION Faire un test du $\chi^2$ permettant de tester le lien entre le tabagisme et l'IMT binarisé.
chi2, pvalue, degrees, expected = stats.chi2_contingency(poids_tabagisme)
print(pvalue)
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
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 corrections,sidak: one-step correctionhs,holm-sidak: step down method using Sidak adjustmentsh,holm: step-down method using Bonferroni adjustmentssh,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.
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.
m = 8
p = [0]*m
for i in range(m):
t, p[i] = stats.ttest_ind(data_men[:,i],data_women[:,i])
print(p)
QUESTION Pour un seuil $\alpha=0.05$, indiquer quels tests sont significatifs.
alpha = 0.05
test_signif = [i for i,x in enumerate(p) if x<=alpha]
print(test_signif)
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.
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)
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.
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)
QUESTION Comparer les deux methodes et conclure.