Bonjour,
J’ai une question qui ne rentre pas véritablement dans le cadre de la
librairie ade4: il s’agit de régression logistique.
L'exemple est tiré de "Applied Logistic Regression, 2nd Edition", David
W.Hosmer et Stanley Lemeshow.
Chapitre 8 : régression logistique multinomiale.
p. 264 : exemple d'une étude sur la mammographie.
Les données peuvent être obtenues à cette adresse:
http://www-unix.oit.umass.edu/~statdata/statdata/stat-logistic.html
Dans un premier temps, les auteurs considèrent un modèle contenant une
seule variable dichotomique codée 0 ou 1 (HIST).
La variable réponse est ME, codée en 0, 1 et 2.
ME=0 est utilisé comme référence.
J’ai utilisé la fonction multinom(nnet) afin de reproduire cet exemple et
ce modèle.
data <- read.delim2("D:/mes donnees/…./data.txt")
# tableau croisé :
table(data$ME, data$HIST)
# 0 1
# 0 220 14
# 1 85 19
# 2 63 11
library(nnet)
mult <- multinom(formula = ME ~ HIST, data = data)
summary(mult)
# Call:
# multinom(formula = ME ~ HIST, data = data)
# Coefficients:
# (Intercept) HIST1
# 1 -0.9509794 1.256531
# 2 -1.2504803 1.009384
# Std. Errors:
# (Intercept) HIST1
# 1 0.1277115 0.3746633
# 2 0.1428926 0.4275097
# Residual Deviance: 792.34
# AIC: 800.34
# Correlation of Coefficients:
# 1:(Intercept) 1:HIST1 2:(Intercept)
# 1:HIST1 -0.3408700
# 2:(Intercept) 0.2490796 -0.0849038
# 2:HIST1 -0.0832534 0.4743673 -0.3342442
J’obtiens les mêmes résultats que les auteurs, concernant les coefficients
et les écarts types, au millième près.
Le calcul des OR à partir de coefficients ne pose pas de problème:
OR1 correspond à ME=1 vs ME=0
OR2 correspond à ME=2 vs ME=0
OR1 <- exp(as.data.frame(summary(mult)$coefficients)$HIST[1])
OR2 <- exp(as.data.frame(summary(mult)$coefficients)$HIST[2])
> OR1
[1] 3.513213
> OR2
[1] 2.743911
Le calcul des intervalles de confiance à 95% non plus:
Int95_OR1p <-
exp((as.data.frame(summary(mult)$coefficients)$HIST[1])+1.96*(as.data.frame(summary(mult)$standard.errors)$HIST[1]))
Int95_OR1m <-
exp((as.data.frame(summary(mult)$coefficients)$HIST[1])-1.96*(as.data.frame(summary(mult)$standard.errors)$HIST[1]))
Int95_OR1 <- round(c(Int95_OR1m, Int95_OR1p), dig=3)
Int95_OR2p <-
exp((as.data.frame(summary(mult)$coefficients)$HIST[2])+1.96*(as.data.frame(summary(mult)$standard.errors)$HIST[2]))
Int95_OR2m <-
exp((as.data.frame(summary(mult)$coefficients)$HIST[2])-1.96*(as.data.frame(summary(mult)$standard.errors)$HIST[2]))
Int95_OR2 <- round(c(Int95_OR2m, Int95_OR2p), dig=3)
> Int95_OR1
[1] 1.686 7.322
> Int95_OR2
[1] 1.187 6.343
J’aimerais obtenir un OR de ME=2 versus ME=1.
Dans le livre ils calculent un coefficient (différence entre les 2
coefficients: 1.009 1.256 = -0.247 donc OR = 0.781).
Ensuite ils calculent un estimateur de la variance de la différence entre
les 2 coefficients, et ceci à partir d’une matrice de covariance estimée,
donnée directement par le logiciel utilisé (en l’occurrence STATA ;
multinom donne une matrice de corrélations). Cet estimateur de la variance
est ensuite utilisé pour calculer l’intervalle de confiance autour du
coefficient et donc de l’OR.
Refaire tout ces calculs ne pose évidemment pas de pb.
Cependant, j’aimerais savoir si l’on peut paramétrer la fonction multinom
afin de pouvoir choisir la classe de référence de la variable réponse, en
l’occurrence pour choisir ME=1 comme référence et obtenir directement un OR
de ME=2 vs ME=1?
Merci
Charline L.
This archive was generated by hypermail 2b30 : Fri Jan 21 2005 - 14:51:56 MET