Calcul d'un score après une regression logistique [Résolu]

[pistol22]

Bonjour,

Je cherche à automatiser le calcul d'un score après une régression logistique.
Dans un 1er temps je récupère les estimations des paramètres sur un échantillon d'apprentissage.
Ensuite je veux calculer un score sur un jeu de données test à partir des estimations de l’échantillon d"apprentissage.

Sur un exemple, voilà comment je procède actuellement :

Code :

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
DATA covar;
input aaa bbb ccc ;
cards;
0 0 1
1 1 1
1 1 0
0 1 1
1 1 1
0 0 2
1 0 0
1 0 1
1 1 0
0 1 2
run;
 
DATA estim;
input name$ estim;
cards;
aaa1 0.879
bbb1 0.578
ccc1 0.937
ccc2 0.376
run;
 
/* Création macro-variables avec l'estimation associée à chq variable */
DATA _null_;
SET estim;
call symputx(name , estim);
run;
 
/* Attribution des valeurs des estimations */
DATA tmp;
SET covar;
IF aaa=0 then V_1=0; IF aaa=1 then V_1=&aaa1;
IF bbb=0 then V_2=0; IF bbb=1 then V_2=&bbb1;
IF ccc=0 then V_3=0; IF ccc=1 then V_3=&ccc1; IF ccc=2 then V_3=&ccc2;
run;
 
/* Calcul d'un score pr chq individu */
DATA score;
SET tmp;
ARRAY V_col[3] V_1-V_3 ;
score=0;
DO i = 1 TO 3 ;
	score=score+V_col[i];
end;
DROP i;
RUN;

Le nombre de variables et le nombre de modalités par variables étant variables, je cherche donc à automatiser cette procédure. Est-ce que vous auriez quelques pistes à me conseiller ? Merci d'avance !

[MEGAMIND2]

L'instruction outmodel te permet de valider le modèle sur un autre èchantillon que celui qui a servi à construire le modèle. C'est plus simple, cela t'évite de passer par des macros.

[olivier.decourt]

Il y a (au moins) 2 solutions en plus de OUTMODEL/INMODEL proposé par Brice :
1) créer une table en empilant les données qui ont servi à construire le modèle, et les données à scorer. Pour ces dernières, remplacer les valeurs de Y par des valeurs manquantes. Mettre cette table en entrée de la proc LOGISTIC et récupérer les scores avec OUTPUT OUT. C'est bête et méchant, mais ça marche sur toutes les versions de SAS.
2) à partir de SAS 9, il y a une instruction SCORE dans la proc LOGISTIC qui permet d'indiquer une table à scorer (séparément de DATA qui reste la table pour construire le modèle).

Et puis si tu es très courageux (courageuse), il y a la récupération des coefficients, suivie d'une proc SCORE ou d'une macro pour faire la combinaison linéaire, puis d'une étape DATA pour appliquer le lien logit inverse.

Bref, les solutions de manquent pas, à toi de voir ce qui est le plus simple selon ta version de SAS (a priori, comme Brice, INMODEL/OUTMODEL).

[MEGAMIND2]

Si tu as la 9.1, voici comment tu dois procédé

Code :

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
DATA apprent;
input Y x;
cards;
1 415
0 48 
0 5
1 89
0 78
0 69
1 104
1 36
;
run;
 
DATA valid;
input Y x;
cards;
1 4
0 4 
0 50
1 9
0 115
0 899
1 904
1 31
;
run;
 
proc logistic DATA=apprent;
model Y(event="1") =x;
score DATA=apprent out=apprent_score;
score DATA=valid out=valid_score;
run;
 
/* cette macro te permet de calculer l'aire sous la courbe sur les données de validation à comparer avec avec celle d'apprentissage pour estimer la robustesse du modèle*/
 
%MACRO auc (TABLE=, variableScore=, variableY=, out=, event=) ;
	PROC SQL ;
		%IF &out NE %THEN %DO ;
			CREATE TABLE &out AS
		%END ;
		SELECT (SUM(t1.&variableScore > t0.&variableScore) 
                 + 0.5*SUM(t1.&variableScore = t0.&variableScore))
               / COUNT(*) AS aire FORMAT = 5.3
		FROM &TABLE (WHERE = (&variableY NE &event)) AS t0,
			 &TABLE (WHERE = (&variableY =  &event)) AS t1
		;
	QUIT ;
%MEND auc ;
 
%auc(TABLE=valid_score, variableScore=P_1, variableY=Y, out=aire_courbe_valid, event=1);
%auc(TABLE=apprent_score, variableScore=P_1, variableY=Y, out=aire_courbe_apprent, event=1);

[pistol22]

Bonjour,

Tout d'abord désolé pour le temps de réponse... Je vous remercie de vos réponses. Je ne connaissais pas les instructions outmodel/inmodel et c'est en effet beaucoup plus simple !

Cela dit, avec cette syntaxe je bute encore sur quelques points.

1- Si je veux construire une échelle de score, avec des sous-scores "ronds" pour chacune des modalités (et ainsi faciliter le calcul d'un score total pour chaque individu), il existe plusieurs méthodes que j'applique comme suit :

Code :

1
2
3
4
5
6
7
8
9
/* Multiplication par 4 du coefficient et arrondi à l'entier le + proche */
DATA estim;
SET estim;
Y=round(4*Y,1);
run;
/* OU division par le + petit coeff et arrondi à l'entier le + proche */
proc sql;
CREATE TABLE estim AS SELECT _NAME_, _LABEL_, round(Y/min(Y),1) AS devenir FROM estim;
quit;

Pour effectuer ces opérations, je ne vois pas d'autre moyen que de passer par le code que j'ai mis dans mon premier post, d'où la question de l'automatisation...


2- Je cherche à calculer un facteur de correction de mes coefficients de régression (comme la qualité apparente d’un modèle est généralement meilleure dans l’échantillon ayant servi à la construction du modèle que dans un autre échantillon), via le bootstrap.
Pour ce faire, j'applique la méthode suivante :

Envoyé par http://www.yaroslavvb.com/papers/steyerberg-application.pdf

1. Take a random bootstrap sample from the original sample, with the same size and
patient records drawn with replacement (see e.g. EFRON, 1993).
2. Estimate the logistic regression coefficients in the bootstrap sample.
3. Calculate the PI for each patient in the original sample. The PI is the linear
combination of the regression coefficients as estimated in the bootstrap sample
with the values of the covariables in the original sample.
4. Estimate the slope of the PI with logistic regression, using the outcomes of the
patients in the original sample and the PI as a single covariable.

J'ai donc besoin de récupérer la "pente" du score dans un modèle de régression logistique du résultat (variable Y) des individus de l’échantillon d’origine en fonction du score. Ce que je fais via la syntaxe suivante (à la suite du code de mon premier post) :

Code :

1
2
3
4
proc logistic DATA=valid outest=slope;
model Y(event="1") = score ;
score DATA=valid out=valid_score;
run;

Est-il possible de récupérer ce paramètre via la syntaxe inmodel/outmodel ?

C'est un peu compliqué, j'espère que j'ai été assez clair... Merci

[pistol22]

Si je reprends sur un exemple, voilà où j'en suis pour le moment.

Code :

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 
DATA valid;
input aaa bbb ccc ;
cards;
0 0 1
1 1 1
1 1 0
0 1 1
1 1 1
0 0 2
1 0 0
1 0 1
1 1 0
0 1 2
run;
 
DATA estim_apprent;
input name$ estim;
cards;
aaa1 0.879
bbb1 0.578
ccc1 1.235
ccc2 0.376
run;
 
DATA estim_apprent2;
SET estim_apprent;
estim=round(4*estim,1);
run;
 
/* Création macro-variables avec l'estimation associée à chq variable */
DATA _null_;
SET estim_apprent2;
call symputx(name , estim);
run;
%put _user_;
 
%let liste_var=aaa bbb ccc;
%LET nb_var=%EVAL(%SYSFUNC(count(%cmpres(&liste.),%STR( )))+1);
 
%macro score();
DATA tmp1;
SET valid;
%do i = 1 %TO &nb_var;
	%let liste = &liste;
	%let sous_score&i = %scan(&liste,&i);
	%let estim&i=&%scan(&liste,&i)1;
		IF &&sous_score&i=0 then sous_score&i=0 ;
		else sous_score&i=&&estim&i;
%end;
run;
 
%end;
run;
 
DATA tmp2;
SET tmp1;
ARRAY sum[&nb_var] sous_score1-sous_score&nb_var ;
score=0;
%DO i = 1 %TO &nb ;
	score=score+sum[&i];
%end;
RUN;
 
%mend score;
 
%score();

Je n'arrive pas à voir comment gérer le cas où une variable a plus de 2 modalités et donc plusieurs coefficients de régression associés.
Sur l'exemple, la variable ccc a 3 modalités : 0,1,2.
Si ccc=0 -> sous_score = 0 c'est ok
Si ccc=1 -> sous_score = &ccc1 = 5 ok
Si ccc=2 -> sous score = &ccc1 = 5 ce n'est pas bon, il faudrait &ccc2

Si quelqu'un voit une solution, merci d'avance !

[MEGAMIND2]

Je n'ai pas trop compris ce que tu cherches à faire mais je te propose ça pour automatiser avant d'avoir plus de détails

Code :

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
DATA valid;
input aaa bbb ccc ;
cards;
0 0 1
1 1 1
1 1 0
0 1 1
1 1 1
0 0 2
1 0 0
1 0 1
1 1 0
0 1 2
run;
 
DATA estim_apprent;
input name$ estim;
cards;
aaa1 0.879
bbb1 0.578
ccc1 1.235
ccc2 0.376
run;
 
DATA estim_apprent2;
SET estim_apprent;
estim=round(4*estim,1);
run;
 
/* Création macro-variables avec l'estimation associée à chq variable */
 
proc contents DATA=valid out=out;
run;
 
proc sort DATA=out;
BY varnum;
run;
proc sql;
SELECT name INTO: liste separated BY " " FROM out;
SELECT estim INTO: estim separated BY " " FROM estim_apprent2;
SELECT count(*) INTO:nb_var FROM out;
quit;
 
 
%put liste=&liste.;
%put estim =&estim .;
%put nb_var=&nb_var.;

[pistol22]

Je cherche à calculer un score pour chaque individu de ma table "valid", qui est la combinaison des coefficients de régression estimés dans la table "estim_apprent2" selon la valeur des covariables de la table "valid".

Plus clairement, pour l'individu 10 par exemple, on a :
aaa=0, bbb=1, ccc=2
Sachant que les estimations des paramètres sont :
aaa1=4, bbb1=2, ccc1=5, ccc2=2
J'obtiens un score pour cet individu égal à :
(0*aaa1 + bbb1 + ccc2) = 0+2+2 = 4

La table "estim_apprent" est une sortie d'une régression logistique. Dans cette table les noms des estimations des coefficients sont de la forme nom_variablenum_modalite, soit aaa1 pour la variable aaa qui a 2 modalités 0 (classe de référence) ou 1, et ccc1 et ccc2 pour la variable ccc qui a 3 modalités 0,1 et 2.

Le code dans mon précédent message permet de calculer correctement le score pour chaque individu dans le cas où je n'ai que des variables à 2 modalités, donc des noms des estimations des coefficients de la forme "nom_var1", d'où la ligne de code

Code :

%let estim&i=&%scan(&liste,&i)1;

Je coince dans le cas où j'ai des variables à j modalités (j>2), donc des noms des estimations des coefficients de la forme "nom_varj".

J'espère que je suis un peu plus clair...

[MEGAMIND2]

cet exemple ne ressemble pas à ce que tu souhaites obtenir?

Code :

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
26
27
28
29
 Data Neuralgia;
      input Treatment $ Sex $ Age Duration Pain $ @@;
      datalines;
   P  F  68   1  No   B  M  74  16  No  P  F  67  30  No
   P  M  66  26  Yes  B  F  67  28  No  B  F  77  16  No
   A  F  71  12  No   B  F  72  50  No  B  F  76   9  Yes
   A  M  71  17  Yes  A  F  63  27  No  A  F  69  18  Yes
   B  F  66  12  No   A  M  62  42  No  P  F  64   1  Yes
   A  F  64  17  No   P  M  74   4  No  A  F  72  25  No
   P  M  70   1  Yes  B  M  66  19  No  B  M  59  29  No
   A  F  64  30  No   A  M  70  28  No  A  M  69   1  No
   B  F  78   1  No   P  M  83   1  Yes B  F  69  42  No
   B  M  75  30  Yes  P  M  77  29  Yes P  F  79  20  Yes
   A  M  70  12  No   A  F  69  12  No  B  F  65  14  No
   B  M  70   1  No   B  M  67  23  No  A  M  76  25  Yes
   P  M  78  12  Yes  B  M  77   1  Yes B  F  69  24  No
   P  M  66   4  Yes  P  F  65  29  No  P  M  60  26  Yes
   A  M  78  15  Yes  B  M  75  21  Yes A  F  67  11  No
   P  F  72  27  No   P  F  70  13  Yes A  M  75   6  Yes
   B  F  65   7  No   P  F  68  27  Yes P  M  68  11  Yes
   P  M  67  17  Yes  B  M  70  22  No  A  M  65  15  No
   P  F  67   1  Yes  A  M  67  10  No  P  F  72  11  Yes
   A  F  74   1  No   B  M  80  21  Yes A  F  69   3  No
   ;
proc logistic data=Neuralgia outest=estim;
      class Treatment Sex;
      model Pain= Treatment Sex Treatment*Sex Age Duration ;
	output out=pred p=phat ;
   run;

[pistol22]

Oui, la probabilité estimée "phat" de la table "pred" est calculée à partir des estimations de la table "estim" suivant les covariables de la table "Neuralgia" (suivant la formule phat=1/[1+exp-(beta0 + beta1 * TreatmentA + beta2 * TreatmentB + beta3 * SexeF + ...)])

C'est aussi ce que je cherche à faire, mais "à la main" et le calcul que je fais est un peu plus simple :
score = (beta1 * TreatmentA + beta2 * TreatmentB + beta3 * SexeF + ...)

[MEGAMIND2]

Il te faut un modèle sans la constante:

Code :

1
2
3
4
5
6
7
8
9
10
proc logistic data=Neuralgia outest=estim;
      class Treatment Sex;
      model Pain= Treatment Sex Treatment*Sex Age Duration/NOINT ;
	output out=pred p=phat ;
run;

data pred;
	set pred;
	score= log((1-phat)/phat);
run;

[pistol22]

Merci pour la réponse, j'avais déjà pensé passer via cette solution (puis je ne sais plus pourquoi j'avais du bloquer et je me suis embarqué dans les macros...).
Par contre les estimations des paramètres ne sont pas les mêmes entre les modèles avec et sans constante, hors je veux calculer mon score à partir des estimations des coefficients du modèle AVEC constante.

Si je ne me trompe pas, ça donne donc ça :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
proc logistic DATA=Neuralgia outest=estim;
      class Treatment Sex;
      model Pain= Treatment Sex Treatment*Sex Age Duration ;
	output out=pred p=phat ;
run;
 
proc sql noprint;
SELECT Intercept INTO: int FROM estim;
quit;
 
DATA pred;
	SET pred;
	score= -log((1-phat)/phat)-∫
run;

Ça devrait pouvoir s'adapter à ce que je veux faire par la suite. Merci bien !

[MEGAMIND2]

J'y ai pensé après avoir posté, par contre il faudrait que tu enlèves l'option NOINT afin d'avoir la constante. Et pourquoi ya un moins devant la log?

[pistol22]

Ok pour le NOINT (j'ai édité le message précédent).

Pour le moins devant le log, si on cherche à isoler le score à partir de l'équation :
phat=1/[1+exp-(beta0 + beta1 * X1)] avec score = (beta1 * X1)

Ca donne bien :

score = -log((1-phat)/phat)-beta0, non ?

[MEGAMIND2]

Au tant pour moi, j'ai omis le moins de l'exponentiel

[pistol22]

En fait je crois que cette solution ne résout pas vraiment mon problème... Je m'explique :
Avec cette solution, on arrive à calculer le score voulu (beta1 * X1 + beta2 * X2 + ... ) pour chacun des individus de l'échantillon d'apprentissage.
(autrement dit : le score est calculé à partir des estimations de l'ech d'apprentissage suivant les covariables l'ech d'apprentissage)

Hors je voudrais calculer ce score pour les individus d'un autre échantillon (éch de validation).
(autrement dit : le score doit être calculé à partir des estimations de l'ech d'apprentissage suivant les covariables l'ech de validation).

[MEGAMIND2]

Ce sont les estimations calculées sur l'échantillon d'apprentissage qui est appliqué sur l'échantillon de validation, ça ne devait pas poser de problèmes.

Code :

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
26
 
 
 
%let seed=123;
DATA apprent valid;
SET neuralgia;
IF ranuni(&seed.)<0.66 then output apprent;
else output valid;
run;
 
proc logistic DATA=apprent outest=estim;
class Treatment Sex;
model Pain(event="yes")= Treatment Sex Treatment*Sex Age Duration ;
score DATA=apprent out=apprent_score;
score DATA=valid out=valid_score;
run;
 
proc sql noprint;
SELECT Intercept INTO: int FROM estim;
quit;
 %put intercept=&int.;
 
DATA valid_score;
	SET valid_score;
	score= -log((1-P_yes)/P_yes)-&int;
run;

[pistol22]

Ok c'est parfait, merci !