Modèle de Point de Rupture (Change Point Model) avec PROC MCMC

Le script met en œuvre un modèle de point de rupture pour identifier un changement structurel dans la relation entre deux variables. Il commence par créer un jeu de données interne, puis utilise PROC MCMC pour estimer les paramètres d'un modèle linéaire par morceaux, y compris la position du point de rupture 'cp'. Les résultats de l'estimation (moyennes a posteriori) sont ensuite stockés dans des macro-variables via un DATA _NULL_ step. Finalement, PROC SGPLOT est utilisé deux fois : une première fois pour visualiser les données brutes, et une seconde fois pour superposer les données, les droites de régression estimées de part et d'autre du point de rupture, et la densité a posteriori de l'emplacement de ce point.

Analyse des données

Type : CREATION_INTERNE

Les données sont entièrement générées au sein du script à l'aide d'une étape DATA avec une instruction DATALINES. Aucune donnée externe n'est nécessaire.

1 Bloc de code

DATA STEP Data

Explication :
Ce bloc crée le jeu de données 'stagnant' à partir de données internes (datalines). L'option '@@' dans l'instruction INPUT permet de lire plusieurs observations sur une même ligne de données.

Copié !

1	title 'Change Point Model';
2	DATA stagnant;
3	INPUT y x @@;
4	ind = _n_;
5	DATALINES;
6	1.12 -1.39 1.12 -1.39 0.99 -1.08 1.03 -1.08
7	0.92 -0.94 0.90 -0.80 0.81 -0.63 0.83 -0.63
8	0.65 -0.25 0.67 -0.25 0.60 -0.12 0.59 -0.12
9	0.51 0.01 0.44 0.11 0.43 0.11 0.43 0.11
10	0.33 0.25 0.30 0.25 0.25 0.34 0.24 0.34
11	0.13 0.44 -0.01 0.59 -0.13 0.70 -0.14 0.70
12	-0.30 0.85 -0.33 0.85 -0.46 0.99 -0.43 0.99
13	-0.65 1.19
14	;
15	RUN;

2 Bloc de code

PROC SGPLOT

Explication :
Visualisation initiale des données brutes sous forme de nuage de points pour inspecter la relation entre x et y.

Copié !

1
2	PROC SGPLOT
3	DATA=stagnant;
4	scatter x=x y=y;
5	RUN;
6

3 Bloc de code

PROC MCMC

Explication :
Cœur de l'analyse. PROC MCMC est utilisée pour ajuster un modèle bayésien. Le code définit les paramètres (alpha, cp, beta1, beta2, s2), leurs distributions a priori, et le modèle linéaire qui stipule que la pente (beta) change en fonction de la position de x par rapport au point de rupture 'cp'. La procédure génère des échantillons de la distribution a posteriori des paramètres.

Copié !

1	PROC MCMC DATA=stagnant outpost=postout seed=24860 ntu=1000
2	nmc=20000;
3	ods select PostSumInt;
4	ods OUTPUT PostSumInt=ds;
5
6	array beta[2];
7	parms alpha cp beta1 beta2;
8	parms s2;
9
10	prior cp ~ unif(-1.3, 1.1);
11	prior s2 ~ uniform(0, 5);
12	prior alpha beta: ~ normal(0, v = 1e6);
13
14	j = 1 + (x >= cp);
15	mu = alpha + beta[j] * (x - cp);
16	model y ~ normal(mu, var=s2);
17	RUN;

4 Bloc de code

DATA STEP

Explication :
Ce DATA step de type _NULL_ (ne crée pas de table) lit les moyennes a posteriori calculées par PROC MCMC (stockées dans la table 'ds') et les assigne à des macro-variables (ex: &cp, &beta1, &beta2) pour une utilisation ultérieure.

Copié !

1	DATA _null_;
2	SET ds;
3	call symputx(parameter, mean);
4	RUN;

5 Bloc de code

DATA STEP Data

Explication :
Crée un jeu de données 'b' pour tracer les lignes de régression du modèle ajusté. Il utilise les macro-variables (&cp, &alpha, &beta1, &beta2) pour calculer les valeurs prédites de y1.

Copié !

1	DATA b;
2	missing A;
3	INPUT x1 @@;
4	IF x1 eq .A THEN x1 = &cp;
5	IF _n_ <= 2 THEN y1 = &alpha + &beta1 * (x1 - &cp);
6	ELSE y1 = &alpha + &beta2 * (x1 - &cp);
7	DATALINES;
8	-1.5 A 1.2
9	;
10	RUN;

6 Bloc de code

PROC KDE Data

Explication :
Utilise la procédure KDE (Kernel Density Estimation) sur la sortie de MCMC pour estimer la distribution de probabilité a posteriori du paramètre 'cp' (le point de rupture).

Copié !

1
2	PROC KDE
3	DATA=postout;
4	univar cp / out=m1 (drop=count);
5	RUN;
6

7 Bloc de code

DATA STEP Data

Explication :
Ce bloc ajuste (met à l'échelle et décale) la valeur de la densité calculée par PROC KDE afin qu'elle puisse être superposée de manière lisible sur le graphique final.

Copié !

1	DATA m1;
2	SET m1;
3	density = (density / 25) - 0.653;
4	RUN;

8 Bloc de code

DATA STEP Data

Explication :
Combine les données originales ('stagnant'), les lignes du modèle ajusté ('b') et les données de densité du point de rupture ('m1') en un seul jeu de données pour le graphique final.

Copié !

1	DATA all;
2	SET stagnant b m1;
3	RUN;

9 Bloc de code

PROC SGPLOT

Explication :
Génère le graphique final complet qui superpose : le nuage de points original, les deux segments de droite du modèle de point de rupture, et la courbe de densité de la distribution a posteriori du paramètre 'cp'.

Copié !

1	PROC SGPLOT DATA=all noautolegend;
2	scatter x=x y=y;
3	series x=x1 y=y1 / lineattrs = graphdata2;
4	series x=value y=density / lineattrs = graphdata1;
5	RUN;

Ce matériel est fourni "tel quel" par We Are Cas. Il n'y a aucune garantie, expresse ou implicite, quant à la qualité marchande ou à l'adéquation à un usage particulier concernant le matériel ou le code contenu dans les présentes. We Are Cas n'est pas responsable des erreurs dans ce matériel tel qu'il existe maintenant ou existera, et We Are Cas ne fournit pas de support technique pour celui-ci.

Informations de Copyright : S A S S A M P L E L I B R A R Y

Retour à la liste