[SAS 프로그래밍 실습] IML로 Weighted Regression 구현하기
※본 포스팅은 선형대수에 대한 선수지식을 요구하는 내용을 포함하고있습니다.
부제 : Weighted regression with IML
- Introduction
지난 시간, 우리는 University Edition을 통해, REG 프록시져를 사용하여 가중 회귀를 구현해보았습니다. (밑의 링크 참조) 간단하게 구현해보았지만, 회귀 모델의 세세한 구현 프로세스를 알기엔 아쉬운 절차였습니다. 이번 시간에는 IML을 통해 조금 더 자세하게 구현해보도록 하겠습니다.
https://blog.naver.com/statpark1014/221521344309
(PROC REG로 가중회귀 구현하기)
- Data
데이터는 이전 포스팅에서 사용한 데이터로 진행하겠습니다.
<code1>
data RegData;
input y x w;
datalines;
2.3 7.4 0.058
3.0 7.6 0.073
2.9 8.2 0.114
4.8 9.0 0.144
1.3 10.4 0.151
3.6 11.7 0.119
2.3 11.7 0.119
4.6 11.8 0.114
3.0 12.4 0.073
5.4 12.9 0.035
6.4 14.0 0
;
- Practice
IML을 사용해보도록 하겠습니다.
우리가 구현하고 싶은 식은 아래와 같습니다.
가중 선형회귀 이론에 대한 자세한 내용은 아래를 참고하시기 바랍니다.
https://blog.naver.com/statpark1014/221519820569
(가중 선형회귀 이론)
<code2>
proc iml;
/* IML을 통한 방정식 계산*/
start PolyRegEst(Y, X, w, deg);
Yw = sqrt(w)#Y; /* 1. w 제곱근에 Y곱하기 */
XDesign = j(nrow(X), deg+1);
do j = 0 to deg; /* 2. 자유도에 따른 다항 행렬 설계 */
Xdesign[,j+1] = X##j;
end;
Xw = sqrt(w)#Xdesign; /* 3. 다항행렬에 가중치 곱하기 */
b = solve(Xw`*Xw, Xw`*Yw); /* 4. 회귀 계수 추정*/
return b;
finish;
use RegData; read all var {Y X w}; close; /* 데이터와 가중치 불러오기*/
맨 첫줄에서 PolyRegEst라는 행렬식을 선언해주었는데요, 이를 사용하여 언제든지 가중회귀 계수를 추정할 수 있습니다.
먼저 Yw를 정의합니다. 여기서 #은 행렬 원소간의 곱을 나타냅니다.(선형대수에 나오는 행렬끼리의 곱셈이 아님) 위에선 11x1인 Y행렬과 11x1인 w제곱근 행렬이 각각 원소끼리 곱해진 것으로 볼 수 있습니다.
(여기서 w가 아닌 w의 제곱근을 곱한 이유는 선형대수식에 의해 w제곱근끼리 곱해져 w 가중치 그대로 산출되기 때문입니다. )
그 밑에 XDesign을 통해 자유도를 포함하는 X식을 만들어줍니다. j(a,b,c) 식은 axb이고 원소가 전부 c 인 행렬을 생성합니다. 여기선 행의 갯수가 nrow(X)이고 열의 개수가 deg+1 인 행렬이 생성되겠군요. 그 바로 밑에는 원소를 넣어주는데 ,X##j란 X의 j승을 말합니다.
XDesign 행렬의 이해를 돕기위해 일반 다항회귀 선형식을 가져왔습니다. 다항회귀는 위와 같은 형식으로 제작되기 때문에 자유도 별로 자승을 해주는 것입니다.
<code3>
/* 중요, 자유도의 개수(추정변수의 개수)에 따른 회귀 계수 추정*/
do deg = 0 to 2;
b = PolyRegEst(Y, X, w, deg);
d = char(deg,1); /* '0', '1', or '2' */
labl = "Estimates (deg=" + d + ")";
print b[L=labl rowname=("b0":("b"+d))]; /* 회귀 계수 출력 */
end;
char함수는 행렬의 원소값을 문자형 변수로 출력해주는 기능을 합니다. 여기선 컬럼이름을 b0, b1, b2로 뽑아내기 위해 char 함수를 사용했습니다.
IML 코드 입력을 완료하면 위와같이 추정계수가 산출됩니다.
이렇게 IML을 통한 가중회귀 구현을 실습해보았습니다.
- Visualization
오늘은 특별히, 각가의 회귀 그래프를 평가 / 비교해보는 시간을 가져보겠습니다 .
/* 회귀식 적합 */
start PolyRegScore(x, coef);
p = nrow(coef);
y = j(nrow(x), 1, coef[p]); /* Y행렬 선언 */
do j = p-1 to 1 by -1;
y = y # x + coef[j];
end;
return(y);
finish;
그래프 산출을 위해 간단하게 계수에 따른 Y햇값을 출력하는 행렬식을 만들어보았습니다.
t = T( do(min(x), max(x), (max(x)-min(x))/25) ); /* X 범위 내의 Uniform Grid 생성 */
Yhat = j(nrow(t), 3);
do d = 0 to 2;
b = PolyRegEst(Y, X, w, d); /* PolyRegEst 불러오기*/
Yhat[,d+1] = PolyRegScore(t, b); /* 산출한 계수를 PolyRegScore에 적용*/
end;
여기서 T에 대한 개념이 조금 어려울 수 있습니다. T는 Uniform Grid를 생성하는 기능을 합니다. Uniform grid란 말그대로 일정한 격자 방안인데요, 즉, 위에서는 x의 범위내에서 일정하게 26개(25개의 경계값)로 쪼개주는 기능을 합니다. 산출해보니 격자 폭은 0.264로 나오네요. 큰 의미는 없습니다.
Z = t || Yhat; /* Deg에 따른 3개의 예측값 행렬 생성*/
create RegFit from Z[c={"t" "Pred0" "Pred1" "Pred2"}];
append from Z;
QUIT;
t값에 따른 y추정값을 Z라는 행렬에 전부 담아보았습니다. 이러면 추정 값들의 비교가 쉽겠지요?
data RegAll; /* 비교용 데이터 생성*/
label w="Weight" Pred0="Weighted Mean"
Pred1="Weighted Linear Fit" Pred2="Weighted Quadratic Fit";
merge RegData RegFit;
run;
기초데이터가 담긴 Regdata와, 적합한 회귀식에서 산출된 예측 값 데이터를 합쳤습니다.
각 회귀 그래프를 비교하기 위해서 입니다.
title "Weighted Regression Models"; /* 그래프 그리기 */
proc sgplot data=RegAll;
scatter x=x y=y / filledoutlinedmarkers markerattrs=(size=12 symbol=CircleFilled)
colorresponse=w colormodel=TwoColorRamp;
series x=t y=Pred0 / curvelabel;
series x=t y=Pred1 / curvelabel;
series x=t y=Pred2 / curvelabel;
xaxis grid; yaxis grid;
run;
PROC SGPLOT에서 Series 문은 그래프를 비교할 떄 주로 사용합니다. 여기에선 deg0,1,2 에따라 산출된 Pred0,1,2 를 비교해보겠습니다. 추정계수의 개수에 따른 그래프 적합성을 확인해 볼 수가 있겠네요.
다항 가중회귀, 가중회귀, 가중평균 그래프 값들이 산출되었습니다. 다항 가중회귀가 참 이쁘게 적합되었군요.
여러분들도 이런 스킬들을 통해 다양한 분석을 실행하고 양질의 결과를 내시길 바랍니다.
5. Summary
세 개의 가중회귀 관련 포스팅을 하면서, 제일 고난이도였던 게시글이 되겠네요. 확실히 PROC REG나 PROC SGPLOT을 이용한 회귀 구현이 쉽게 느껴집니다. 하지만 우리는 때로 고급분석을 진행하며, 프로세스 내에서 섬세한 파라미터 조정이나 수학식 변형을 적용하는 경우가 있습니다. 이런식으로 IML을 통해 한 단계씩 구현하는 역량을 키우면 고급 분석 역량을 키우는 데 많은 도움이 될 것입니다.
어쨋건, 이번에 배워본 가중회귀는 실데이터를 분석할 때 유용하게 쓰일 수 있습니다. 가중치가 클수록 추정값에 많은 영향을 끼치는 것 또한 알아 볼 수 있 습니다. 다음시간에는 Kernel Regression 구현을 다뤄보겠습니다.
