Les statistiques bayésiennes forment une trousse d'outils à garder dans votre pack sack.
En deux mots: modélisation probabiliste. Un approche de modélisation probabiliste se servant au mieux de l'information disponible. Pour calculer les probabilités d'une variable inconnu en mode bayésien, nous avons besoin:
De manière plus formelle, le théorème de Bayes (qui forme la base de l'analyse bayéseienne), dit que la distribution de probabilité des paramètres d'un modèle (par exemple, la moyenne ou une pente) est proportionnelle à la mutliplication de la distribution de probabilité estimée des paramètres et la distribution de probabilité émergeant des données.
Plus formellement,
$$P\left(\theta | y \right) = \frac{P\left(y | \theta \right) \times P\left(\theta\right)}{P\left(y \right)}$$,
où $P\left(\theta | y \right)$ $-$ la probabilité d'obtenir des paramètres $\theta$ à partir des données $y$ $-$ est la distribution de probabilité a posteriori, calculée à partir de votre a prioti $P\left(\theta\right)$ $-$ la probabilité d'obtenir des paramètres $\theta$ sans égard aux données, selon votre connaissance du phénomène $-$ et vos données observées $P\left(y | \theta \right)$ $-$ la probabilité d'obtenir les données $y$ étant donnés les paramètres $\theta$ qui régissent le phénomène. $P\left(y\right)$, la probabilité d'observer les données, est appellée la vraissemblance marginale, et assure que la somme des probabilités est nulle.
Avec la notion fréquentielle de probabilité, on teste la probabilité d'observer les données recueillies étant donnée l'absence d'effet réel (qui est l'hypothèse nulle généralement adoptée). La notion bayésienne de probabilité combine la connaissance que l'on a d'un phénomène et les données observées pour estimer la probabilité qu'il existe un effet réel. En d'autre mots, les stats fréquentielles testent si les données concordent avec un modèle du réel, tandis que les stats bayésiennes évaluent la probabilité que le modèle soit réel.
Le hic, c'est que lorsqu'on utilise les statistiques fréquentielles pour répondre à une question bayésienne, on s'expose à de mauvaises interprétations. Si l'on désire, par exemple, évaluer la probabilité de l'existance de vie sur Mars, on devra passer par le bayésien, car avec les stats fréquentielles, l'on devra plutpot conclure si les données sont conforme ou non avec l'hypothèse de la vie sur Mars (exemple tirée du blogue Dynamic Ecology).
Bien que la formule du théorème de Bayes soit plutôt simple, calculer une fonction a posteriori demandera de passer par des algorithmes de simulation, ce qui pourrait demander une bonne puissance de calcul, et des outils appropriés.
En particulier, le module générique pymc3 permet de générer une panoplie de modèles bayésiens. Le module bambi offre un interface convivial pour générer des modèles bayésiens mixtes.
Empruntons un exemple du livre Introduction to WinBUGS for Ecologists: A Bayesian Approach to Regression, ANOVA and Related Analyses, de Marc Kéry et examinons la masse de faucons pélerins. Mais alors que Marc Kéry utilise WinBUGS, un logiciel de résolution de problème en mode bayésien, nous utiliserons PyMC3.
Pour une première approche, nous allons estimer la masse moyenne d'une population de faucons pélerins.
À titre de données, générons des nombres aléatoires. Cette stratégie permet de valider les statistiques en les comparant aux paramètre que l'on impose. Ici, nous imposons une moyenne de 600 grammes et un écart-type de 30 grammes. Générons une séries de données avec 20 échantillons.
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
import scipy.stats as stats
%matplotlib inline
In [2]:
np.random.seed(63123) # fait en sorte que les mêmes nombres aléatoires sont tirés à chaque fois que la cellule est exécutée
y20 = np.random.normal(loc=600, scale=30, size=20)
y200 = np.random.normal(loc=600, scale=30, size=200)
plt.figure()
plt.subplot(121)
sns.distplot(y20)
plt.title('20 écchantillons')
plt.subplot(122)
sns.distplot(y200)
plt.title('200 écchantillons')
Out[2]:
Calculons les statistiques de manière classique.
In [3]:
def confidence_interval(x, on='devidation', distribution = 't', confidence=0.95):
m = np.mean(x)
se = np.std(x)
n = len(x)
if distribution == 't':
error = se * stats.t.ppf((1+confidence)/2., n-1)
elif distribution == 'normal':
error = se * stats.norm.ppf((1+confidence)/2.)
if on == 'error':
error = error/np.sqrt(n)
return((m-error, m, m+error))
In [4]:
print('Déviation, 95%:', np.round(confidence_interval(y20, on='deviation', confidence=0.95), 2))
print('Erreur, 95%:', np.round(confidence_interval(y20, on='error', confidence=0.95), 2))
print('Écart-type:', np.round(np.std(y20), 2))
In [5]:
print('Déviation, 95%:', np.round(confidence_interval(y20, on='deviation', distribution = 'normal', confidence=0.95), 2))
print('Erreur, 95%:', np.round(confidence_interval(y20, on='error', distribution = 'normal', confidence=0.95), 2))
print('Écart-type:', np.round(np.std(y20), 2))
En faisant cela, vous prenez pour acquis que les données sont distribuées normalement. En fait, nous savons qu'elles devraient l'être pour de grands échantillons, puisque nous avons nous-même généré les données. Par contre, comme observateur par exemple de la série de 20 données générées, la distribution est définitivement asymétrique. Sous cet angle, la moyenne, ainsi que l'écart-type, pourraient être des paramètres biaisés.
Par ailleurs, il est instructif de visualiser jusqu'à quel point un petit échantillon peut aboutir à des distributions différentes. La cellule suivante peut être exécutées plusieurs fois.
In [6]:
sns.distplot(np.random.normal(loc=600, scale=30, size=20))
Out[6]:
Nous pouvons justifier le choix d'une loi normale par des connaissances a priori des distributions de masse parmi des espèces d'oiseau. Ou bien transformer les données pour rendre leur distribution normale.
PyMC3 isole la classe Model, qui initie la définition du problème, dans un segment de code indenté sous un contexte with. À l'intérieur de ce segment indenté, on définit la connaissance a priori sous forme de variables aléatoires non-observées (ou libres) selon une distribution. Les fonctions de vraissemblance sont quant à elles définient sous forme de variables aléatoires observées.
Prenons l'exemple des faucons pélerins.
In [7]:
#!export MKL_THREADING_LAYER=GNU
import pymc3 as pm
Initions l'instance model de pm.Model() avec un a priori sur la masse moyenne. Disons que nous ne savons pas à quoi ressemble la moyenne du groupe a priori. Nous pouvons utiliser un a priori vague, où la masse moyenne peut prendre n'importe quelle valeur entre 0 et 2000 grammes, sans préférence: nous lui imposons donc un a priori selon une distribution uniforme.
In [8]:
np.random.seed(17109)
with pm.Model() as model:
mass_m = pm.Uniform('mass_m', lower=0, upper=2000)
mass_sd = pm.Uniform('mass_sd', lower=0, upper=100)
mass_pop = pm.Normal('mass_pop', mu=mass_m, sd=mass_sd) # vérifier si ça concorde avec l'analyse fréquentielle
obs = pm.Normal('obs', mu=mass_m, sd=mass_sd, observed=y20)
trace = pm.sample(1000, njobs=4)
In [9]:
pm.traceplot(trace)
Out[9]:
In [10]:
pm.plot_posterior(trace[:500], varnames=['mass_m', 'mass_sd', 'mass_pop']) # burn in [:500]
Out[10]:
Population: entre 536 et 675. Braquette plus large que celle de l'analyse fréquenntielle. (uninformative prior).
In [11]:
print('Déviation, 95%:', np.round(confidence_interval(y20, on='deviation', confidence=0.95), 2))
print('Erreur, 95%:', np.round(confidence_interval(y20, on='error', confidence=0.95), 2))
print('Écart-type:', np.round(np.std(y20), 2))
Rappelons les résultats obtenus en mode fréquentiels.
Même avec de vagues a prioris et peu de données, nous obtenons des résultats conformes à l'analyse fréquentielle. Néanmoins, les résultats peuvent être interprétés de manière différente. En ce qui a trait à la moyenne:
Fréquentiel. Il y a une probabilité de 95% que mes données aient été générées à partir d'une moyenne se situant entre 595 et 622 grammes.
Bayésien. Étant donnée mes connaissances (vagues) de la moyenne et de l'écart-type avant de procéder à l'analyse (a priori), il y a une probabilité de 95% que la moyenne de la masse de la population se situe entre 594 et 623 grammes.
In [12]:
pm.traceplot(trace)
Out[12]:
Nous avons maintenant une idée de la distribution de moyenne de la population. Mais, rarement, une analyse s'arrêtera à ce stade. Il arrive souvent que l'on doive comparer les pparamètres de deux, voire plusieurs groupes. Par exemple, comparer des populations vivants dans des écosystèmes différents, ou comparer un traitement à un placébo. Ou bien, comparer, dans une même population de faucons pélerins, l'envergure des ailes des mâles et celle des femelles.
Pour comparer des groupes, on exprime généralement une hypothèse nulle, qui typiquement pose qu'il n'y a pas de différence entre les groupes. Puis, on choisit un test statistique pour déterminer si les distributions des données observées sont plausibles dans si l'hypothèse nulle est vraie.
En d'autres mots, le test statistique exprime la probabilité que l'on obtienne les données obtenues s'il n'y avait pas de différence entre les groupes.
Par exemple, si
cela signifie qu'il y a une probabilité de 5% que vous ayiez obtenu ces données s'il n'y avait en fait pas de différence entre les groupe. Il serait donc peu probable que vos données euent été générées comme telles s'il n'y avait en fait pas de différence.
Exemple des envergures des ailes depuis Kéry (2010).
In [13]:
n_f = 30
moy_f = 105
n_m = 20
moy_m = 77.5
sd_fm = 2.75
np.random.seed(213156)
envergure_f = np.random.normal(loc=moy_f, scale=sd_fm, size=n_f)
envergure_m = np.random.normal(loc=moy_m, scale=sd_fm, size=n_m)
sns.distplot(envergure_f, label='Femelle')
sns.distplot(envergure_m, label='Mâle')
plt.legend()
Out[13]:
In [14]:
print(pd.DataFrame(envergure_f).describe())
print(pd.DataFrame(envergure_m).describe())
Évaluer s'il y a une différence significative peut se faire avec un test de t (ou de Student).
In [15]:
stats.ttest_ind(envergure_f, envergure_m)
Out[15]:
La probabilité que les données ait été générées de la sorte si les deux groupes n'était semblables est très faible. On obtiendrait sensiblement les mêmes résultats avec une régression linéaire.
In [16]:
import statsmodels.api as sm
import statsmodels.formula.api as smf
In [17]:
envergure = pd.DataFrame({'is_female': np.concatenate((np.repeat(1, n_f), np.repeat(0, n_m))),
'envergure': np.concatenate((envergure_f, envergure_m))})
linmod_envergure = smf.ols('envergure ~ is_female', data=envergure).fit()
linmod_envergure.summary()
Out[17]:
Les valeurs de t sont identiques, mais le modèle linéaire est plus informatif. Il nous apprend que l'envergure des ailes des femelles est plus en moyenne plus élevée de 28.8 cm que celle des mâles, avec un intervalle de confiance d'un niveau de 95% entre 27.4 et 30.1.
En mode bayésien,
Utilisons l'information dérivée de statistiques classiques dans nos a priori. Oui-oui, on peut faire ça. Mais attention, un a priori trop précis ou trop collé sur nos données orientera le modèle vers une solution préalablement établie: ce qui e constituerait aucune avancée par rapport à l'a priori. Nous allons utiliser a priori pour les deux groupes la moyenne des deux groupes, et comme dispersion la moyenne le double de l'écart-type. Rappelons que cet écart-type est l'a priori de écart-type sur la moyenne, non pas de la population.
Procédons à l'instanciation d'un modèle pm.Model() dans l'objet model. Nous lui ajouterons des objets par bloc. D'abord, utilisons les tendances centrales des envergures des mâles et des femelles.
In [18]:
prior_moy = envergure['envergure'].mean()
prior_moy_sd = envergure['envergure'].std()
with pm.Model() as model:
env_f_moy = pm.Normal('env_f_moy', prior_moy, sd=prior_moy_sd)
env_m_moy = pm.Normal('env_m_moy', prior_moy, sd=prior_moy_sd)
Quant à l'écart-type de la population, un a priori vague et uniforme peut être utilisé, s'étalant uniformément entre 5 et 30 cm.
In [19]:
prior_sd_ll = 1
prior_sd_ul = 10
with model:
env_f_sd = pm.Uniform('env_f_sd', lower=prior_sd_ll, upper=prior_sd_ul)
env_m_sd = pm.Uniform('env_m_sd', lower=prior_sd_ll, upper=prior_sd_ul)
Il serait possible d'utiliser une distribution normale pour comparer les moyennes, ce qui serait l'équivalent de procéder à une ANOVA. Toutefois, comparons notre modèle bayésien au test de t avec une distribution de t (ou de Student). La distribution de t comprend aussi un paramètre $\nu$ qui fait converger vers la courbe normale un échantillon dont le nombre augmente. Kruschke (2012), cité dans un tutoriel de test de t de PyMC3, suggère une loi exponentielle comme a priori pour ce paramètre. Posons 20 comme espérance de la distribution exponentielle.
In [20]:
with model:
nu_f = pm.Exponential('nu_f_minus_one', 1.0/(20-1)) + 1
nu_m = pm.Exponential('nu_m_minus_one', 1.0/(20-1)) + 1
Dans PyMC3, la distribution de t est définie non pas selon un écart-type, mais selon une précision (lam), égale à l'inverse de la variance.
In [21]:
with model:
precision_f = env_f_sd**-2
precision_m = env_m_sd**-2
female = pm.StudentT('female', nu=nu_f, mu=env_f_moy,
lam=precision_f, observed=envergure_f)
male = pm.StudentT('male', nu=nu_m, mu=env_m_moy,
lam=precision_m, observed=envergure_m)
In [22]:
with model:
diff_of_means = pm.Deterministic('difference of means', env_f_moy - env_m_moy)
diff_of_stds = pm.Deterministic('difference of stds', env_f_sd - env_m_sd)
effect_size = pm.Deterministic('effect size',
diff_of_means / np.sqrt((env_f_sd**2 + env_m_sd**2) / 2))
In [23]:
with model:
trace = pm.sample(2000, init=None, njobs=4)
In [24]:
pm.traceplot(trace)
Out[24]:
In [25]:
pm.plot_posterior(trace[500:],
varnames=['env_f_moy', 'env_m_moy', 'env_f_sd', 'env_m_sd', 'nu_f_minus_one', 'nu_m_minus_one']);
In [26]:
pm.plot_posterior(trace[1000:], varnames=['difference of means', 'difference of stds', 'effect size'], ref_val=0);
In [27]:
pm.forestplot(trace[500:], varnames=[v.name for v in model.vars])
Out[27]:
In [28]:
pm.summary(trace[500:], varnames=['difference of means', 'difference of stds', 'effect size'])
Out[28]:
In [ ]:
Prenons l'exemple l'évolution de la population de renette faux-grillon (trouver les données)
In [29]:
np.random.seed(56428)
intercept = 40
pente = -1.5
variance_residuelle = 25
annees = np.arange(1, 17, 1)
erreur = np.random.normal(loc=0, scale=np.sqrt(variance_residuelle), size=len(annees))
population = intercept + pente * annees + erreur
population_df = pd.DataFrame({'Année': annees,
'Population': population})
population_df.plot('Année', 'Population')
Out[29]:
In [30]:
import statsmodels.api as sm
import statsmodels.formula.api as smf
linmod_pop = smf.ols('Population ~ Année', data=population_df).fit()
linmod_pop.summary()
Out[30]:
In [31]:
population_df['Prédit'] = linmod_pop.predict()
population_df.plot('Année', ['Population', 'Prédit'])
Out[31]:
In [32]:
np.random.seed(10573)
with pm.Model() as model:
intercept = pm.Normal('intercept', 40, sd=5)
pente = pm.Normal('pente', 0, sd=5)
erreur_sd = pm.HalfCauchy('erreur_sd', beta=10, testval=1.0)
pop = pm.Normal('population', mu=intercept + pente*annees,
sd=erreur_sd, observed=population)
# pm.glm.GLM.from_formula('population ~ annees', population_df) # class GLM non trouvée
trace = pm.sample(2000, init=None, njobs=4)
In [33]:
pm.traceplot(trace);
Il n'y a pas de manière unique d'aborder ce genre de problème. On doit choisir la manière de représenter les données (brutes, transformées, catégorisation, etc.) les données à inclure (valeurs aberrantes, données manquantes, effets aléatoires, etc.), le modèle utilisé (linéaire, non-linéaire) ainsi que les critères de significativité (p-value <= 0.05).
Les études scientifiques basent souvent ces choix arbitrairement, par conformité avec une littérature établie (Johnson 1999) REF. Cette approche est à haut risque de mauvaise interprétation des résultats (Goodman 1999) REF.
In [ ]:
In [ ]:
In [ ]:
Les livres de Mark Kéry, bien que rédigés pour les calculs en langage R et WinBUGS, offre une approche bien structurée et applicable directement.
Le site web de pymc3 contient la documentation pour traduire en Python les codes R/WinBUGS écrits par Mark Kéry. Le livre Probabilistic Programming & Bayesian Methods for Hackers de Cameron Davidson-Pilon (biostatisticien renommé), disponible gratuitement en ligne, offre une approche plus éducative: il s'agit d'une ressource incontournable pour les statistiques bayésiennes en Python.