projet_final.Rnw

\documentclass[12pt]{article}

%ackages
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[a4paper,left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}
\usepackage[french]{babel}
\usepackage{graphicx}
\pagestyle{plain}
\usepackage{fancyhdr}
\usepackage{fancybox}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{empheq}
\usepackage{float}
\usepackage{mathrsfs}
\usepackage{bigints}

\usepackage{hyphenat}
\usepackage{listings}
\usepackage{booktabs}


\pagestyle{fancy}
\fancyhead[C]{\rule{.5\textwidth}{4\baselineskip}}% Add something BIG in the header
\setlength{\headheight}{63pt}
\linespread{1.2}
\setlength{\parindent}{0cm}
\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex}
\newcommand{\hsp}{\hspace{20pt}}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\usepackage{fourier} 
\usepackage{array}
\usepackage{makecell}


\begin{document}
\SweaveOpts{concordance=TRUE}
<<options, echo = FALSE>>=
options(prompt = " ", continue = " ", width = 85)
@

\begin{titlepage}

  \begin{center}

    %\includegraphics[scale=0.5]{logoInsa.png}~\\[1cm]
    
    \vfill
    \textsc{ \Large Compte rendu de projet }\\[1cm]
    %\vfill
    \HRule \\[0.4cm]
    {\LARGE \textbf{Design d'essai clinique en méta-analyse : risque cardiovasculaire et traitement anti-diabétique}
\\ [0.2cm] }

 \HRule \\[1cm]
\vfill

%trouver une image à mettre    \includegraphics[scale=0.2]{pacman.jpg}

%\vfill
%    \textsc{"Not all things that you learn are taught to you, but many things that you learn you realize you have taught yourself.
%”\\C. JoyBell}\\[1cm]
    
%\vfill
   
    \textsc{CHASSAGNOL Bastien, DE NAILLY Paul, PISTER Alexis, PRIN Amaury}\\[1cm]
\vfill 

    \begin{minipage}{0.4\textwidth}
      \begin{flushleft} 
        Bio-informatique et modélisation\\ 
        Année scolaire : 2018-2019\\
      \end{flushleft}
    \end{minipage}
    \begin{minipage}{0.4\textwidth}
      \begin{flushright} \large
        Encadrant : \\Mr \textsc{Fabien SUBTIL}\\
      \end{flushright}
    \end{minipage}

  \end{center}
  
\end{titlepage}


\section{Contexte biologique}
L'objectif de ce projet est de planifier une étude clinique dans le but de tester, par une étude de non-supériorité, un nouveau traitement anti-diabétique par rapport à un traitement de référence. Afin d'estimer la puissance nécessaire pour réaliser cette étude, nous allons utiliser l'approche bayésienne à partir d'études anciennes qui ont porté sur le traitement de référence. À travers ce projet, l'objectif est d'évaluer l'influence de l'approche bayésienne dans une étude clinique, en particulier sur la manière d'ajouter l'information {\it apriori} dans l'étude. Pour ce projet, nous essaierons de résumer l'information de 5 études précédentes en essayant de modéliser le délai avant un évènement cardiovasculaire. Nous simulerons également des essais cliniques, de diffentes manières, en intégrant ou non l'information {\it apriori}.

\section{Analyse des études antérieures}
Dans un premier temps, nous avons construit un modèle sous \texttt{Jags} pour estimer le temps moyen d'apparition d'un événement cardiovasculaire (noté ici $\theta$) à partir des 5 études antérieures. On suppose que la moyenne de ce paramètre est globalement la même pour toutes les études (celles-ci ayant recouru au même type de traitement de référence), mais en lui laissant pour chaque étude la possibilité de s'écarter un peu de la moyenne. En effet, selon la prise en charge des patients par l'hôpital, le profil des patients sélectionnés par chaque étude ou encore le fait que chaque étude ne consitue qu'un échantillon de la population globale, on est obligé d'introduire de la variabilité pour chaque expérience. Celle-ci est introduite avec le paramètre $ \sigma$. On a donc recouru à certaines hypothèses simplificatrices pour la construction de notre modèle, en supposant notamment que chaque patient au sein d'une même étude réagissait de la même façon au traitement (pas de variabilité intra groupe). 

 \vspace{2\baselineskip}
Concrètement, pour chaque patient de chaque étude, nous allons modéliser le temps d'apparition d'un évènement cardiovasculaire grâce à une loi exponentielle. On se sert souvent de cette loi pour déterminer le temps d'apparition d'une maladie puisque le paramètre de cette loi est $\theta$ et son espérance est $ \dfrac{1}{\theta}$, correspondant au temps moyen de l'apparition d'un événement cardiovasculaire. De plus, la loi exponetielle est dite \textbf {sans mémoire}, ce qui signiefie que la probabilité qu'un phénomène se produise après s + t heures sachant qu'il a déjà duré t heures est la même que la probabilité de durer s heures à partir du début du traitement. Elle implique donc une autre hypothèse simplificatrice en supposant que l'état des malades varie de façon constante, et ne dépend pas de ce qui a précédé. Si elle permet de simplifier amplement l'approche du problème, elle peut aussi être constestée, l'état du patient dépendant forcément des antécédents. Une autre hypothèse apportée est l'indépendance de la date d'apparition des événements entre les patients d'une même étude. Celle-ci est normalement vérifiée puisque l'état d'un patient n'influe normalement pas l'état des autres patients, d'autant plus que les maladies cardiovasculaires ne sont pas contagieuses. Ces différentes hypothèses sont indispensables pour la construction de notre modèle, dans lequel on suppose que chaque patient apporte la même quantité d'informations, et constitue une réalisation indépendante et identiquement distribuée de notre loi exponentielle. Une autre difficulté vient du fait que nous ne suivons pas sur une vie entière les patients, puisqu'assez logiquement le temps d'étude est fixé pour une certaine durée. Pour prendre en compte ces individus censurés qui sont majoritaires dans notre étude (96\% des patients), on utilise la fonction jags \texttt{dinterval} qui vaut 0 pour une donnée non censurée et 1 sinon. On cherche ensuite à comparer ces valeurs avec celles estimées par notre modèle. Comme nous ne disposons pas du temps d'étude sur chaque patient, on suppose volontairement ce temps de censure très grand (ici fixé dans notre modèle à 10 000, mais il suffit de le prendre supérieure à la valeur maximale du temps d'étude).Comme le vrai temps d'apparition n'est pas connu, pour tous les patients censurés, on indique à jags qu'on ne connait pas ce temps. Ce modèle peut être représenté par le graphe acyclique suivant:  


\begin{figure}[H]

\begin{center}
\includegraphics{dag_modele1_non_optimise.png} 
\end{center}

\end{figure}

et sera codé de cette façon dans jags: 
<<modele_non_optimise>>=
#Lecture données et variables
require(rjags)
data = read.table("dataprior_projet.txt", header=T)
event=data$event
time=data$time
etude=data$etude
vec = c("a","b","c","d","e")
etude2 = match(etude,vec)
censure = 1-event 

#calcul Ti =temps apparition d'un évènement cardiovasculaire et C temps de censure observé
C = vector(le = length(event))
C[event==0] = time[event==0]
C[event==1] = 3
Ti = vector(le = length(event))
Ti[event==0] = NA
Ti[event==1] = time[event==1]

#ecriture du modele sous jags
ms = "
model{
  for(i in 1:5){
  log_teta[i] ~ dnorm(mu,tau)
  teta[i]=exp(log_teta[i])
  }

  for (j in 1:Np){
  censure[j] ~ dinterval(time[j], C[j])
  time[j] ~ dexp(teta[etude[j]]) 
}
#lois uniformes non informatives
#ici on tire le temps moyen d'apparition selon une échelle logarithmique, 
#puisque la loi utilisant ce paramètre est exponentielle. Pour avoir donc
#une représentation la moins informative possible de ce paramètre, 
#on utilise une loi uniforme sur son logarithme.
  tau = 1/(sigma**2)
  sigma ~ dunif(0,20)
  mu ~ dunif(-20,10)
}
"
@

<<etude_parametres_modele_basique,echo=False,figure=True>>=
#Lecture données et variables
#Codification des données
data.jags = list(Np=nrow(data), etude=etude2, time = Ti,C=C,censure=censure)


#Initialisation des chaines de Markov 
#Définition des valeurs initiales des noeuds entrants pour chaque chaine
init_params = list(
  list(sigma=2, mu=-15),
  list(sigma=4, mu = 0)
)

#Simulation des chaines de Markov
#Fonction qui va vérifier l'intégrité du modèle codé.
#m = jags.model(file=textConnection(ms),  data = data.jags, inits = init_params,  n.chains = 2)
#update(m, 3000)
#mcmc2 <- coda.samples(m, c("mu","sigma"), n.iter = 8000)
@


Ce premier modèle nous donne un intervalle de confiance à 95 \% plus resseré pour la détermination du paramètre $ log(\theta) \in [-4.454,-3.921] \Leftrightarrow \theta \in [0.012,0.020]$. Mais le temps de calcul est très long, ce qui posera problème par la suite pour estimer la puissance du modèle. 



Ce premier modèle n'étant pas optimisé, alors qu'il est possible de calculer directement à la main une loi de probabilité maximisant la vraisemblance, on a décidé d'en développer un autre s'appuyant sur le \textbf{modèle du zéro trick}. En effet, les statistiques bayésiennes s'appuient sur la formule suivante: $$ P(y|\theta)= \dfrac {P(\theta|y)P(y)}{P(\theta)}\Leftrightarrow P(y|\theta)\propto P(\theta|y)\times P(y)$$ avec dans notre cas $ P(y)$ qui représente la distribution a priori de notre paramètre (ici suivant une loi uniforme) et $ P(\theta|y)$ qui représente la vraisemblance de nos données (on cherche à déterminer une distribution collant au plus près de nos données). 

Pour une loi exponentielle, la vraisemblance se calcule de la façon suivante: 
$$L(x_{1},...x_{n},\theta)= \prod_{i=1}^{n}\theta \exp ^{-\theta x_{i}} \Leftrightarrow log(L(x_{1},...x_{n},\theta))= n(event=1)log(\theta) -\theta \times \sum_{i=1}^{n}x_{i}$$, avec ici n le nombre de personnes par étude, et subtilité introduite par la censure, la somme des temps correspond pour chaque patient soit à son temps d'observation total si censuré, soit au temps d'apparition de la maladie si non censuré. On a supposé encore ici que chaque réalisation était indépendante et identiquement distribué, sinon ce calcul ne fait pas sens. Or, on peut réinterpréter cette vraisemblance comme une probabilité globale combinant les probabilités individuelles de chaque individu d'apparition de la maladie. De plus, la loi de Poisson: $$ p(k)=P(X=k)=\dfrac{\lambda ^{k} \exp ^{-\lambda}}{k!}$$ a pour paramètre $ \lambda$ qui correspond au nombre d'événements attendus pour un certain laps de temps.Si on associe le nombre de ces évènements au nombre de patients malades au moment de la censure, alors on peut considérer que la probalité qu'aucun patient ne soit malade au sein d'une étude: $(P(X=0)=\exp^{-\lambda})$ suit une loi de Poisson de paramètre $ \lambda$ égale à moins la vraisemblance déterminée précédemment.En prenant un autre point de vue, on peut assimiler en discret cette expérience à une loi binomiale de paramètres (p,n) avec $p=\dfrac{Tobs}{Ttotal}$ la probabilité pour chaque chaque patient qu'il développe une maladie avant la fin de la censure, et n le nombre de patients total. Le paramètre $\lambda$ est égale à $ p\times n$ et aussi égale à l'espérance de la loi de Poisson. Or on attend bien en moyenne $n\times \dfrac {Tobs}{Ttotal} $ ou autrement dit le nombre total de patients multiplié par la probabilité pour chaque patient d'être malade au moment de la censure (supposé être égal pour chacun). Et la vraisemblance ne fait rien que représenter cette distribution. En voici le DAG résumé: 

\begin{figure}[H]

\begin{center}
\includegraphics{modele_zero_trick.png} 
\end{center}

\end{figure}



Et son codage sous R:
<<>>=
#Nombre de patients par étude
n1 = tapply(event,etude,function(x) {sum(x)})
##somme des temps(censure ou apparition) par étude
sum_ti = tapply(time,etude,function(x) {sum(x)})
zero = rep(0,5)

ms_zerotrick = " 
model{
C = 10000
for(i in 1:nb_etude){

zero[i] ~ dpois(-logvrai[i] +C)
logvrai[i] <- n1[i]*log(teta[i]) - teta[i]*sum_ti[i]

logteta[i] ~ dnorm(mu,tau)
teta[i] = exp(logteta[i])
}

tau = 1/(sigma**2)
sigma ~ dunif(0,20)
mu ~ dunif(-20,10)
}
"
@


Cette méthode présente l'avantage de prendre en compte le cas spécifique de la censure, tout en nous évitant de prendre une valeur fixée pour la vraisemblance. On aurait en effet pu construire un modèle dans lequel on a une estimation de $ \theta$ correspondant directement au maximum de vraisemblance, mais dans ce cas, on n'est plus dans un cadre bayésien mais inférentiel. Par rapport à la méthode précédente, elle est aussi plus stable et indépendante du choix des paramètres initiaux.  Comme on peut le voir sur les graphiques ci dessous, le critère de Gelman et Rubin, ou \textbf{indice de réduction de la variance}, calculé pour chaque paramètre comme rapport entre la variance globale sur la variance intra chaînes, indique une bien meilleure convergence de notre chaîne, même pour des vlaurs plus faibles d'itérations avec notre nouvelle méthode.  Il doit notamment être proche de 1 puisque la variance inter-chaîne doit être négligeable par rapport à la variance intra-chaîne. Numériquement, on constate qu'il est de 1.07 contre 1 pour le modèle zéro-trick; graphiquement, on observe une convergence plus rapide autour de 1, et ce dès les premières itérations avec le modèle zéro-trick. De plus, l'autocorrélation est plus faible avec le modèle zéro-trick, et surtout peut être réglé avec un décalage ou \textbf {thin} de 4 (signifiant que l'on ne conserve qu'une valeur sur quatre pour l'estimation des paramètres finaux).On peut notamment le voir avec le graphe d'autocorrélation, mais surtout la fonction raftery qui indique à la fois le nombre minimum d'itérations à réaliser pour suffisament converger, la taille efficace ou limite inférieure(en enlevant les données trop corrélés) et le facteur de dépendance à appliquer. 


<<echo=False>>=
#Codification des données
data.jags = list(zero=zero,n1=n1,sum_ti=sum_ti,nb_etude=5)


#Initialisation des chaines de Markov 
#Définition des valeurs initiales des noeuds entrants pour chaque chaine
init_params = list(
  list(sigma=2, mu=-15),
  list(sigma=4, mu = 0),
  list(sigma=6, mu =5)
)

#Simulation des chaines de Markov
#Fonction qui va vérifier l'intégrité du modèle codé.
m_1 = jags.model(file=textConnection(ms_zerotrick),
                 data = data.jags,
                 inits = init_params,
                 n.chains = 3)


update(m_1, 3000)
mcmc1 <- coda.samples(m_1, c("mu","sigma"), n.iter = 8000)
@

<<plot_comparaison,echo=False,fig=True>>=
par(mfrow = c(1, 2))
require (lattice)
#densityplot(mcmc2)
densityplot(mcmc1)
@

<<autocorrelation_comparaison,echo=False,fig=True>>=
par(mfrow = c(1, 2))
#autocorr.plot(mcmc2)
autocorr.plot(mcmc1)
#raftery.diag(mcmc2)
raftery.diag(mcmc1)
summary(mcmc1)
@




<<echo=False>>=
update(m_1, 3000)
mcmc1 <- coda.samples(m_1, c("mu","sigma"), n.iter = 16000,thin=4)
par(mfrow = c(1, 1))
print("autocorrélation après application d'un décalage de 4")
autocorr.plot(mcmc1)

@


On peut aussi effectuer d'autres types de calcul, comme la détermination du temps moyen d'apparition d'un événement cardiovasculaire égale à $ \dfrac {1}{\theta}$. Pour cela, puisque les valeurs de la chaine mcmc correspondent à la distribution réelle de $\theta$, il suffit de calculer la variable notre moyenne pour chaque valeur de notre chaîne, puis de déterminer un quantile pour avoir un intervalle de confiance à 95\%. On suit le même raisonnement pour obtenir une estimation (médiane de la distribution) de $ \sigma$ et $\theta$. 

<<>>=
mcmtot1=as.data.frame(as.matrix(mcmc1))
mu_simule=mcmtot1[,"mu"]
#estimation de log theta moyen
mu_mediane=quantile(mu_simule,probs=0.5) 
#intervalle de confiance de theta
distribution_theta=exp(mu_simule)
print(quantile(distribution_theta,probs=c(0.025,0.975)))
#intervalle de confiance de 1/theta
distribution_temps_moyen=1/distribution_theta
print(quantile(distribution_temps_moyen,probs=c(0.025,0.5,0.975)))
#intervalle de sigma moyen
sigma_simule=mcmtot1[,"sigma"]
print(quantile(sigma_simule,probs=c(0.025,0.975)))
@
Ainsi, en moyenne, le patient va développer la maladie cardiovasculaire en 66 ans. Cette estimation montre ici bien l'importance de la prise en compte de la censure (nombreuses données manquantes), mais aussi l'utopie de considérer avoir réellement uen distribution complète des valeurs du paramètre en poursuivant l'étude, puisque de toute façon nombre de patients seront déjà morts avant même de développer la maladie.  

On peut de plus comparer l'information apportée à priori avec l'information apportée par les données (ici volontairement très faible). On peut voir ci dessous que les données sont très informatives et permettent de resserrer très nettement les lois
a priori. Pour dessiner les lois à prioro, 2 possibilités s'offrent à nous: on peut décider de faire tourner le modèle à vide (mais la fonction density a tendance à ne pas représenter correctement les distribtions uniformes), ou plus simple dans le cas d'une loi uniforme tracer la droite d'équation $ y=\dfrac {1}{b-a}$, avec (b-a) les bornes de défintion de notre loi.
<<fig=True,echo=False>>=
par(mfrow = c(1, 2))

hist(mcmtot1[,"mu"], main = "Apport des données pour la distribution de logtheta", xlab = names(mcmtot1)[1], freq = FALSE,font.main=4,cex.main=0.8)
abline(h=1/30, col = "blue", lwd = 2)
hist(mcmtot1[,"sigma"], main = "Apport des données pour la distribution de sigma", xlab = names(mcmtot1)[2], freq = FALSE,font.main=4,cex.main=0.9)
abline(h=1/20, col = "blue", lwd = 2)
@


Enfin, une dernière possibilité est d'utiliser la fonction diffdic pour comparer l’ajustement aux données avec le critère DIC (Deviance Information Criterion), notamment pour déterminer le modèle le plus ajusté aux données, tout en évitant le "overfitting" car prenant en compte le nombre de paramètres nécessaires pour faire tourner le modèle. Voici le code qui permet sur notre exemple de calculer le DIC classique sur 2000 itérations. Plus le DIC est faible, plus le modèle est ajusté aux données pondéré par le nombre d eparamètres sélectionnés. 
<<>>=
#dic_ancien=dic.samples(m, n.iter = 2000,thin=4)
dic_nouveau=dic.samples(m_1, n.iter = 16000,thin=4)
print(dic_nouveau)
#diffdic(dic_ancien,dic_nouveau)
@

En conclusion, que ce soit par la rapidité des calculs ou par la mesure de qualité des résultats obtenus, nous allons nour concentrer sur le modèle zero trick pour le reste du projet. 

\section{Simulation d'une expérience}
Dans cette partie, nous avons simulé les résultats d'une nouvelle étude en incluant n  patients pour chaque bras de traitement : un expérimental et un de référence. Une censure administrative à 3 ans est considérée à partir du premier individu inséré dans l'étude. Pour être au plus proche de la réalité, on suppose que l'intégration des 200 patients s'effectue de façon aléatoire entre les 2 bras de traitement, et que leur arrivée est échelonnée uniformément sur la première année. Ensuite, au bout des 3 années de suivi, tous les patients ayant développé une maladie cardiovasculaire sont considérés comme positifs (événement 1), leur temps étant celui de l'incubation de la maladie depuis leur entrée à l'hôpital. Pour les autres, ils sont considérés comme censurés (événement 0), et le temps indiqué est alors celui de l'observation totale sur l'individu (temps censure -date entrée).\\
De plus, pour une étude considérée,tous les patients sont générés à partir du même theta tirés de la loi normale de paramètres $ \mu \sigma$ déterminés à partir de l'estimation médiane de ces paramètres dans l'étude globale antérieure. Ainsi, comme auparavant, on suppose ce dernier constant pour une population donnée, supposant ainsi que le choix des patients et la qualité des traitements fournis est homogène dans un même hôpital. Le $ \sigma$ est encore ici présent, mais pour indiquer une variabilité possible entre des études différentes, et non au sein d'une même population. Ce choix peut poser problème:ainsi, si le theta est biaisé dès le départ, même en augmentant le nombre de patients, on n’aura une estimation plus précise mais jamais juste: on sera ainsi peut-être plus précis autour d’une valeur qui est « fausse ».\\
Toutefois, ce qui va nous intéresser plus particulièrement ici est la capacité à distinguer une différence significative entre le $ \theta_{exp} $ et le $ \theta_{ref}$, mesure donnée par le hazard ratio: $ HR= \dfrac{\theta_{exp}}{\theta_{ref}}$. Ainsi, on suppose que $ \theta$ est similaire entre les deux études (puisque le recrutement est similaire), la seule différence étant porté par une éventuelle variation du hazard ratio lié à l'effet supplémentaire ajouté par le traitement expérimental. On distingue deux cas :
\begin {itemize}
\item Sous $ H_{0}$, on suppose que les différences entre les deux modèles sont significatives. ce seuil est considéré comme tel si le hazard ratio est supérieure ou égale à 1.3. Dans ce cas, on conclut systématiquement que la différence de développement de maladies est signficative entre les deux traitements. On teste alors: $$ HR=1.3$$.
\item Sous $ H_{1}$, on suppose que les différences entre les deux modèles ne sont pas significatives. Dans ce cas, Hr doit être inférieure à 1.3, ce qui revient à déterminer la probabilité que $$ HR \leq 1.3$$. Mais pour le calcul de puissance, il est nécessaire de fixer une valeur de HR bien définie (ici on prend un HR de 1). Il serait éventuellement possible de considèrer une petite incertitude sur la valeur attendue sous H1 ;mais pour des soucis de simplification, on calculera cette puissance en supposant que les deux bras de traitement ont le même $\theta $, et donc que HR est égale à un. 
\end{itemize}


\section {Description du modèle hiérarchique et du modèle classique}
Une fois différentes siumlations lancées (dont on peut faire varier l'hypothèse sous jacente, le nombre de patients dans chaque bras de traitement ou encore l'écart type), on va chercher à retrouver la valeur du hazard ratio, et surtout à déterminer un intervalle de confiance de cette valeur. 

\begin {itemize}
\item \textbf{Modèle hiérarchique} Deux possibilités nous sont sont offertes: soit on n'utilise pas les données à priori ( utilisation des lois vaguement informatives) en considérant que le temps moyen d'incidence (paramètre $ \theta $) est similaire pour les deux études et seule la vriation du hazard ratio va expliquer les différences d'effectifs constatés entre les deux groupes. En effet, si on a suivi une randomisation correcte de l'échantillon de départ, nous sommes normalement en présence de deux groupes homogènes entre les deux bras de traitement. Ou on peut aussi utiliser les données de l'ensemble des 6 études (les 5 études précédentes ajoutés à notre étude personelle), en ajoutant une variable indiquant s'il s'agit d'un traitement expérimental ou de référence. Pour calculer HR, on suppose que :$$\theta[f_{exp}] = \theta[f_{ref}] \times HR$$.
\item \textbf{Modèle non hiérarchique} Dans ce cas, on intégre directement les études antérieures dans la distribution à priori de $ log(\theta)$ qui va donc suivre une loi normale de paramètre $ \mu$ et $\sigma$. 
\end{itemize}

Voici les deux DAGs correspondant à ces deux modèles:

\begin{figure}[H]

\begin{center}
\includegraphics{modele_non_hierarchique.png} 
\end{center}

\end{figure}


\begin{figure}[H]

\begin{center}
\includegraphics{modele_hierarchique.png} 
\end{center}

\end{figure}

Suite à ces simulations, on peut réaliser différentes études pour comparer ou justifier la validité d'un modèle. On peut ainsi récupérer les chaînes MCMC, afin de simuler la distribution prédictive des données et visualiser la qualité d’ajustement du modèle aux données, en comparant les valeurs prédites (issues de la distribution a posteriori) aux valeurs observées.
Ci-dessous, on compare les intervalles de crédibilité du nombre de patients attendus à la censure en fonction du bras de traitement, que l'on va comparer avec les données observées. Pour cela, on concatène les informations de 5 simulations réalisées sous $ H_{0}$. 

<<echo=FALSE>>=
simulation_clinique <- function(h0=FALSE,sigma=0.2,mu=mu_mediane,n=100,temps_administratif=3.0) {
  
  #si on est sous h0, on suppose que les traitements ont un impact significatif sur le hazard ratio
  #de 1.3 sinon pas d'impact, ce qui est modelise par une meme probabilite de taux d'evenement
  if (h0) {
    mu_reference=exp(rnorm(1,mu,sd=0.2))
    mu_experience=mu_reference*1.3
    
  }
  else {
    #tirage du taux moyen d'apparition d'evenements=parametre de la loi exponentielle
    mu_reference=exp(rnorm(1,mean=mu,sd=0.2))
    mu_experience=mu_reference
  }

  
  #tirage de l'entree des patients dans l'etude au cours de la première annee
  entree_reference=runif(n,min=0,max=1)
  entree_experimental=runif(n,min=0,max=1)
  
  #debut etude 
  debut_etude=min(c(entree_reference,entree_experimental))
  fin_etude=debut_etude+temps_administratif
 
  
  
  #simulation loi exponentielle de parametre exponentielle de mu
  evenement_reference=entree_reference+rexp(n,mu_reference)
  evenement_experimental=entree_experimental+rexp(n,mu_experience)

  #creation du data frame type de jeu de données
  #0 pour experimental et 1 pour ancien traitement
  data<- data.frame(etude=c(rep("fancien",n),rep("fnew",n)),time=c(evenement_reference,evenement_experimental),event=rep(0,2*n),bras_traitement=c(rep(0,n),rep(1,n)),temps_entree=c(entree_reference,entree_experimental),stringsAsFactors = FALSE)
  

 #mise a jour des vecteurs d'evenement (0 si censuré, 1 sinon)
  data$event[data$time<=fin_etude] = 1
 
  
  #mise a jour des temps sur la base du data source
  #temps=temps-temps_entree si evenement, sinon temps_administratif+debut inclusion=duree totale etude-temps d'introduction correspondant ainsi au temps d'observation
  data$time[data$event==0]=fin_etude-data$temps_entree[data$event==0]
  data$time[data$event==1]=data$time[ data$event==1]-data$temps_entree[ data$event==1]
 
  
 
  return(data[1:4])
}
@

<<echo=FALSE>>=
modele_hierarchique<-function(nb_personne_par_simulation=20000,h0=FALSE,a_priori=NULL,sigma=0.2,type="hazard_ratio",a_vide=FALSE) {
  
  if (!is.null(a_priori)) {
    #on concatene nos deux tables
    bras_temp=data.frame(bras_traitement=rep(0,nrow(a_priori)), stringsAsFactors = FALSE)
    a_priori=cbind(a_priori,bras_temp)
  }
 
    donnne_par_simulation=simulation_clinique(h0=h0,sigma=sigma,n=nb_personne_par_simulation)
    print(nrow(donnne_par_simulation))
  if (is.null(a_priori)) {
      #on repete notre modele pour chaque simulation
      effectif_bras_traitement=tapply(donnne_par_simulation$event,donnne_par_simulation$etude,function(x) {sum(x)})
      sum_ti = tapply(donnne_par_simulation$time,donnne_par_simulation$etude,function(x) {sum(x)})
      dcomplet <- list( sum_ti=sum_ti,effectif_par_experience=effectif_bras_traitement, zero=c(0,0),sigma=sigma,nb_etude=2,type_traitement=c(0,1),indice=c(1,1))
      if (a_vide) {
        dcomplet <- list( sum_ti=sum_ti,effectif_par_experience=effectif_bras_traitement,sigma=sigma,nb_etude=2,type_traitement=c(0,1),indice=c(1,1))
      }
      
      
    }
   else {
     
     
     data_complete=rbind(a_priori,donnne_par_simulation)
     effectif_par_experience=tapply(data_complete$event,data_complete$etude,function(x) {sum(x)})
     sum_ti = tapply(data_complete$time,data_complete$etude,function(x) {sum(x)})
     dcomplet <- list( sum_ti=sum_ti,effectif_par_experience=effectif_par_experience, zero=rep(0,7),sigma=sigma,nb_etude=7,type_traitement=c(rep(0,6),1),indice=c(seq(1,6),6))
     if (a_vide) {
       dcomplet <- list( sum_ti=sum_ti,effectif_par_experience=effectif_par_experience,sigma=sigma,nb_etude=7,type_traitement=c(rep(0,6),1),indice=c(seq(1,6),6))
     }
    
   }

      ms3 = "
      model{
     C = 10000
     for(i in 1:nb_etude){
      temp[i]~dnorm(mu,tau)
      logteta[i]=temp[i]*(1-type_traitement[i])+temp[indice[i]]*type_traitement[i]
     zero[i] ~ dpois(-logvrai[i] + C)
     teta[i]=exp(logteta[i])*hazard_ratio*type_traitement[i] + exp(logteta[i])*(1-type_traitement[i])
     logvrai[i] =effectif_par_experience[i]*log(teta[i]) - teta[i]*sum_ti[i]
     }
     #lois a priori mais avec sigma fixee
     #on ne cherche pas a estimer la variabilite intra etude, on la suppose egale a 0.2 aussi
     tau = 1/(sigma**2)
     mu ~ dunif(-20,10)
     
     
     log_hazard_ratio~dunif(-10,10)
     hazard_ratio =exp(log_hazard_ratio)
     
   }
     "
    ini <- list(list( mu=5,log_hazard_ratio=-5),list(mu=-5,log_hazard_ratio=5),list(mu=2,log_hazard_ratio=-7))
    m3<-jags.model(file = textConnection(ms3), data = dcomplet, inits = ini, n.chains = 3,quiet=TRUE)
    update(m3, 3000)
    mcmc3 <- coda.samples(m3, type, n.iter = 20000, thin = 5)
    
    mcmtot1=as.data.frame(as.matrix(mcmc3))
  return(mcmtot1)
}


@

<<echo=FALSE>>=
modele_non_hierarchique<-function(nb_personne_par_simulation=100000,h0=FALSE,sigma=0.2,mu=mu_mediane,a_priori=NULL,type="hazard_ratio",a_vide=FALSE) {
  
    #on repete notre modele pour chaque simulation
    donnne_par_simulation=simulation_clinique(h0=h0,sigma=sigma,n=nb_personne_par_simulation)
    n1=tapply(donnne_par_simulation$event,donnne_par_simulation$etude,function(x) {sum(x==1)})
    sum_ti = tapply(donnne_par_simulation$time,donnne_par_simulation$etude,function(x) {sum(x)})
    dcomplet <- list( sum_ti=sum_ti,n1=n1, nb_etude=2, zero=c(0,0),sigma=sigma,mu=mu_mediane,type_etude=c(1,0))
    if (a_vide) {
      dcomplet <- list( sum_ti=sum_ti,n1=n1, nb_etude=2,sigma=sigma,mu=mu_mediane,type_etude=c(1,0))
    }
    
  
  ms_zerotrick_nonhierarchical_simulation = "
  model{
  C = 10000
  for(i in 1:nb_etude){
  zero[i] ~ dpois(-logvrai[i] + C)
  logvrai[i] =n1[i]*log(teta[i]) - teta[i]*sum_ti[i]
  teta[i] = exp(logteta)*type_etude[i]+exp(logteta)*(1-type_etude[i])*hazard_ratio
  }
  #lois a priori
  tau = 1/(sigma**2)
  logteta~dnorm(mu,tau)

  log_hazard_ratio~dunif(-10,10)
  hazard_ratio =exp(log_hazard_ratio)
  }
  "

ini <- list(list( logteta=-15,log_hazard_ratio=-5),list(logteta=0,log_hazard_ratio=-3),list(logteta =5,log_hazard_ratio=4))
m3<-jags.model(file = textConnection(ms_zerotrick_nonhierarchical_simulation), data = dcomplet, inits = ini, n.chains = 3,quiet=TRUE)
update(m3, 3000)
mcmc3 <- coda.samples(m3, c(type), n.iter = 20000, thin = 5)
mcmtot1=as.data.frame(as.matrix(mcmc3))
  
  return(mcmtot1)
}


@



<<comparaison_modele,echo=False>>=
#####################calcul puissance par simulation########
creation_table_puissance <- function(fonction=modele_hierarchique,nb_personne_par_simulation=1000,sigma=0.2,mu=mu_mediane,nb_simulations=4,a_priori=NULL) {
  
  #test sous H0
  for (i in seq (nb_simulations)) {
    simulation=fonction(nb_personne_par_simulation=nb_personne_par_simulation,h0=TRUE,a_priori=data) 
    if (i==1) {
      hazard_ratio_simule=simulation
    }
    else {
      hazard_ratio_simule=cbind(hazard_ratio_simule,simulation)
    }
  }
  
  #recuperation du quantile à 95 % pour chaque colonne
  quantile_95_H0=apply(hazard_ratio_simulation,2,function(x) {quantile(x,probs=0.95)})
  #determination d'un intervalle de confiance
  largeur_intervalle_confiance_H0=log(apply(hazard_ratio_simulation,2,function(x) {quantile(x,probs=0.975)}))-log(apply(hazard_ratio_simulation,2,function(x) {quantile(x,probs=0.025)}))
  precision_H0=mean(largeur_intervalle_confiance_H0)
  #determination du biais, sous H0, c'est le log ecart entre le hazard ratio estime et le hazard ratio prevu (ici 1.3)
  biais_introduit_H0=log(hazard_ratio_simulation/1.3) #revient a appliquer pour chauqe element de la mtrice l'operation log(HR)-log(Hrprevu)=log(HR/HRprevu)
  quantile_biais_h0=apply(biais_introduit_H0,2,function(x) {quantile(x,probs=0.5)})
  moyenne_biais_H0=mean(quantile_biais_h0)
  false_negatif=sum(quantile_95_H0<1.3)
  true_positif=sum(quantile_95_H0>1.3)
  
 
  
  
  #test sous H1
  for (i in seq (nb_simulations)) {
    simulation=fonction(nb_personne_par_simulation=nb_personne_par_simulation,h0=FALSE,a_priori=data) 
    if (i==1) {
      hazard_ratio_simulation2=simulation
    }
    else {
      hazard_ratio_simulation2=cbind(hazard_ratio_simulation2,simulation)
    }
  }
  #recuperation du quantile à 95 % pour chaque colone
  quantile_95_H1=apply(hazard_ratio_simulation2,2,function(x) {quantile(x,probs=0.95)})
  #determination precision H1
  largeur_intervalle_confiance_H1=log(apply(hazard_ratio_simulation2,2,function(x) {quantile(x,probs=0.975)}))-log(apply(hazard_ratio_simulation2,2,function(x) {quantile(x,probs=0.025)}))
  precision_H1=mean(largeur_intervalle_confiance_H1)
  #determination du biais, sous H1 (revient a calculer le log car log(1)=0)
  biais_introduit_H1=log(hazard_ratio_simulation2) #revient a appliquer pour chauqe element de la mtrice l'operation log(HR)-log(Hrprevu)=log(HR/HRprevu)
  quantile_biais_h1=apply(biais_introduit_H1,2,function(x) {quantile(x,probs=0.5)})
  moyenne_biais_H1=mean(quantile_biais_h1)
  true_negative=sum(quantile_95_H1<1.3)
  false_positive=sum(quantile_95_H1>1.3)
  
  
  matrice_globale=matrix(data = c(true_positif,false_negatif,false_positive,true_negative),nrow=2,ncol=2,dimnames = list(c("H0 predits","H1 predits"),c("H0 reels","H1 reels")))
  table_contingence=as.table(matrice_globale)
  
  #accuracy mesure le nombre total de bonnes predictions sur le nombre total d'observations
  accuracy=(true_negative+true_positif)/(nb_simulations*2)
  
  #precision donne le pourcentage de vrais bonnes détections parmi les détections postives détectées par notre modele d'apprentissage
  precision=true_positif/(true_positif+false_positive)
  
  #sensibilite mesure la probabilite que l'on conclut à un HR significatif si HR est reellement significatif
  sensibilite=true_positif/(true_positif+false_negatif)
  
  #specificite mesure la probabilite de conclure a une egalite de l'effet des traitements alors qu'effectivement les traitements
  #ont le meme effet sur le taux d'incidence de maladies cardio
  specificite=true_negative/(true_negative+false_positive)
  
  #la f-measure est utilise en machine learning, prend en compte la mesure des deux types d'erreur possibles

  #fausse detection d'effet ou non detection de l'effet
  F_measure=(precision*sensibilite*2)/(precision+sensibilite)
  
  statistiques_generales=list(table_contingence=table_contingence,accuracy=accuracy,sensibilite=sensibilite,specificite=specificite,F_measure=F_measure,precisionH1=precision_H1,precisionH0=precision_H0,moyenne_biais_H0=moyenne_biais_H0,moyenne_biais_H1=moyenne_biais_H1,precision=precision)
  
  return(statistiques_generales)
}

creation_table_puissance(nb_simulations=10,a_priori = data,nb_personne_par_simulation=2000)
@

\section {Comparaison entre le modèle hiérarchique et le modèle non hiérarchique}
Pour l'étude de l'effet du choix de la variance ou de la taille du bras sur la précision des paramètres en sortie, on se référera à l'étude auparavant réalisée par nos collègues. Nous allons donc fixer dans toute notre étude le nombre de patients dans chaque bras à 1000 et la variance à 0.2, valeur précédemment estimée à l'aide des données a priori. De plus, puisque dans le modèle non hiérarchique, on inclut directement la distribution a priori, on fera ici de même pur le modèle hiérarchique afin d'avoir un réel élément de comparaison. \\
Pour cela, on va utiliser 6 critères différents:
\begin {itemize}
\item Le risque de première espèce. C'est la probabilité de se tromper et d'affirmer que l'on est sous $ H_{1}$, alors qu'en réalité on est sous $ H_{0}$. Autrement dit, c'est la probabilté de considérer que les traitements ont le même effet sur le taux d'incidence des maladies alors qu'en réalité, les traitements sont significativement différents. La \textbf {sensibilité} n'est autre qu'une mesure de la capacité du modèle à détecter une différence significative lorsqu'elle existe. Elle est égale à $ 1-\alpha$. Ce critère est essentiel, et justifie la forme atypique des hypothèses du test. En effet, si on se trompe sur cette hypothèse, on peut éventuellement conclure à l'absence d'effet du traitement sur le taux d'incidence, avec des conséquences graves. On pourrait en effet distribuer un traitement qui s'avère finalement plus dangereux que l'ancien en termes de bénéfice/risque et éventuellement conduisant à une surmortalité au sein de notre population traitée. Celle-ci, lié au relativement haut temps moyen d'apparition de la maladie ne pourrait bien être détectée que des années après le début du traitement. En effet, la moyenne d'une loi exponentielle est de $ \dfrac {1}{\theta}$, une augmentation de $ \theta$ impliquerait donc une diminution de ce temps moyen et in fine, une apparition plus rapide et importante du nombre de maladies cardiovasculaires. Généralement, en statistiques, on considère un seuil de 5 \% pour réfuter une hypothèse (même si dans le domaine de la santé, où les conséquences sont éventuellement plus importantes, on tend à diminuer ce dernier). De plus, avec les chaines mcmc, nous avons la distribution réelle de la densité de probabilité de $ \theta$. On calcule alors le quantile à 95 \%: si ce dernier est inférieure à 1.3 c'est que la probabilité que $ \theta$ soit inférieure 1.3 est supérieure à 95 $ \%$ ou autrement dit la probabilité d'être sous $ H_{0}$ est inférieure au seuil de conservation de $ 5\%$. Dans ce cas on réfute $ H_{0}$ alors qu'on est bien sous $H_{0}$, on est donc en présence d'un faux négatif. Le risque de première espèce se détermine donc ainsi: $ puissance= \dfrac {true postive}{nombre total d'experiences sous h0}$


\item La puissance ou spécificité. C'est la capacité à détecter que l'on est sous $ H_{0}$ alors qu'en réalité on est bien sous $ H_{0}$. Autrement dit, elle reflète la capacité de notre modèle à détecter lorsque les traitements sont siginifativent égaux en terme d'incidence des mladies alors qu'ils le sont réellement. Le risque de seconde espèce, nommé $ \beta$ correspond à 1-spécificité. Cette mesure a aussi son importance, puisqu'en la rejetant alors qu'elle est vraie, on pourrait potentiellement rejeter un nouveau traitement soulageant plus efficacement les malades, sous un prétexte erroné qu'elle augmenterait le taux d'incidence des maladies cardiovasculaires. Pour son calcul, on suit le même raisonnement que précédent: si le quantile à 95  \% est inférieure à 1.3, c'est que la probabilité d'être sous $H_{0}$ est inférieure à 5\%, on conclut donc à l''hypothèse $ H_{1}$. On est en présence d'un vrai négatif. Pour déterminer la puissance, il suffit de calculer le ratio suivant : $$ puissance= \dfrac {\text{true negative}}{\text {true negative + false positive}} \Leftrightarrow puissance= \dfrac {\text{true negative}}{\text {true negative + false positive}} $$. Plus le nombre de simulations sera importante, tous les autres paramètres étant fixés par ailleurs, plus notre calcul de puissance sera robuste. Pour illuster notre propos, cf graphique ci dessous.




<<fig=True,echo=True>>=
nb_simulations=10
  for (i in seq (nb_simulations)) {
    simulation=modele_hierarchique(nb_personne_par_simulation=200,h0=FALSE,a_priori=data) 
    if (i==1) {
      hazard_ratio_simule=simulation
    }
    else {
      hazard_ratio_simule=cbind(hazard_ratio_simule,simulation)
    }
  }
 
  quantile_95_H1=apply(hazard_ratio_simule,2,function(x) {quantile(x,probs=0.95)})
  indices <- which(quantile_95_H1<1.3)
  for (indice in indices) {
    hist(hazard_ratio_simule[,indice],main="distribution de probabilité du hazard ratio",freq = FALSE,xlab="hazard ratio")
    abline(v=1.3,col="red")
    
    print( paste("la probabilite de se tromper est egale a ",toString(sum(hazard_ratio_simule[,indice]>1.3)/nrow(hazard_ratio_simule))))
    }
@



\item L'accuray, plutôt utilisée en machine learning, mesure le nombre de bonnes détections sur le nombre total de simulations. Elle traduit donc la probabilité, toutes hypothèses confondues, que le modèle détermine la bonne hypothèse.

\item On peut aussi trouver la f-measure, qui combine la sensibilité avec un critère nommé précision en data mining (à ne pas confondre avec la précision décrite ci après). Elle se mesure ainsi: $$ F-measure= \dfrac {precision*sensibilite*2}{precision+sensibilite}$$, et est une mesure de la qualité globale de notre modèle à détecter les expériences sous $ H_{0}$, sans en ajouter.

\item La précision, ici au sens statistique, est le logarithme de la largeur de l'intervalle à 95 \% de notre paramètre $ \gamma$.  Pour cela, on détermine la largeur de cet intervalle pour chaque simulation, et on récupére la médiane de toutes les simulations. Plus la précision est faible, plus la largeur de l'intervalle est faible, donc plus on se ressere autour de la vraie valeur du hazard ratio, et plus les risques de première ou seconde espèce sont réduits.

\item Le biais est l'écart logarithmique médian attendu entre la valeur prédite (au sens médian) par notre modèle et la valeur attendue se lon notre hypothèse. Plus ce biais est faible, meilleur est notre modèle, puisque fidèle à la réalité.

\item La justesse d'un modèle, que l'on ne détermine pas directement est une combinaison de la précision et du biais. Si est on est précis (faible largeur de l'intervalle) et fidèle, alors meilleur et plus juste est notre modèle, et plus les chances d'aboutir à une conclusion correcte sont importantes. 

\end {itemize}


Voici les résultats obtenus avec 1 000 simulations, pour n=200 personnes puis n=1 500 personnes.
<<echo=False>>=
creation_table_puissance(fonction=modele_hierarchique,nb_personne_par_simulation=200,nb_simulations=1000,a_priori=data) 
print("resultats obtenus pour 200 personnes avec le modèle hiérarchique") 
creation_table_puissance(fonction=modele_hierarchique,nb_personne_par_simulation=1500,nb_simulations=1000,a_priori=data) 
print("resultats obtenus pour 1500 personnes avec le modèle hiérarchique") 
creation_table_puissance(fonction=modele_non_hierarchique,nb_personne_par_simulation=200,nb_simulations=1000) 
print("resultats obtenus pour 200 personnes avec le modèle non hiérarchique") 
creation_table_puissance(fonction=modele_non_hierarchique,nb_personne_par_simulation=1500,nb_simulations=1000) 
print("resultats obtenus pour 1500 personnes avec le modèle non hiérarchique") 
@


On voit notamment que selon les différents critères de distinction d'un modèle, nos choix vont se porter sur le modèle hiérarchique. Il présente notamment un écart moyen avec la valeur recherchée plus faible, une largeur de l'intervalle moins importante conduisant à une sensibilité et spécificté globalement meilleure. Cela peut s'expliquer notamment par le fait que dans un modèle non hiérarchique, l'essentiel de l'information est apportée par les données a priori et ne prend pas en compte les éventuelles spécificités de notre étude. Dans le même temps, avec un modèle hiérarchique, toutes les études sont au même niveau (auparavant le grand nombre de patients étudiés restreignait beaucoup les possibilités de variation de notre étude), ce qui nous permet de mieux prendre en compte la variabilité de réaction au traitement de référence, et donc d'avoir in fine une meilleure estimation de notre hazard ratio. On va donc plutôt porter nos choix sur un modèle hiérarchique pour mener une étude de puissance, ce dernier réduisant un peu plus les risques de première et seconde espèce. 



Une autre façon de le voir est de faire tourner les deux modèles à vide, et de comparer la distribution a priori avec la distribution a posteriori.
<<echo=False>>=

par(mfrow = c(1, 2))
mcmctot=modele_non_hierarchique(a_vide = FALSE,nb_personne_par_simulation=200,type="logteta") 
mcmctot0=modele_non_hierarchique(a_vide = TRUE,nb_personne_par_simulation=200,type="logteta") 
hist(mcmctot[,1], main = "Comparaison distribtution a priori et a posterio du HR avec le modèle non hiérarchique", xlab = "logteta",xlim=c(-5,-3), freq = FALSE)
lines(density(mcmctot0[,1]), col = "blue", lwd = 2)



mcmctot=modele_hierarchique(nb_personne_par_simulation=200,a_priori=data,sigma=0.2,type="logteta") 
mcmctot0=modele_hierarchique(nb_personne_par_simulation=200,a_priori=data,sigma=0.2,type="logteta",a_vide = TRUE) 
hist(mcmctot[,6], main = "Comparaison distribtution a priori et a posterio du HR avec le modèle hiérarchique", xlab = "logteta",xlim=c(-5,-3), freq = FALSE)
lines(density(mcmctot0[,6]), col = "blue", lwd = 2)
   


@

Les chaînes MCMC peuvent enfin être utilisées dans R pour simuler la distribution prédictive des données et construire ainsi des graphes appropriés à chaque cas pour visualiser la qualité d’ajustement du modèle aux données, en comparant les valeurs prédites par le modèle avec les valeurs observées. 




<<validation_modele_non hierarchique,fig=TRUE,echo=False>>=
#stockage des effectifs observes
nb_simulations=100
effectif_observe=c()
theta_obs=c()
for (i in seq (1,nb_simulations)) {
  simulation=simulation_clinique(h0=TRUE,n=1000,temps_administratif=3.0)
  effectif_bras_traitement=tapply(simulation$event,simulation$etude,function(x) {sum(x==1)})
  effectif_observe=c(effectif_observe,effectif_bras_traitement)
  theta_obs=c(theta_obs,c(1,1.3))
  
}

effectif_observe=as.numeric(effectif_observe)


prediction des effectifs
#stockage des quantiles `a 2.5, 50 et 97.5 %
 patients_predits=c()
qinf <- c();qsup <- c()
 theta_simule=c(1,1.3)
 for (j in seq(1,2)) {
   for (i in 1:nb_simulations) {
     if (j==1) {
       teta_distribution=modele_non_hierarchique(nb_personne_par_simulation=1000,h0=FALSE,type="teta")[,1]

     }
     else {
       teta_distribution=modele_non_hierarchique(nb_personne_par_simulation=1000,h0=TRUE,type="teta")[,2]
     }
   teta_predit=quantile(teta_distribution,0.5)
   #simulation temps attente patients
   temps_patients=rexp(1000,teta_predit)
   #on suppose un temps moyen d'oberservation de 2.5 ans
   nb_patients_malades=sum(temps_patients<2.5)
   patients_predits=c(patients_predits,nb_patients_malades)
   }
   qinf=c(qinf,quantile(patients_predits,0.025));qsup=c(qsup,quantile(patients_predits,0.975))
 }
 #representation graphique
 par(mfrow = c(1, 1))
 print(qinf)
 print(qsup)
 plot(effectif_observe ~ theta_obs,  pch = 16, xlab = "taux d'incidence",
 ylab = "nombre de patients malades a la cesure",main="validation du modele par comparaison avec les donnees observées",ylim=c(0,100))
 for (i in 1:2) {
 segments(theta_simule[i], qinf[i], theta_simule[i], qsup[i], col = "red", lwd = 2)
 }



#avec modele hierarchique
effectif_observe=c()
theta_obs=c()
for (i in seq (1,nb_simulations)) {
  simulation=simulation_clinique(h0=TRUE,n=1000,temps_administratif=3.0)
  effectif_bras_traitement=tapply(simulation$event,simulation$etude,function(x) {sum(x==1)})
  effectif_observe=c(effectif_observe,effectif_bras_traitement)
  theta_obs=c(theta_obs,c(1,1.3))
}

effectif_observe=as.numeric(effectif_observe)
#prediction des effectifs
##stockage des quantiles `a 2.5, 50 et 97.5 %
patients_predits=c()
qinf <- c();qsup <- c()
theta_simule=c(1,1.3)
for (j in seq(1,2)) {
  for (i in 1:nb_simulations) {
    if (j==1) {
      teta_distribution=modele_hierarchique(nb_personne_par_simulation=1000,h0=FALSE,type="teta",a_priori = data)[,6]
      print(head(teta_distribution))
     
    }
    else {
      teta_distribution=modele_hierarchique(nb_personne_par_simulation=1000,h0=TRUE,type="teta",a_priori = data)[,7]
    }
  teta_predit=quantile(teta_distribution,0.5)
  #simulation temps attente patients
  temps_patients=rexp(1000,teta_predit)
  #on suppose un temps moyen d'oberservation de 2.5 ans
  nb_patients_malades=sum(temps_patients<2.5)
  patients_predits=c(patients_predits,nb_patients_malades)
  }
  qinf=c(qinf,quantile(patients_predits,0.025));qsup=c(qsup,quantile(patients_predits,0.975))
}


#representation graphique
par(mfrow = c(1, 1))
print(qinf)
print(qsup)
plot(effectif_observe ~ theta_obs,  pch = 16, xlab = "taux d'incidence",
ylab = "nombre de patients malades a la cesure",main="validation du modele par comparaison avec les donnees observées",ylim=c(0,100))
for (i in 1:2) {
segments(theta_simule[i], qinf[i], theta_simule[i], qsup[i], col = "red", lwd = 2)
}



@


Comme par exemple dans la simulation ci dessus, où l'on va chercher à comparer les valeurs prédites des effectifs pour un nombre important de simulations avec les valeurs observées.


\end{document}