Générer un échantillon de Monte Carlo de taille \(50\) pour une loi normale de moyenne \(m=2\) et d’écart-type \(s=3\). Calculer les estimateurs de Monte-Carlo de la moyenne et de l’écart-type sur ce \(50\)-échantillon. Générer ensuite un échantillon de taille \(20\,000\) et calculer l’estimateurs de Monte-Carlo de la moyenne sur ce \(20\,000\)-échantillon. Que remarque-t-on ? Quel résultat théorique illustre-t-on ici ?
En reprenant l’exemple du cours, programmer un estimation de Monte Carlo du nombre \(\pi\approx 3,1416\)
Programmer une fonction roulette_coord qui prend un seul argument ncases (représentant le nombre de valeurs possible sur la roulette utilisée) valant 35 par défaut, permettant de générer les coordonnées d’un point (entre \(0\) et \(35\)). On utilisera la fonction R sample(dont on pourra consulter l’aide grâce à la commande ?sample). La fonction renvera le vecteur des 2 coordonnées ainsi générées.
En utilsant la formule de la distance entre 2 points, programmer une fonction calculant la distance à l’origine (qui a ici pour coordonnées \((\frac{ncases}{2}, \frac{ncases}{2})\)) : \(d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}\) et vérifiant si celle-ci est inférieure ou égale au rayon du cercle unité (\(R = \frac{ncases}{2}\)). Cette fonction, qu’on appellera par exemple inside_disk() prendra deux arguments : d’une part le vecteur p des deux coordonnées du point évalué, d’autre part ncases. Cela renverra une valeure booléenne (TRUE ou FALSE) indiquant si le point est effectivement à l’intérieur du disque.
Comme le ratio entre la surface du disque de rayon \(\frac{ncases}{2}\) et celle du carré de coté \(ncases\) est égal à \(\frac{\pi}{4}\), la probabilité d’échantillonner un point à l’intérieur du disque (sachant qu’il est dans le carré) est de \(\frac{\pi}{4}\). Appuyez vous sur ce résultat pour programmer une fonction calculant l’estimation de Monte-Carlo de \(pi\) à partir d’un vecteur booléen de taille \(n\) (le nombre de points échantillonnés) qui vaut TRUE si le point se trouve effectivement à l’intérieur du cerlce, et FALSE sinon.
À l’aide du code ci-dessous, représenter les données générées puis afficher l’estimation de Monte-Carlo de \(\pi\) correspondante. Faire varier npoints et commenter. Comment améliorer l’estimation ?
# Taille de la grille (résolution)
ncases <- 35
# Taille de l'échantillon de Monte Carlo
npoints <- 200
# Génération des points
pp <- matrix(NA, ncol = 2, nrow = npoints)
for (i in 1:nrow(pp)) {
pp[i, ] <- roulette_coord(ncases)
}
# Estimateur de Monte-Carlo de Pi
in_disk <- apply(X = pp, MARGIN = 1, FUN = inside_disk, ncases = ncases)
piMC(in_disk)
# Dessin on commence par initialiser un plot vide, de la bonne taille avec
# l'argument `type = 'n'`
plot(x = pp[, 1], y = pp[, 2], xlim = c(0, ncases), ylim = c(0, ncases), axes = 0,
xlab = "x", ylab = "y", type = "n")
## on gradue les axes x puis y de 0 à ncases
axis(1, at = c(0:ncases))
axis(2, at = c(0:ncases))
## on ajoute un carré encadrant le dessin
box()
## on trace (en pointillés grâce à l'argument `lty = 3`) la grille sur laquelle
## sont échantillonnés nos points
for (i in 0:ncases) {
abline(h = i, lty = 3)
abline(v = i, lty = 3)
}
## on ajoute les points échantillonés
lines(x = pp[, 1], y = pp[, 2], xlim = c(0, ncases), ylim = c(0, ncases), xlab = "x",
ylab = "y", type = "p", pch = 16)
## on ajoute la représentation du cercle
x.cercle <- seq(0, ncases, by = 0.1)
y.cercle <- sqrt((ncases/2)^2 - (x.cercle - ncases/2)^2)
lines(x.cercle, y = y.cercle + ncases/2, col = "red")
lines(x.cercle, y = -y.cercle + ncases/2, col = "red")
## enfin on colorie en rouge les points échantillonnés à l'intérieur du cercle
lines(x = pp[in_disk, 1], y = pp[in_disk, 2], xlim = c(0, ncases), ylim = c(0, ncases),
xlab = "x", ylab = "y", type = "p", pch = 16, col = "red", cex = 0.7)Programmer un algorithme d’échantillonnage selon la loi Exponentielle de paramètre \(\lambda\) grâce à la méthode de la fonction inverse et la fonction R runif. Comparer la densité de probabilité de l’échantillon obtenu avec cet algorithme par rapport à celle obtenue en utilisant la fonction R rexp.
En reprenant l’exemple historique, programmer un algorithme de Metropolis-Hastings indépendant pour estimer la loi a posteriori du paramètre \(\theta\) (la probabilité d’avoir une fille à la naissance). On prendra comme loi instrumentale la loi a priori de \(\theta\), et on utilisera dans un premier temps un a priori uniforme. On considerera les \(493\,472\) naissances observées à Paris entre 1745 et 1770, dont \(241\,945\) furent des filles.
Programmer une fonction calculant le numérateur du posterior, dont on rappel l’expression : \(p(\theta|n,S)\propto\theta^S(1-\theta)^{n-S}\) où \(S = 241\,945\) et \(n = 493\,472\) (prévoir un argument permettant de renvoyer - ou non - le logarithme).
## [1] -445522.1## [1] -354063.6Programmer l’algorithme de Metropolis-Hastings correspondant qui renvoie un vecteur de taille \(n\) échantillonné selon la loi a posteriori. On prendra soin également de renvoyer le vecteur des probabilités d’acceptation. Que se passe-t-il si l’on ne calcul pas la probabilité d’acceptation dans l’échelle log ?
Comparer la densité a posteriori obtenu avec l’algorithme M-H sur 2000 itérations à celle théorique. On prendra soin d’enlever les 500 premières itérations afin d’obtenir la convergence de la chaîne de Markov. Commenter ces résultats, notamment au vu des différentes probabilités d’acceptation calculées au cours de l’algorithme et des différentes valeurs échantillonnée pour \(\theta\).
Considérer maintenant que l’on n’observe que \(100\) naissances dont \(49\) filles, et prendre un a priori de loi \(\text{Beta}(\alpha = 3, \beta=3)\). Programmer l’algorithme correspondant et étudier à nouveau les résultats ainsi obtenus (on pourra faire \(10\,000\) du nouvel algorithme de M-H).
En utilisant les données de l’example historique et avec un a priori de loi \(\text{Beta}(\alpha = 3, \beta=3)\), programmer un algorithme de Metropolis-Hastings à marche aléatoire (on pourra prendre un pas aléatoire gaussien d’écart-type \(0,02\)). Étudier à nouveau les résultats ainsi obtenus (on pourra faire varier la taille du pas aléatoire).
BUGS & JAGSLe projet BUGS (Bayesian inference Using Gibbs Sampling) a été initié en 1989 par l’unité de Biostatistique du MRC (Medical Research Council) de l’Université de Cambridge (au Royaume-Uni) afin de proposer un logiciel flexible pour l’analyse bayésienne de modèles statistiques complexes à l’aide d’algorithmes MCMC. Son implémentation la plus connue est WinBUGS, un logiciel clic-bouton disponible sous le système d’exploitation Windows. OpenBUGS est une implémentation fonctionnant sous Windows, Mac OS ou Linux. JAGS (Just another Gibbs Sampler) est une autre implémentation plus récente qui s’appuie également sur le langage BUGS. Enfin, il faut également noter le logiciel STAN, récemment développé à l’Université de Columbia qui n’est similaire à BUGS que dans son interface, s’appuyant sur des algorithmes MCMC innovants, comme par exemple les approches de Monte-Carlo Hamiltonien ou l’approche variationnelle. Une ressource très utile est le manuel de l’utilisateur de JAGS.
Afin de se familiariser avec le logiciel JAGS (et son interface depuis R via le package rjags), nous nous intéressons ici à l’estimation a posteriori de la moyenne et de la variance d’un échantillon que l’on va modéliser avec une loi normale.
rjags.Un model en BUGS a 3 composants :
.txt) selon la syntaxe propre à BUGSGénerer \(N=50\) observations selon une loi normale de moyenne \(m=2\) et d’écart-type \(s=3\).
Lire l’aide du package rjags, puis enregistrer dans un fichier .txt le modèle en BUGS défini à l’aide du code suivant :
# Modèle
model{
# Vraisemblance
for (i in 1:N){
obs[i]~dnorm(mu,tau)
}
# A priori
mu~dnorm(0,0.0001) # propre mais très plat (faiblement informatif)
tau~dgamma(0.0001,0.0001) # propre mais très plat (faiblement informatif)
# Variables d'interet
sigma <- pow(tau, -0.5)
}Chaque fichier de spécification de modèle doit commencer par l’instruction model{ qui indique à JAGS qu’il s’apprête à recevoir la spécification d’un modèle. Ensuite il faut mettre en place le modèle en bouclant (boucle for) sur l’ensemble des données. Ici, nous voulons déclarer que qu’il y a N données, et que pour chacune d’entre elle obs[i] suit une loi normal (caractérisée avec la commande dnorm) de moyenne mu et de précision tau. Attention : dans BUGS la distribution Normale est paramétrée par sa précision, qui est l’inverse de sa variance (\(\tau = 1/\sigma^2\)). Ensuite, il s’agit de définir les lois a priori pour nos deux paramètres mu et tau, qui sont identiques pour chaque données. Pour mu, nous prenons un a priori selon la loi Normale de moyenne \(0\) et de précision \(10^{-4}\) (donc de variance \(10\,000\) : cela correspond à un a priori faiblement informatif et relativement étalé au vu de l’échelle de nos données. Pour tau nous prenons la loi a priori conjuquée dans un modèle Gaussien, c’est-à-dire la loi Gamma (avec des paramètres très faibles, là encore pour rester le moins informatif possible). Enfin, nous donnons une définition déterministe de paramètre d’intérêt supplémentaire, ici l’écart-type sigma.
À l’aide de la fonction jags.model(), créer un objet jags dans R.
À l’aide de la fonction coda.samples(), obtenir un échantillon de taille \(2\,000\) des distibutions a posteriori des paramètres de moyenne et d’écart-type
En étudiant le résultat renvoyé par la fonction coda.samples, calculer les estimateurs de la moyenne a posteriori et de la médiane a posteriori pour mu et sigma et donner un intervalles de crédibilité à 95%.
Charger le package coda. Il s’agit d’un package contenant un certain nombre de fonctions pour le diagnostic de convergence des algorithmes MCMC et le traitement des résultats.
Afin de pouvoir diagnostiquer la convergence d’un algorithme MCMC, il est nécessaire de générer plusieurs chaînes de Markov différentes, partant de valeurs initiales différentes. Recréer un objet jags en spécifiant l’utilisation de 3 chaines de Markov parallèles et en initialisant mu et tau à \((0; -10; 100)\) et \((1; 0,01; 0.1)\) respectivement. Commenter.
À l’aide de la fonction gelman.plot(), tracer la statistique de Gelman-Rubin.
À l’aide des fonctions autocorr.plot() et acfplot() évaluer l’autocorrélation des paramètres étudiés.
À l’aide de la fonction cumuplot() évaluer les quantiles cumulés des paramètres étudiés. Comment les interpréter ?
À l’aide de la fonction crosscorr.plot() évaluer les corrélations entre les paramètres étudiés. Comment les interpréter ?
À l’aide de la fonction hdi du package R HDInterval donner les intervalles de crédibilités de plus haute densité a posteriori à 95%, et les comparer à ceux obtenu à partir des quantiles \(2,5\)% et \(97,5\)%.
L’essai clinique randomisé EOLIA (Combes et al. 2018) a évalué un nouveau traitement du syndrome respiratoire aigu sévère (SDRA) en comparant le taux de mortalité après 60 jours chez 249 patients répartis de manière aléatoire parmi un groupe contrôle (traité de manière traditionnelle, c’est-à-dire par ventilation mécanique) et un groupe traitement recevant une oxygénation extracorporelle par membrane (ECMO), le nouveau traitement étudié. Une analyse fréquentiste des données a conclu à un Risque Relatif de décès de \(0,76\) dans le groupe ECMO par rapport au groupe témoin (en intention de traiter), avec \(IC_{95\%} =[0,55 ; 1,04]\) et la p-valeur associée de \(0,09\).
(Goligher et al. 2018) ont ré-analysé ces données avec une approche bayésienne, explorant différentes façon de quantifier et de résumer les résultats apportés par cette études.
Écrire le modèle Bayésien utilisé par (Goligher et al. 2018).
Ecrire le modèle BUGS correspondant.
| Contrôle | ECMO | |
|---|---|---|
| \(n\) observé | 125 | 57 |
| nombre de décédés à 60 jours | 124 | 44 |
Commencer par créer 2 vecteurs de données binaires (1 ou 0) ycontrol et yecmo encodant les observations individuelles résumées dans la table ci-dessus. Utiliser ensuite les fonctions jags.model() et coda.samples() pour reproduire les estimations de (Goligher et al. 2018) (ProTip: utiliser la fonction window() afin de retirer les observations correspondant à la phase de chauffe de coda.samples fonction).
Évaluer la convergence, puis commenter les résultats et les estimations (ProTip: la fonction R effectiveSize() permet de calculer la taille d’echantillon effective).
Changer pour une loi a priori plus informative en utilisant une distribution Gaussienne pour le \(log(RR)\), centrée en log(0.78) et avec un écart-type de 0.15 (dans l’échelle du log(RR), ce qui correspond à une précision \(\approx 45\)). Commenter les resultats. Essayer d’autres distributions a priori.
En 2014, Crins et al. ont publiés une méta-analyse pour évaluer l’incidence des rejets aigu (RA) avec ou sans traitement avec recepteur antagonistes à l’Interleukin-2. Dans cet exercice nous allons recréer cette analyse.
Charger le package R bayesmeta (Röver 2017) et les données de (Crins et al. 2014) avec la commande R: data("CrinsEtAl2014").
library(bayesmeta)
data(CrinsEtAl2014)Explorer l’application web shiny accompagnat le package et disponbible à l’adresse : http://ams.med.uni-goettingen.de:3838/bayesmeta/app/. Commenter les sorties et resultats, essayer différentes loi a priori, etc.
Maintenant depuis R, avec la fonction escalc() du package metafor, calculer les log odds ratios estimés à partir des 6 études considérées avec les variances d’échantillonnage associées (ProTip: lire la section Measures for Dichotomous Variables de l’aide de la fonction escalc()). Vérifier la cohérence avec l’application shinyen ligne (ProTip: sigma est l’écart-type, i.e. la racine carré de la variance d’échantillonnage vi)
library("metafor")
crins.es <- escalc(measure = "OR", ai = exp.AR.events, n1i = exp.total, ci = cont.AR.events,
n2i = cont.total, slab = publication, data = CrinsEtAl2014)
crins.es[, c("publication", "yi", "vi")]| publication | yi | vi |
|---|---|---|
| Heffron (2003) | -2.3097026 | 0.3593718 |
| Gibelli (2004) | -0.4595323 | 0.3095760 |
| Schuller (2005) | -2.3025851 | 0.7750000 |
| Ganschow (2005) | -1.7578579 | 0.2078161 |
| Spada (2006) | -1.2584610 | 0.4121591 |
| Gras (2008) | -2.4178959 | 2.3372623 |
Les log odds ratios (OR) sont symmétriques autour de zéros et ont leur distribution d’échantillonnage est plus proche d’une distribution Gaussienne que dans l’échelle narturelle des OR. C’est pourquoi ils sont généralement préférés dans le cadre des méta-analyses. Leur variances d’échantillonnage sont calculées comme la somme de l’inverse de tous les comptes dans le tableau de contingence \(2 \times 2\) associé (Fleiss and Berlin 2009).
Effectuer une méta-analyse avec un effet-aléatoire sur ces données avec la fonction bayesmeta() du package R bayesmeta, depuis R. Utiliser une loi a priori uniforme sur \([0,4]\) pour \(\tau\) (la précision) et une loi a priori Gaussienne centrée autour de \(0\) et avec un écart-type de 4 pour \(\mu\) la moyenne.
Ecrire le modèle Bayésien correspondant à cette méta-analyse à effet aléatoire.
Utiliser rjags pour estimer ce même modèle. Commenter.
Nous allons commencer par une lecture critique de l’article de (Kaguelidou et al. 2016).
Lister les éléments méthodologiques qui manquent dans cet article.
Lire et discuter la Table 3.
Charger le package R bcrm, et mener un essai CRM imaginaire en utilisant les commandes suivantes :
library(bcrm)
sdose <- c(1, 1.5, 2, 2.5, 3)
dose.label <- c(5, 10, 15, 25, 40, 50, 60)
bcrm(stop = list(nmax = 42), p.tox0 = p.tox0, dose = dose.label, ff = "power", prior.alpha = list(1,
1, 1), target.tox = 0.3, start = 1)On veut comparer deux moyennes.
Commencer par générer 2 échantillons chacun de taille \(n=15\), de même variance \(s^2 = 9\) et de moyenne respective \(10\) et \(8\). On les notera échantillons a et b.
Rappeler les hypothèses du t-test de Student, et effectuer un test unilatérale de l’hypothèse alternative unilatérale \(H_1 : \mu_a >\mu_b\) au niveau de risque de l’erreur de type I \(\alpha=5%\)
Calculer la borne inférieure de l’intervalle de confiance bootstrap de la différence de moyenne à partir de \(20\;000\) échantillons bootstrap. Calculer également la p-valeur associée au test unilatéral. Enfin, représenter la distribution des différences bootstrap, et faire figurer la valeur observée.
Grâce au bootstrap, estimer l’écart-type de la différence de moyenne, ainsi que d’autres quantités d’intérêt telles que la borne inférieure de l’intervalle de confiance et la p-valeur. Commenter (notamment en comparantl’estimation de l’écart-type au résultat de la formule asymptotique \(\sqrt{\frac{Var(a)}{n} + \frac{Var(b)}{n}}\)).
Recommencer avec \(n=200\) (au lieu de \(n=15\)) et commenter à nouveau.
Utiliser maintenant la méthode de Monte-Carlo pour estimer la variabilité de la différence de moyenne étudiée (on appelle également cette technique le bootstrap paramétrique). Puis utiliser à nouveau Monte-Carlo pour calculer la p-valeur associée. Comparer aux résultats précédents et commenter.
À l’aide de la méthode de Monte-Carlo verifier que le t-test de Student contrôle bien l’erreur de type I au niveau \(\alpha = 5%\) voulu lorsque \(n=15\)?
À l’aide de la méthode de Monte-Carlo , représenter la puissance du t-test de Student pour différentes valeur de \(n\) ?
BONUS : utiliser un test par permutation, et comparer aux résultats précédents. On pourra lire l’article de (Phipson and Smyth 2010) à propos du calcul de la p-valeur pour un test par permutation.
Combes, Alain, David Hajage, Gilles Capellier, Alexandre Demoule, Sylvain Lavoué, Christophe Guervilly, Daniel Da Silva, et al. 2018. “Extracorporeal Membrane Oxygenation for Severe Acute Respiratory Distress Syndrome.” New England Journal of Medicine 378 (21): 1965–75. https://doi.org/10.1056/NEJMoa1800385.
Crins, Nicola D., Christian Röver, Armin D. Goralczyk, and Tim Friede. 2014. “Interleukin-2 Receptor Antagonists for Pediatric Liver Transplant Recipients: A Systematic Review and Meta-Analysis of Controlled Studies.” Pediatric Transplantation 18 (8): 839–50. https://doi.org/10.1111/petr.12362.
Fleiss, Joseph L., and Jesse A. Berlin. 2009. “Effect Sizes for Dichotomous Data.” In The Handbook of Research Synthesis and Meta-Analysis, 2nd Ed, 237–53. New York, NY, US: Russell Sage Foundation.
Goligher, Ewan C., George Tomlinson, David Hajage, Duminda N. Wijeysundera, Eddy Fan, Peter Jüni, Daniel Brodie, Arthur S. Slutsky, and Alain Combes. 2018. “Extracorporeal Membrane Oxygenation for Severe Acute Respiratory Distress Syndrome and Posterior Probability of Mortality Benefit in a Post Hoc Bayesian Analysis of a Randomized Clinical Trial.” JAMA 320 (21): 2251. https://doi.org/10.1001/jama.2018.14276.
Kaguelidou, Florentia, Corinne Alberti, Valerie Biran, Olivier Bourdon, Caroline Farnoux, Sarah Zohar, and Evelyne Jacqz-Aigrain. 2016. “Dose-Finding Study of Omeprazole on Gastric pH in Neonates with Gastro-Esophageal Acid Reflux Using a Bayesian Sequential Approach.” Edited by Imti Choonara. PLOS ONE 11 (12): e0166207. https://doi.org/10.1371/journal.pone.0166207.
Phipson, Belinda, and Gordon K Smyth. 2010. “Permutation P-Values Should Never Be Zero: Calculating Exact P-Values When Permutations Are Randomly Drawn.” Statistical Applications in Genetics and Molecular Biology 9 (1): 1544–6115. https://doi.org/10.2202/1544-6115.1585.
Röver, Christian. 2017. “Bayesian Random-Effects Meta-Analysis Using the Bayesmeta R Package.” arXiv Preprint 1711.08683. http://www.arxiv.org/abs/1711.08683.