Discussion :

[Emmanuel R]

Bonjour,
_____

J'ai déjà posté ce sujet sur le forum CIRAD, sans succès jusque là. J'espère que vous pourrez m'aider. Si j'ai une réponse sur un des deux sujets je ferai part de la réponse dans tous les cas.
Merci par avance pour votre aide.
_____

Je suis en train de travailler sur une régression logistique et j'ai détecté une anomalie lorsque je fais une validation croisée. Je poste donc ici dans l'espoir que vous m'aidiez à cerner le problème.

J'utilise 100 observations de deux variables générées aléatoirement selon des lois normales de paramètres différents.
Par la suite je crée un modèle qui a l'air assez cohérent (voir la courbe ROC à la fin).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
#### Génération des données ####
z = c(rep(0,50),rep(1,50))
x = c(rnorm(50,0,4),rnorm(50,2,4))
y = c(rnorm(50,0,1),rnorm(50,2,1))
plot(x,y,col=z+1)
 
#### Création du modèle et prédiction par GLM ###
res.glm = glm(formula=z~x+y,family=binomial)
predict.res.glm = predict(res.glm, type="response")

J'effectue maintenant la validation croisée sur quatre groupes. Le principe étant de prédire les données de chaque groupe à partir des trois autres.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#### Validation croisée ####
 
# Définition de 4 groupes. Trois d'entres eux construisent le modèle, on prédit sur le dernier.
g = sample(rep(1:4,25))
 
# Initialisation des prédictions
vect.proba = rep(0,100)
 
for (k in 1:4){
  # Echantillon d'apprentissage sur 3 groupes
  appr = as.data.frame(cbind(x[g!=k], y[g!=k]))
  z.appr = z[g!=k]
  res.glm.kfold = glm(formula=z.appr~appr[,1]+appr[,2],family=binomial)
 
  # Echantillon de vérification
  verif = as.data.frame(cbind(x[g==k], y[g==k]))
  vect.proba[g==k] = predict(res.glm.kfold, newdata = verif, type="response") #https://stat.ethz.ch/R-manual/R-patched/library/stats/html/predict.lm.html
}

Alors ici, plusieurs choses :
--> J'obtiens des warnings (8 exactement) au niveau de la prédiction dans la boucle.

1: 'newdata' had 25 rows but variables found have 75 rows
2: In vect.proba[g == k] = predict(res.glm.kfold, newdata = verif, :
number of items to replace is not a multiple of replacement length

J'avoue ne pas comprendre car la définition de newdata est, selon le manuel d'utilisation :

optionally, a data frame in which to look for variables with which to predict. If omitted, the fitted linear predictors are used.

Je fournis bien le même nombre de variables, mais en effet il n'y a pas le même nombre d'observations ce qui est tout à fait normal. J'ai vu plusieurs scripts sur le net qui utilisent ce paramètre pour faire une validation croisée (ici par exemple).
D'où la question : Est-ce que c'est la bonne chose à faire ?

--> Ensuite : mes prédictions par validation croisée sont foireuses. Je vous indique un extrait des prédictions entre la 40e et la 60e ligne (à la 50e, on passe de 0 à 1)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Nb  Z     Predic_GLM  Predic_ValCroisee
40  0     0.319408438 0.145276517
41  0     0.205839858 0.119031341
42  0     0.012653655 0.029117270
43  0     0.181420164 0.749168854
44  0     0.020543180 0.914641430
45  0     0.442173491 0.645873783
46  0     0.157567167 0.131166378
47  0     0.045872366 0.102201020
48  0     0.085971756 0.060728705
49  0     0.489664348 0.021215917
50  0     0.013036409 0.036993982
51  1     0.984022606 0.007688740
52  1     0.918224041 0.039215022
53  1     0.944663266 0.201269337
54  1     0.297424130 0.701688379
55  1     0.803391279 0.897638029
56  1     0.946804556 0.210754072
57  1     0.260940236 0.097831911
58  1     0.996113114 0.130593783
59  1     0.919221343 0.071754892
60  1     0.508894219 0.071091520

Vous pouvez voir que dans certains cas les prédictions ne concordent pas du tout. J'ai testé les différents modèles en modifiant les valeurs de k à la main pour tester. Les paramètres ne sont pas excessivement différents, et du coup je ne comprends pas pourquoi il y a un tel écart.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
> res.glm
Coefficients:
(Intercept)            x            y  
    -2.0162       0.1303       2.1463  
 
> glm(formula=z.appr.test.1~appr.test.1[,1]+appr.test.1[,2],family=binomial)
Coefficients:
     (Intercept)  appr.test.1[, 1]  appr.test.1[, 2]  
          -2.121             0.224             2.437  
 
> glm(formula=z.appr.test.2~appr.test.2[,1]+appr.test.2[,2],family=binomial)
Coefficients:
     (Intercept)  appr.test.2[, 1]  appr.test.2[, 2]  
         -2.3885            0.1578            2.3182  
 
 
> glm(formula=z.appr.test.3~appr.test.3[,1]+appr.test.3[,2],family=binomial)
Coefficients:
     (Intercept)  appr.test.3[, 1]  appr.test.3[, 2]  
        -2.02831           0.09401           2.26203  
 
> glm(formula=z.appr.test.4~appr.test.4[,1]+appr.test.4[,2],family=binomial)
Coefficients:
     (Intercept)  appr.test.4[, 1]  appr.test.4[, 2]  
        -1.73949           0.09701           1.78679

Je finis en indiquant les courbes ROC, où vous pouvez bien voir que les prédictions obtenues par ma validation croisée sont très mauvaises.



Le code pour info

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
library("ROCR")
par(mfrow=c(1,2))
 
roc.pred.glm <- prediction(predict.res.glm,z)
roc.perf.glm <- performance(roc.pred.glm,"tpr","fpr")
plot(roc.perf.glm,colorize=T,lwd=3,xlab = "Taux de faux positifs",ylab = "Taux de vrais positifs",main = "Courbe ROC - Prédiction par GLM")
 
roc.pred.croisee <- prediction(vect.proba,z)
roc.perf.croisee <- performance(roc.pred.croisee,"tpr","fpr")
plot(roc.perf.croisee,colorize=T,lwd=3,xlab = "Taux de faux positifs",ylab = "Taux de vrais positifs",main = "Courbe ROC - Validation croisée")

Est-ce que vous avez une idée du problème pour faire ma validation croisée de manière cohérente ?
J'espère que le pavé ne posera pas de souci, mais vu que je m'arrache les cheveux depuis quelques heures sur cela, je tenais à être le plus précis possible.

Merci par avance pour votre aide.

Cordialement,

Emmanuel Rambeau

[dev_ggy]

Bonjour,

Je découvre avec toi la validation croisée. Après avoir lu un doc sur le sujet. Je m'interroge. Pourquoi prends-tu le reste des données et ne travailles-tu pas sur les groupes générés. Pour faire plus simple pourquoi n’utilises-tu pas une différence au lieu d’une égalité dans ta validation croisée ?

Cordialement.

[Emmanuel R]

Bonjour,
L'idée de la validation croisée est de faire le modèle avec un nombre restreint de groupe (ici 3, mais ça peut être plus ou moins) et de valider sur le groupe restant. Ainsi, lorsqu'on fait la prédiction du dernier groupe, ce dernier n'a pas participé à la création du modèle. Ce dernier n'est donc pas influencé. Si les prédictions sont cohérentes, cela permet de valider le modèle.

Dans mon code, tu peux voir que j'utilise à tour de rôle l'échantillon d'apprentissage et l'échantillon de vérification.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
for (k in 1:4){
  # Echantillon d'apprentissage sur 3 groupes
  appr = as.data.frame(cbind(x[g!=k], y[g!=k]))
  z.appr = z[g!=k]
  res.glm.kfold = glm(formula=z.appr~appr[,1]+appr[,2],family=binomial)
  
  # Echantillon de vérification
  verif = as.data.frame(cbind(x[g==k], y[g==k]))
  vect.proba[g==k] = predict(res.glm.kfold, newdata = verif, type="response")
}

Je précise que le vect.proba est le vecteur où je range toutes les prédictions (une ligne par observation).

J'utilise bien donc une différence puis ensuite une égalité.
Était-ce bien ta question ?

Je reste dispo si besoin de plus d'infos.