6.6 Applications
1) Le tableau 6.1 présente les résultats d’une étude de cohorte visant à déterminer, entre autres, les avantages de l’utilisation, comme outil de dépistage, de mesures de test de stress à l’exercice (EST). Dans de tels tests, un résultat de réussite / échec peut être obtenu lors du diagnostic de maladies coronariennes (CAD).
Tableau 6.1. Screening test for physical stress
| CAD | ||||
|---|---|---|---|---|
| Healthy | Ill | Total | ||
| EST | Negative | 327 | 208 | 535 |
| Positive | 115 | 815 | 930 | |
| Total | 442 | 1 023 | 1 465 | |
On supposera qu’il n’y a pas de biais de vérification:
a)
Sur la base de cette matrice de confusion, indiquez les valeurs suivantes (avec des intervalles de confiance à 95%): sensibilité et spécificité, valeur prédictive positive et négative.
b)
Quelle est la valeur de l’aire sous la courbe pour les données actuelles ?
La création de la table de données peut se faire comme indiqué ci-après, à partir d’une table créée avec la commande matrix() :
>tab <-as.la table (matrice(c(815,115,208,327), nrow =2, byrow=TRUE,
dimnames=list(EST=c(« +”, »−”),
CAD=c(« +”, »-”))))
>tab
Il est à noter que le tableau a été réorganisé de manière à correspondre aux notations habituelles, avec des événements « positifs” sur la première ligne (typiquement, l’exposition) et la première colonne (typiquement, la maladie) du tableau.
Il serait alors possible de calculer toutes les quantités requises en utilisant très peu d’opérations arithmétiques avec R. Cependant, le package epiR contient tous les outils nécessaires pour répondre aux questions épidémiologiques des tableaux 2 × 2. Ainsi, la commande epi.tests() fournit des valeurs de prévalence, de sensibilité /spécificité, ainsi que des valeurs prédictives positives et négatives:
>library (epiR)
>epi.tests (tab)
Concernant l’aire sous la courbe, une commande externe peut également être utilisée, roc.de.table(), dans le paquet epicalc:
>bibliothèque (epicalc)
>roc.de.tableau (tab, graphique = FAUX)
2) Le tableau 6.2 résume la proportion d’infarctus du myocarde observée chez les hommes âgés de 40 à 59 ans et pour lesquels le niveau de pression artérielle et le taux de cholestérol ont été mesurés, considérés sous forme de classes ordonnées.
Tableau 6.2. Infarctus et pression artérielle
| Cholestérol (mg / 100 ml) | |||||||
|---|---|---|---|---|---|---|---|
| BP | & lt; 200 | 200 – 209 | 210 – 219 | 220 – 244 | 245 – 259 | 260 – 284 | & gt; 284 |
| < 117 | 2/53 | 0/21 | 0/15 | 0/20 | 0/14 | 1/22 | 0/11 |
| 117 – 126 | 0/66 | 2/27 | 1/25 | 8/69 | 0/24 | 5/22 | 1/19 |
| 127 – 136 | 2/59 | 0/34 | 2/21 | 2/83 | 0/33 | 2/26 | 4/28 |
| 137 – 146 | 1/65 | 0/19 | 0/26 | 6/81 | 3/23 | 2/34 | 4/23 |
| 147 – 156 | 2/37 | 0/16 | 0/6 | 3/29 | 2/19 | 4/16 | 1/16 |
| 157 – 166 | 1/13 | 0/10 | 0/11 | 1/15 | 0/11 | 2/13 | 4/12 |
| 167 – 186 | 3/21 | 0/5 | 0/11 | 2/27 | 2/5 | 6/16 | 3/14 |
| > 186 | 1/5 | 0/1 | 3/6 | 1/10 | 1/7 | 1/7 | 1/7 |
The data are available in the file hdis.dat sous la forme d’un tableau avec quatre colonnes indiquant, respectivement, la pression artérielle (huit catégories, notées 1 à 8), le taux de cholestérol (sept catégories, notées 1 à 7), le nombre d’infarctus du myocarde et le nombre total d’individus. L’association entre la pression artérielle et la probabilité de souffrir d’un infarctus du myocarde est le point d’intérêt:
a)
Calculer les proportions d’infarctus du myocarde pour chaque niveau de pression artérielle et les représenter sous forme de tableau et sous forme graphique;
b)
Exprimer les proportions calculées en (a) sous forme logit;
c)
Sur la base d’un modèle de régression logistique, déterminer s’il existe une association significative au seuil α = 0,05 entre la pression artérielle, traitée comme une variable quantitative en considérant les centres de classe et la probabilité d’avoir une crise cardiaque;
d)
Exprimer en unités logit les probabilités d’infarctus prédites par le modèle pour chacun des niveaux de pression artérielle;
e)
Afficher sur le même graphique les proportions empiriques et la courbe de régression logistique en fonction des valeurs de pression artérielle (centres de classe).
Importation des données hdis.la dat est obtenue comme suit:
>bp < – lire.tableau (« hdis.dat », header=TRUE)
>str(bp)
Comme cela peut être vérifié, il n’y a pas d’étiquette associée aux niveaux des variables d’intérêt (bpress pour la pression artérielle, chol pour le taux de cholestérol). Pour générer et associer des étiquettes, les commandes suivantes peuvent être insérées:
> blab < – c(« <117”, »117-126”, »127-136”, »137-146”,
« 147-156”, »157-166”, »167-186”, »>186”)
> clab < – c(« <200”, »200-209”, »210-219”, »220-244”,
« 245-259”, »260-284”, »>284”)
> bp$bpress < – factor(bp$bpress, labels = blab)
> bp$chol <-factor(bp$chol, labels=clab)
La dernière commande convertit les variables d’origine en variables qualitatives et associe à leurs modalités les étiquettes définies par blab et clab. Pour vérifier que la base de données est maintenant sous la forme souhaitée, les commandes suivantes sont utiles :
>str(bp)
>résumé(bp)
Le tableau des fréquences relatives exposées dans l’instruction peut maintenant être reproduit:
>round(xtabs(hdis / total bpress + chol, data = bp), 2)
Comme nous allons nous concentrer uniquement sur la relation entre la pression artérielle et l’infarctus du myocarde, il est nécessaire d’agréger des données sur le taux de cholestérol. En d’autres termes, il est nécessaire d’accumuler les valeurs pour chaque niveau de pression artérielle, quel que soit le taux de cholestérol. Il faut également saisir l’opportunité de renommer les niveaux de la variable bpress en utilisant les centres des intervalles de classe :
>blab2 <-c(111.5, 121.5, 131.5,141.5,151.5,161.5,176.5,191.5)
> bp$bpress <-rep(blab2, each=7)
>dfrm< – aggregate(bp,
list(bpress=bp), sum)
C’est la dernière commande, aggregate(), qui permet d’agréger les données : elle additionne toutes les valeurs (somme) de la variable chol qui ne sont pas incluses dans la liste des variables à conserver pour analyse. En même temps, les résultats sont stockés dans une nouvelle base de données nommée dfrm. Notez que dans ce cas, nous n’utilisons pas la notation de formule mais celle de la syntaxe alternative pour aggregate() (variable de réponse suivie d’une liste composée d’une ou plusieurs variables explicatives). Un aperçu des données agrégées peut être obtenu en donnant head() comme entrée :
>head(dfrm, 5)
Lorsqu’une proportion p est disponible, sa valeur sur une échelle dont les unités sont des logits est donnée par la relation log(p/(1-p)). De là, les commandes suivantes peuvent être dérivées pour convertir les proportions, calculées comme hdis/ total en unités logit:
> logit < – function(x) log(x/(1-x))
> dfrm$prop < – dfrm$idh/dfrm$total
> dfrm$logit < – logit(dfrm$idh/dfrm$total)
Il doit être observé qu’une petite fonction a été définie pour la conversion des valeurs de x, ce qui est supposé être une des proportions, dans son équivalent log(x/(1 − x)). De même, on pourrait écrire:
>log((dfrmhdhdis/dfrmtotaltotal)/(1-dfrmhdhdis/dfrmtotaltotal))
Les valeurs calculées ont été ajoutées directement à la trame de données sous le nom prop et logit.
Le modèle de régression logistique s’écrit de la manière suivante:
>glm(cbind(hdis, total-hdis) bpress, data=dfrm, family=binomial)
La formulation utilisée, cbind(hdis, total-hdis) bpress, prend en compte le fait que nous accédons à des données groupées et non à des réponses individuelles. La commande glm() avec l’option family=binomial correspond à une régression logistique qui, sans trop entrer dans les détails, utilise par défaut la fonction logit comme lien canonique.
Les coefficients de régression peuvent être affichés en utilisant summary():
>summary(glm(cbind(hdis, total-hdis) bpress, data=dfrm,
family=binomial))
Le résultat précédent comprend les informations essentielles pour répondre à la question de la signification statistique de l’association entre la pression artérielle et la probabilité de crise cardiaque: le coefficient de régression (sur l’échelle log-odds) est égal à 0,024 et est significatif au seuil habituel de 5% (voir colonne Pr (>/z/)).
La probabilité d’avoir une crise cardiaque en fonction des différents niveaux de pression artérielle considérés est obtenue comme suit:
>glm.res<-glm(cbind(hdis, total-hdis) bpress, data=dfrm,
famille=binomial)
>predict(glm.res)
Il est à noter que les résultats intermédiaires générés par R ont été stockés sous le nom de glm.res avant d’utiliser la commande predict(). Les prédictions générées par R sont exprimées sous forme de logit et il est possible de comparer les logits observés et prédits :
>cbind(dfrm, logit.predit = prédire (glm.res))
Pratiquement tous les éléments sont disponibles pour représenter les proportions d’infarctus du myocarde observées et prédites graphiquement en fonction du niveau de pression artérielle. Les prédictions du modèle exprimées sous forme de proportions et non sur l’échelle logodds sont manquantes. Il peut ainsi être souhaitable de tracer la courbe de prédiction, c’est-à-dire la probabilité d’avoir une crise cardiaque en fonction de la pression artérielle, sans limiter cette dernière aux huit valeurs observées pour la variable bpress. Voici une solution possible :
>dfrmpropprop.predit < – predict(glm.res, type= »response »)
>f <- function(x)
1/(1 + exp(-(coef(glm.res) + coef(glm.res)*x)))
> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,
xlab = « Blood pressure”,ylab = « Infarction likelihood”,
panel = function(x, y, …) {
panel.xyplot(x, y, col = « gray30”, pch = 19, …)
panel.curve(f, lty = 3, col = « gray70”)
panel.points(x, dfrm$prop.predit, col= »gray70″, })
})
3) Une enquête cas-témoins a porté sur la relation entre la consommation d’alcool et celle du tabac et du cancer de l’œsophage chez l’homme (étude « Ille et Villaine”). Le groupe de cas comprenait 200 patients atteints d’un cancer de l’œsophage, diagnostiqués entre janvier 1972 et avril 1974. Au total, 775 hommes témoins ont été sélectionnés sur les listes électorales. Tableau 6.3 montre la répartition de la totalité des sujets en fonction de leur consommation quotidienne d’alcool, en considérant qu’une consommation supérieure à 80 g est considérée comme un facteur de risque.
Tableau 6.3. Consommation d’alcool et cancer de l’œsophage
| Consommation d’alcool (g/ jour) | |||
|---|---|---|---|
| ≥ 80 | & lt; 80 | Total | |
| Cas | 96 | 104 | 200 |
| Contrôle | 109 | 666 | 775 |
| Total | 205 | 770 | 975 |
a)
Quelle est la valeur du rapport de cotes et de son intervalle de confiance à 95% (méthode Woolf) ? Est-ce une bonne estimation du risque relatif?
b)
La proportion de consommateurs à risque est-elle la même parmi les cas et parmi les témoins (considérons α = 0,05)?
c)
Construire le modèle de régression logistique permettant de tester l’association entre la consommation d’alcool et le statut des individus. Le coefficient de régression est-il significatif ?
d)
Trouvez la valeur du rapport de cotes observé, calculé en (b) et son intervalle de confiance sur la base des résultats de l’analyse de régression.
Comme les données brutes (individuelles) ne sont pas disponibles, il est nécessaire de travailler directement avec la table de fréquences fournie dans l’énoncé du problème:
> alcohol < – matrix(c(666,104,109,96), nr = 2,
dimnames = list(c(« Control”, »Case”),
c(« <> = 80”)))
> alcohol
Regarding the odds ratio, the command oddsratio() from the package vcd will be used:
>library(vcd)
>oddsratio(alcohol, log=FALSE)
L’option log=FALSE garantit que le résultat retourné correspond correctement à un rapport de cotes et non au journal du rapport de cotes. Pour obtenir un intervalle de confiance asymptotique, on utilisera :
>confint(oddsratio(alcohol, log= FALSE))
De même, summary(oddsratio(alcohol)) pourrait être utilisé pour effectuer un test d’hypothèse sur le rapport de cotes log(H0: log(θ) = 0).
L’accessoire de commande.test() peut être utilisé pour tester l’hypothèse que la proportion de personnes ayant une consommation quotidienne est ≥ 80 g et est identique dans les cas ainsi que dans les témoins, indiquant les valeurs observées à partir de la tabulation croisée donnée dans l’énoncé du problème:
>prop.test (c(96,109), c(200,775), correct = FALSE)
Ce test est exactement équivalent à un test Z pour tester la différence entre deux proportions estimées à partir des données (si aucune correction de continuité n’est utilisée).
En ce qui concerne le modèle de régression logistique, le tableau de contingence doit être transformé en un tableau de données dans lequel les deux variables qualitatives (maladie et exposition) sont représentées correctement, c’est-à-dire avec une trame de données. Cela peut être fait en utilisant la commande melt() du package reshape2 qui permet de transformer les données sous forme de tableau en format « long”. A ce stade, les niveaux de la maladie variable doivent être recodés avec 0 (témoin) et 1 (cas), bien que cela ne soit pas vraiment nécessaire et que le niveau 0 soit considéré comme la catégorie de référence (ce qui facilite l’interprétation des résultats):
>bibliothèque (reshape2)
>alcool.df < – melt(alcool)
>noms (alcool.df) <-c(”maladie », ”exposition », ”n »)
> niveaux (alcool.dfdiseasemaladie) <-c(1,0)
> alcool.dfdiseasedisease < – relevel(alcool.dfdiseasedisease, « 0 »)
Le modèle de régression logistique s’écrit alors de la manière suivante :
>glm.res < – glm (exposition à la maladie, données = alcool.df,
famille=binomial, poids=n)
> résumé (glm.res)
Le résultat d’intérêt est la ligne associée à l’exposition >=80 car elle fournit des informations sur la valeur du coefficient de régression associé à la variable d’exposition estimée par R, avec son erreur type, ainsi que la valeur de la statistique de test. Ici, le coefficient de régression est interprété comme le journal du rapport de cotes. Notez que l’on obtiendrait exactement les mêmes résultats en échangeant le rôle des variables dans la formulation précédente, maladie d’exposition.
Le rapport de cotes calculé ci-dessus peut être trouvé à partir du coefficient de régression associé au facteur d’intérêt (exposition), ainsi que de son intervalle de confiance à 95%, de la manière suivante:
>exp(coef(glm.res))
>exp(confint(glm.res))
La deuxième solution consiste à considérer le nombre de cas et le nombre total d’individus:
>alcohol2 < – données.si vous avez besoin d’un cadre, vous pouvez le faire en utilisant la méthode suivante : <=80″), case=c(104,96),
total=c(104 + 666, 96 + 109))
> résumé (glm(cbind(case, total-case) expos, data=alcohol2,
famille=binomial))
4) Ces données proviennent d’une étude visant à établir la validité pronostique de la concentration de créatine kinase dans l’organisme sur la prévention de la survenue d’un infarctus du myocarde.
Les données sont disponibles dans le fichier sck.dat: la première colonne correspond à la créatine kinase variable (ck), la seconde à la présence de la maladie variable (pres) et la dernière à l’absence variable de maladie (abs).
a)
Quel est le nombre total de sujets ?
b)
Calculez les fréquences relatives malades/saines et représentez leur évolution en fonction des valeurs de la créatine kinase à l’aide d’un nuage de points (points + segments reliant les points).
c)
Sur la base d’un modèle de régression logistique qui vise à prédire la probabilité de tomber malade, calculez la valeur de CK à partir de laquelle ce modèle prédit que les personnes souffrent de la maladie en considérant un seuil de 0,5 (si P (malade) ≥ 0,5 alors malade = 1).
d)
Représentent graphiquement les probabilités d’être mal prédites par ce modèle ainsi que les proportions empiriques selon les valeurs ck.
e)
Établissez la courbe ROC correspondante et rapportez la valeur de l’aire sous la courbe.
Comme le nom des variables n’est pas dans le fichier de données, il sera nécessaire de les attribuer immédiatement après l’importation des données :
>sck < – read.tableau (« sck.dat”, header=FALSE)
>noms(sck)<-c(« ck”, »pres”, »abs”)
>résumé(sck)
Le total le nombre de sujets correspond à la somme des fréquences pour les deux variables pres et abs, c’est-à-dire :
>sum(sck)
ou, de manière équivalente: sum (sck$pres) + sum(sckabsabs) (mais pas sckprepres + sckabsabs).
Pour calculer les fréquences relatives de ces deux variables, il est nécessaire de connaître les grandeurs totales par variable. Celles-ci peuvent être obtenues en utilisant la commande apply et en opérant par colonnes :
>ni <- apply(sck, 2, sum)
Sur cette base, il suffit de diviser chaque valeur des variables pres et abs par les valeurs calculées précédemment. Les valeurs obtenues seront stockées dans deux nouvelles variables dans la même table de données :
>sck$pres.je ne sais pas si vous avez besoin de l’aide d’un fournisseur de services de sécurité, mais je ne sais pas si vous avez besoin de l’aide d’un fournisseur de services de sécurité.prop < -sck$abs/ni
Il est possible de vérifier que les calculs sont corrects : la somme des valeurs pour chaque variable doit maintenant être égale à 1 :
>apply(sck, 2, sum)
Alors les proportions obtenues sont simplement représentées dans un seul graphique, en considérant les valeurs de la variable ck en abscisses :
>xyplot(pres.accessoire + abs.il est possible de créer un fichier de données de type c(« b”, »g”),
auto.key=TRUE, ylab= »Frequency”)
L’instruction type=c(« b”, « g”) signifie que nous voulons que les points soient affichés reliés par des lignes (« b” = « o” + « l’) ainsi qu’une grille (« g”).
Le modèle de régression pour les données groupées est écrit comme suit :
>glm.res< – glm(cbind(pres, abs) ck, data=sck, family=binomial)
>résumé (glm.res)
Les prédictions, exprimées sous forme de probabilités et non sur l’échelle log-odds, sont obtenues au moyen de la commande predict en spécifiant l’option type= »response” comme suit :
>glm.pred < – predict(glm.res, type= »response »)
> noms (glm.pred)<-sck$ck
En considérant que les probabilités ≥ 0,5 désignent des individus malades, on obtient ainsi la distribution suivante :
>glm.pred
Il est conclu que les personnes seront considérées comme malades, selon ce modèle, pour des valeurs ck de 80 ou plus. Ceci est facilement vérifié avec un graphique où les probabilités prédites en fonction des valeurs de la variable ck sont représentées:
>sck$sick < – sckprepres/(sckprepres + sckabsabs)
> div> xyplot(glm.pred ~sck$ck, type= »l »,
ylab= »Probabilité”, xlab= »ck”,
panel=function( function) {
panel.xyplot( panel)
panneau.xyplot(sck$ck, scksicksick, pch=19, col= »grey »)
})
Dans un premier temps, il est nécessaire de ”décompresser » les données groupées et de créer un tableau avec deux colonnes : la première représentant la variable ck et la seconde représentant la présence ou l’absence de maladie. Nous utiliserons les comptages par sous-groupe précédemment calculés et disponibles dans la variable ni :
>sck.développez < – données.si vous avez besoin d’une table de rep (sck$ck, sckprepres),
rep(sck$ck, sckabsabs)),
sick =c(rep(1, ni), rep(0, ni)))
> table (sck.je n’ai pas besoin de le faire, mais je ne peux pas le faire.expand, tapply(sick, ck, sum))
Les deux dernières commandes sont destinées à s’assurer que les mêmes comptages sont rencontrés et que la répartition des personnes malades en fonction de la valeur de ck est correcte. On peut également vérifier que les mêmes résultats sont obtenus concernant la régression logistique et en même temps les valeurs prédites dans le tableau de données précédent peuvent être ajoutées (la même procédure que précédemment pourrait être suivie et les prédictions pourraient être répliquées, mais c’est plus simple de cette façon):
>glm.res2 < – glm(ck malade, data=sck.si vous avez besoin d’un nom de famille ou d’un nom de famille, n’hésitez pas à nous contacter.expandpredictionprediction < – ifelse(predict(glm.res2,
type= »réponse ») >=0.5,
1, 0)
> avec (sck.expand, table(sick, prediction))
La dernière commande permet d’afficher une matrice de confusion dans laquelle les diagnostics réels et ceux prédits par le modèle de régression sont croisés. Les taux de classification corrects peuvent être comparés lorsque nous faisons varier le seuil de référence :
>classif.tab < – avec (sck.expand, table(sick, prediction))
>sum(diag(classif.tab)) / sum(classif.pour afficher la courbe ROC, nous utiliserons le paquet ROCR :
>library(ROCR)
>pred < – prediction(predict(glm.res2, type= »réponse »),
sck.expand$sick)
> perf < – performance(pred, « tpr”, »fpr”)
> plot(perf, ylab = « Sensitivity”, xlab = « 1-Specificity”)
> grid()
> abline(0, 1, lty = 2)
The value of the area under the curve is obtained as follows:
> performance(pred, « auc”)@ »y.values”