J'ai un jeu de données contenant des mesures de concentrations en différentes molécules, ainsi qu'un statut Cas ou Témoin pour mes échantillons.
Afin de trouver le(s) meilleur(s) prédicteurs du statut Cas/Témoin, j'ai construit un modèle logistique comme suit :
Code : Tout sélectionner
#Aperçu des données
quantif[36:40,2:ncol(quantif)]
case Gluc Asc Gsh Lac GPC Gly Succ bHB Val Leu.Iso
36 0 4.063753 2.2016668 0.458457233 2.398714 1.236988112 11.445137 0.2232978 0.4887412 0.04124332 0.19572592
37 0 4.252199 2.2591980 0.563522342 2.330555 4.127798786 12.219406 0.3387899 0.5608153 0.15179670 0.03442021
38 0 4.696948 1.8506493 0.500991625 2.186188 4.383154903 11.681825 0.2496662 0.5350908 0.14228403 0.30912699
39 1 6.226774 0.8334419 0.009612858 22.858587 0.001401988 8.529733 7.2276462 1.8062908 0.83424979 1.18558873
40 1 10.223897 1.2005618 0.054661308 4.327522 0.023142013 8.550414 2.9679294 4.5692782 1.34797602 2.32337683
# Etape 1 : modèle vide
glm0 <- glm(case~1, data=quantif, family=binomial(link=logit))
# Etape 2 : Procédure "stepAIC" pour l'ajout des variables une à une en respectant le principe de Parcimonie
stepAIC(glm0, .~Gluc+Asc+Gsh+Lac+GPC+Gly+Succ+bHB+Val+Leu.Iso, data=quantif, family=binomial())
## Start: AIC=70.59
## case~ 1
##
## Df Deviance AIC
## + Succ 1 23.136 27.136
## + Gly 1 39.416 43.416
## + bHB 1 48.049 52.049
## + GPC 1 51.068 55.068
## + Val 1 53.598 57.598
## + Leu.Iso 1 55.970 59.970
## + Lac 1 56.849 60.849
## + Asc 1 62.614 66.614
## + Gsh 1 65.989 69.989
## <none> 68.593 70.593
## + Gluc 1 67.261 71.261
##
## Step: AIC=27.14
## case ~ Succ
##
## Df Deviance AIC
## + Gly 1 0.000 6.000
## + Lac 1 10.095 16.095
## + Gluc 1 15.700 21.700
## + Asc 1 18.379 24.379
## + GPC 1 20.189 26.189
## <none> 23.136 27.136
## + bHB 1 21.693 27.693
## + Leu.Iso 1 22.505 28.505
## + Val 1 22.855 28.855
## + Gsh 1 23.131 29.131
## - Succ 1 68.593 70.593
##
## Step: AIC=6
## case ~ Succ + Gly
##
## Df Deviance AIC
## <none> 0.000 6.000
## + GPC 1 0.000 8.000
## + Gsh 1 0.000 8.000
## + Lac 1 0.000 8.000
## + Asc 1 0.000 8.000
## + Val 1 0.000 8.000
## + Gluc 1 0.000 8.000
## + Leu.Iso 1 0.000 8.000
## + bHB 1 0.000 8.000
## - Gly 1 23.136 27.136
## - Succ 1 39.416 43.416
J'ai donc décidé à ce stade d'évaluer les deux modèles.
Modèle 1 : Succ + Gly
Code : Tout sélectionner
glmsuccgly <- glm(formula = case ~ Succ + Gly, family = binomial(link=logit), data = quantif)
summary(glmsuccgly)
Call:
glm(formula = case ~ Succ + Gly, family = binomial(link = logit),
data = quantif)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.194e-04 -2.100e-08 2.100e-08 2.100e-08 1.480e-04
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 419.04 99295.13 0.004 0.997
Succ 321.72 72991.84 0.004 0.996
Gly -58.72 13419.51 -0.004 0.997
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6.8593e+01 on 49 degrees of freedom
Residual deviance: 4.1559e-08 on 47 degrees of freedom
AIC: 6
Number of Fisher Scoring iterations: 25
# Prédictions vs. Observations :
quantif$pred_succgly <- predict(glmsuccgly,type="response", newdata=quantif)
quantif$pred_succgly01 <- ifelse(quantif$pred_succgly>=0.5, 1, 0)
tab2 <- table(quantif$case, quantif$pred_succgly01, dnn = c("Expérience", "Prédiction"))
rownames(tab2) <- ifelse(rownames(tab2)==0, "Témoin", "Cas")
colnames(tab2) <- ifelse(colnames(tab2)==0, "Témoin", "Cas")
tab2
Prédiction
Expérience Témoin Cas
Témoin 38 0
Cas 0 12
... étrange, pourquoi les p-value sont-elles si grandes ? Alors que tous mes échantillons sont bien classés ?
Je me suis dit que j'avais loupé quelque chose, et ai continué avec le modèle plus simple :
Modèle 2 : Succ
Code : Tout sélectionner
lmsucc <- glm(formula = case ~ Succ, family = binomial(link=logit), data = quantif)
summary(glmsucc)
Call:
glm(formula = case ~ Succ, family = binomial(link = logit),
data = quantif)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.15615 -0.30840 0.00006 0.07010 2.02193
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.275 1.681 -3.137 0.00171 **
Succ 8.025 2.745 2.924 0.00346 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 68.593 on 49 degrees of freedom
Residual deviance: 23.136 on 48 degrees of freedom
AIC: 27.136
Number of Fisher Scoring iterations: 8
OK ici je trouve ça bien mieux : ma variable Succ est significativement liée au modèle (mais du coup je comprends encore moins pourquoi ce n'était pas le cas dans le modèle 1...).
J'ai enfin une dernière question pour l'interprétation des coefficients :
Puis-je écrire que logit(case) = -5.275 + 8.025*Succ ?
J'espère avoir été claire dans mes questions,
Merci d'avance si quelqu'un peut m'aider !