Code R et WinBUGS : différentes solutions
La solution à l'exemple de l'essai thérapeutique précédent : dans R ou dans WinBUGS, principe général
simulations de valeurs issues des lois Beta() respectives des deux groupes
itérations MCMC ou échantillonnage aléatoire dans la loi
calcul des indices de taille d'effet sur chaque itération de la ou des chaine(s)
test de ces indices par rapport à une valeur de référence
Exemple :
échantillonnage dans les lois des proportions
M<-100000
p1<-rbeta(M,110+1,90+1)
p2<-rbeta(M,80+1,110+1)
round(p1[1:10],3)
round(p2[1:10],3)
# à titre informatif : les 10 premières valeurs des 2 vecteurs
> round(p1[1:10],3)
[1] 0.570 0.611 0.552 0.579 0.530 0.586 0.554 0.546 0.549 0.549
> round(p2[1:10],3)
[1] 0.406 0.438 0.476 0.457 0.406 0.430 0.398 0.438 0.463 0.357
Puis calcul des mesures de taille d'effet :
DR <- p1-p2
RR <- p1/p2
RR<- (p1*(1-p2))/(p2*(1-p1))
NNT<-1/DR
Le principe : R calcule, pour chaque valeur des vecteurs les valeurs des indices et en fait un vecteur final.
Donc, DR = (DR[1],DR[2], etc ..., DR[M])
et DR[i] = p1[i]-p2[i], etc idem pour tous les indices.
Complément : En pratique :
> round(quantile(DR,probs=c(0.025,0.05,0.25,0.50,0.75,0.95,0.975)),3)
2.5% 5% 25% 50% 75% 95% 97.5%
0.029 0.045 0.095 0.128 0.161 0.209 0.225
> round(quantile(RR,probs=c(0.025,0.05,0.25,0.50,0.75,0.95,0.975)),3)
2.5% 5% 25% 50% 75% 95% 97.5%
1.124 1.200 1.463 1.678 1.923 2.345 2.505
> round(quantile(OR,probs=c(0.025,0.05,0.25,0.50,0.75,0.95,0.975)),3)
2.5% 5% 25% 50% 75% 95% 97.5%
1.124 1.200 1.463 1.678 1.923 2.345 2.505
> round(quantile(NNT,probs=c(0.025,0.05,0.25,0.50,0.75,0.95,0.975)),3)
2.5% 5% 25% 50% 75% 95% 97.5%
4.347 4.717 6.158 7.757 10.479 20.893 29.437
En outre, possibilité de faire des tests :
> sum(DR>0)/M
[1] 0.99473
> sum(DR>0.15)/M
[1] 0.3275
> sum(RR>1)/M
[1] 0.99473
> sum(RR>2)/M
[1] 7e-05
> sum(RR>1.4)/M
[1] 0.24961
La même chose mais en WinBUGS : variante (1)
model{
for (i in 1:n1)
{ TA[i] ~ dbern(p1) }
p1 ~ dbeta(1,1);
for (i in 1:n2)
{ TB[i] ~ dbern(p2) }
p2 ~ dbeta(1,1);
propdiff <- p1-p2
rr <- p1/p2
or<- (p1*(1-p2))/(p2*(1-p1))
}
# Data
list(TA = c(0,1, ...), TB = c(0,1,...),n1 = 200, n2 = 200)
Le code dans WinBUGS : variante (2)
le principe de l'estimation de p est identique dans les deux groupes
[\(\to\)] en tirer partie pour simplifier le code
avec dans le code une boucle sur le \(i\) de \(p[i]\)
et calcul des indices
en utilisant le fait que p est le paramètre d'une loi binomiale
Donc...
model{
for(i in 1:2){ ## boucle
r[i] ~ dbin(p[i],n[i]) ## p : à estimer
}
## priors
p[1] ~ dbeta(1,1)
p[2] ~ dbeta(1,1)
## compute quantities of interest
DR <- p[1] - p[2] ## difference de probs
DRsup0 <- step(DR) ## diff > 0 ?
}
# Données
list(r=c(110,85), n=c(200,200))
# Valeurs initiales
list(p=c(0.5,0.5))
La cas d'une table sans aucune marge fixée
une table de dimension \(I \times J\)
soit \( I\times J\) cases
une loi multinomiale, généralisation de la loi binomiale
\(m\) valeurs possibles et effectifs total \(N\)
la densité de probabilité de la loi multinomiale :
\(\Pr(N_1= n_1, \dots, N_m=n_m) = \frac{N !}{n_1 ! \dots n_m !}\times p_{1}^{n_1} \dots p_{m}^{n_m}\)
La loi conjuguée est une loi de Dirichlet
Exemple d'une table \(2 \times 2\) :
la loi de Dirichlet a priori : \(Di(\alpha_1,\alpha_2,\alpha_3,\alpha_4)\)
les effectifs de la table : \(n_1, n_2, n_3, n_4\)
la loi de Dirichlet a posteriori : \(Di(\alpha_1+n_1,\alpha_2+n_2,\alpha_3+n_3,\alpha_4+n_4)\)
inférence --> simulations directement dans la loi a posteriori
dans R, package MCMCpack :
rdirichlet(M,paramètres)
Pour une table sans aucune marge fixée :
## ouvrir le package ci-dessous :
library(MCMCpack)
# choisir le nombre d'échantillons dans la loi = nb d'itérations
M<-100000
tablo<-rdirichlet(M, c(36,56,88,64) )
# accessoirement :
rowSums(tablo[1:10,])
OR<-tablo[,1]*tablo[,4]/(tablo[,2]*tablo[,3])
round(quantile(OR,probs=c(0.025,0.05,0.25,0.50,0.75,0.95,0.975)),3)
2.5% 5% 25% 50% 75% 95% 97.5%
0.271 0.296 0.387 0.465 0.558 0.725 0.789