Ce script modélise la croissance de la population américaine. Il démontre d'abord comment implémenter la méthode des moindres carrés repondérés itérativement (IRLS) en utilisant la fonction 'biweight' de Tukey directement dans PROC NLIN grâce à la variable spéciale `_weight_`. Il compare ensuite cette approche avec la procédure standardisée PROC ROBUSTREG utilisant un M-estimateur.
Analyse des données
Type : CREATION_INTERNE
Les données de population (uspop) sont créées directement dans le script via une étape DATA et des instructions DATALINES.
1 Bloc de code
DATA STEP Data
Explication : Création du jeu de données 'uspop'. L'instruction INPUT utilise '@@' pour lire plusieurs observations sur une même ligne de données. Les variables 'year' et 'yearsq' (année au carré) sont générées pour le modèle quadratique.
Copié !
title 'U.S. Population Growth';
data uspop;
input pop :6.3 @@;
retain year 1780;
year = year+10;
yearsq = year*year;
datalines;
3929 5308 7239 9638 12866 17069 23191 31443 39818 50155
62947 75994 91972 105710 122775 131669 151325 179323 203211
226542 248710
;
Explication : Exécution d'une régression non-linéaire. Le script implémente manuellement une régression robuste (IRLS) en définissant la variable spéciale `_weight_`. Les poids sont recalculés à chaque itération en fonction des résidus (`resid`) selon la fonction biweight de Tukey.
Copié !
title 'Beaton/Tukey Biweight Robust Regression using IRLS';
proc nlin data=uspop nohalve;
parms b0=20450.43 b1=-22.7806 b2=.0063456;
model pop=b0+b1*year+b2*year*year;
resid = pop-model.pop;
sigma = 2;
k = 4.685;
if abs(resid/sigma)<=k then _weight_=(1-(resid / (sigma*k))**2)**2;
else _weight_=0;
output out=c r=rbi;
run;
1
title 'Beaton/Tukey Biweight Robust Regression using IRLS';
2
PROC NLINDATA=uspop nohalve;
3
parms b0=20450.43 b1=-22.7806 b2=.0063456;
4
model pop=b0+b1*year+b2*year*year;
5
resid = pop-model.pop;
6
sigma = 2;
7
k = 4.685;
8
IF abs(resid/sigma)<=k THEN _weight_=(1-(resid / (sigma*k))**2)**2;
9
ELSE _weight_=0;
10
OUTPUT out=c r=rbi;
11
RUN;
3 Bloc de code
DATA STEP Data
Explication : Étape DATA pour recalculer ou vérifier les poids finaux basés sur les résidus (rbi) stockés dans la table de sortie de la PROC NLIN.
Copié !
data c;
set c;
sigma = 2;
k = 4.685;
if abs(rbi/sigma)<=k then _weight_=(1-(rbi /(sigma*k))**2)**2;
else _weight_=0;
run;
1
DATA c;
2
SET c;
3
sigma = 2;
4
k = 4.685;
5
IF abs(rbi/sigma)<=k THEN _weight_=(1-(rbi /(sigma*k))**2)**2;
6
ELSE _weight_=0;
7
RUN;
4 Bloc de code
PROC PRINT
Explication : Affichage des résultats de la régression manuelle (table 'c').
Copié !
proc print data=c;
run;
1
PROC PRINTDATA=c;
2
RUN;
5 Bloc de code
PROC ROBUSTREG Data
Explication : Utilisation de la procédure PROC ROBUSTREG dédiée à la régression robuste. Elle utilise ici la méthode d'estimation M. Les résultats sont stockés dans la table 'weights'.
Copié !
proc robustreg data=uspop method=m(scale=2);
model pop = year year*year;
output out=weights weight=w;
run;
1
PROC ROBUSTREGDATA=uspop method=m(scale=2);
2
model pop = year year*year;
3
OUTPUT out=weights weight=w;
4
RUN;
6 Bloc de code
PROC PRINT
Explication : Affichage des poids calculés par la procédure ROBUSTREG pour comparaison.
Copié !
proc print data=weights;
run;
1
PROC PRINTDATA=weights;
2
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.
SAS et tous les autres noms de produits ou de services de SAS Institute Inc. sont des marques déposées ou des marques de commerce de SAS Institute Inc. aux États-Unis et dans d'autres pays. ® indique un enregistrement aux États-Unis. WeAreCAS est un site communautaire indépendant et n'est pas affilié à SAS Institute Inc.
Ce site utilise des cookies techniques et analytiques pour améliorer votre expérience.
En savoir plus.