상세 컨텐츠

본문 제목

[ SAS ] - 단순 회귀 모형을 분석하는 프로그램 작성 (2)

Language/SAS

by bing_su 2023. 5. 8. 16:56

본문

반응형
SMALL

저번 시간에는 "단순 회귀 모형(Simple Regression Model)"에 관한 통계 개념에 대하여 학습한 후, 가상 데이터와 실제 데이터를 기반으로 한 단순 회귀 모형을 분석하는 프로그램을 작성했다. 이번 시간에는 저번 실습 시간에 이어서 단순 회귀 모형을 분석하는 여러 프로그램을 작성해보려고 한다.
 
프로그램을 작성하기 이전에, 프로그램 분석에 필요한 통계 개념을 간단하게 짚고 넘어가려고 한다. 단순 회귀 모형에서 조금 더 추가된 내용이니 단순 회귀 모형에 대한 자세한 설명이 필요하신 분은 이전 게시물 보고 오시면 될 듯!

[ 오차(error), 잔차(residual), 예측값(predicted value) ]

1. 오차(error)

모집단에서 실제 관찰된 y값과 회귀 직선에 의해 예측된 Y값 차이.

2. 잔차(residual)

표본에서 실제 관찰된 y값($y_i$)과 회귀 직선에 의해 예측된 Y값($Y_i$) 사이의 차이.
오차는 실제로 측정이 불가능하기 때문에 회귀 분석 시에는 오차 대신 잔차를 이용.

3. 예측값(predicted values)

특정 독립 변수에서 회귀 직선에 의해 예측된 Y값. 즉, 예측값은 회귀 직선 위에 존재하는 Y값이라고 할 수 있다.

- 단순 회귀 모형(Simple Regression Model) ]

단순 회귀 모형에서 $X_i, Y_i$ 변수 사이의 선형적 관계를 표시하면 아래와 같다.

$Y_i=\alpha+\beta X_i+u_i$

  • $Y_i$: i번째 측정된 종속 변수(반응 변수) Y의 값
  • $X_i$: i번째 독립 변수(설명 변수) X의 값
  • $u_i$: i번째 측정된 Y의 오차항(error term)
  • $\alpha, \beta$: 회귀 계수($\alpha$는 절편 회귀 계수, $\beta$는 기울기 회귀 계수)

매개 변수 $\alpha, \beta$의 값을 OLS를 이용하여 $a, b$로 추정한 단순 회귀 모형은 아래와 같다.

$Y_i=a+bX_i+e_i$

  • $a$: OLS로 추정한 $\alpha$ 값
  • $b$: OLS로 추정한 $\beta$ 값
  • $e_i$: 잔차(residual). 표본(sample)으로 추정한 회귀식과 실제 관측값의 차이

이때, $a+bX_i$ 값은 predicted value가 되고, $e_i$는 residual이므로 $Y_i$는 observed value(실제 관측한 값)이 된다고 할 수 있다. 

[ 자유도(DF. Degree of Freedom) ]

지난 회귀 분석 실습 프로그램에서 통계적 내용을 설명할 때, 자유도에 대한 개념을 제대로 짚지 않고 넘어갔다. 이번 강의에서 자유도에 대한 개념이 언급되어 있기 때문에, 원론 수준에서 설명을 해 보려고 한다. (아직 학부생에 불과해서 자유도에 대한 개념이 미숙할 수 있습니다. 이건 좀 감안해주세요🥲)
 
자유도는 원론 수준에서 '값이 자유롭게 변할 수 있는 원소의 개수' 정도로 생각하면 된다. 이 개념을 몇몇 통계적 개념에서 쓰였던 자유도를 토대로 이해해보자.

1. 표본 분산(표준편차)를 구할 때의 자유도:  n-1

표본 분산(편차 제곱의 평균)을 구할 때 자유도는 (n-1)이다. (엄밀히 말하면 평균의 정의도 전체 개수로 나누는 것이 아닌, 자유도로 나누는 것이다.) 표본 분산에 대한 자세한 설명이 필요하다면 이전 게시글을 참고하자.

표본 분산 $s^2=\frac{1}{n-1}\sum_{i=1}^n e_i^2$ ($e_i$: 잔차)

이를 이해하기 위해 n개의 원소로 이루어진 표본을 전체 모집단에서 추출했다고 가정하자. 이때, 개별 원소의 값은 모르지만 표본 평균을 알고 있다고 가정하자. 이러한 상황은 아래와 같이 표현될 수 있다.

  • $x_1+x_2+x_3+\cdots+x_n=n \times \bar{X}$
  • $x_k$: 표본 내 $k$번째 원소
  • $\bar{X}$: 표본 평균

표본 평균이 주어진 상황에서 1~(n-1)번째 원소의 값이 주어지면 나머지 $x_n$ 값은 자동으로 결정되기 때문에 자유롭게 변할 수 있는 원소의 개수, 즉, 자유도는 (n-1)이 된다.

2. 잔차의 분산(표준편차)를 구할 때의 자유도: n-2

표본 분산을 구할 때와 동일한 원리로 잔차의 분산에서 자유도가 (n-2)로 계산되는 과정 역시 이해할 수 있다. 먼저 주어진 '표본'을 회귀 분석하여 회귀 직선을 그렸다고 가정하자. 회귀 직선 자체는 아래와 같은 형태로 나온다.

$Y=a+b \times X$

개별 원소로만 따지면 $Y_i=a+bX_i+e_i$의 형태를 띌 것이다. 전체 표본의 수(X의 수)가 n이라고 가정하면, 잔차의 자유도(독립적인 잔차의 수)는 x의 자유도보다 1이 적을 것이다. 왜냐하면, 개별 y값, 개별 x값, 잔차의 평균은 주어져 있기 때문이다. (잔차의 평균은 OLS의 First Normal Equation에 의해 0이 됨을 증명할 수 있다. 자세한 증명은 이 링크를 참고하도록 하자.) 이에 대해 조금 더 설명하기 위해 예시를 하나 살펴보자. 

No.(obs) y x e
1 10 29 0
2 11 30 1
3 12 31 2
4 13 32 3
5 14 33 -3
6 15 34 -2
7 16 35 -1
mean 13 32 0

순차적으로 살펴보자. (개별 y값은 주어진다고 가정.) 첫 행부터 순차적으로 x값(29)과 e값(0)을 구한다. 두번째 행도 무엇이 나올지 모르니 또 결과를 받아봐야 한다. x값(30)과 잔차(1)을 구할 수 있다. 비슷한 방식으로 다섯번째 행까지 결과를 뽑는다.
 
여섯번째, 일곱번째 행의 x값을 뽑는 상황을 생각해보자.

  1. $x_1 \sim x_6$와 표본 평균을 알고 있기 때문에 $x_7$ 값은 뽑지 않아도 자동 결정된다. 따라서 $x$의 분산을 구할 때 자유도는 $7-1=6$으로 구할 수 있다.
  2. $x_6$의 결과를 받는 순간, 회귀 직선을 통해 $e_6$값은 정해진다. $x_6$값을 알면 $x_7$값은 자동 결정되고, 따라서 $e_7$의 값도 정해진다. 따라서 $e_6, e_7$값은 따로 결과를 받아볼 필요 없이 정해진다고 할 수 있다. 따라서 $e$의 분산을 구할 때 자유도는 $7-2=5$로 구할 수 있다. 

※ 2 - 보충 설명

회귀 직선에서 개별 y값과 표본 평균이 주어질 때를 전제로 $x$와 $e$의 분산을 구할 때 자유도를 살펴봤다. 이때, 개별 y값을 알고 있기 때문에 x값을 알면 e값을 자동으로 알 수 있는 것이 아니냐는 의문이 들 수 있다. 그러나 $x_1 \sim x_5$까지는 자유롭게 변할 수 있는 값이기 때문에 각 x값에 따라 e값도 자유롭게 변할 수 있다.
 
$e_6$의 경우에는 $x_6$값을 받음과 동시에 결정되는 값이기 때문에 변할 수 있는 값이라고는 할 수 없다. 또한, $x_7$ 역시 $x_6$값을 받음과 동시에 결정되는 값이므로, $e_7$ 역시 정해진 값이라고 할 수 있다. 따라서 $e$의 분산을 구할 때 자유도는 $7-2=5$ 꼴로 구할 수 있게 된다.
 
위의 개념들을 토대로 이제 이번 실습시간에 배울 프로그램에 대해 분석해 보려고 한다. 

[ 가상 데이터를 단순 회귀 모형으로 분석하는 프로그램 ]

DATA fun;
	seed=1234;
	alpha=10;  /* 절편 설정(true value) */
	beta=1;  /* 기울기 설정(true value) */
	
	DO i=1 TO 200;
		x=i;  /* 독립 변수 값 */
		u=3*RANNOR(seed);  /* 오차항. Var(u)=9 */
		y=alpha+beta*x+u;  /* 종속 변수 값 */
		OUTPUT;  /* DO Loop 내부의 변수 값들을 모두 기억 */
	END;
RUN;
/* regression */
PROC REG DATA=fun;
	MODEL y=x;  /* regression model */
    /* regression 결과 Predicted(P)=pred로, Residual=resid로
    표기하여 이 변수들을 DATA 바구니 out1에 담는다.*/
	OUTPUT OUT=out1 P=pred R=resid;
RUN;

DATA fun;
	/* out1 DATA 바구니를 fun 테이블에 담는다. */
	SET out1;
    /* Observed value(prre)=Predicted value(pred)+Residual(resid) */
	prre=pred+resid;
RUN;
/* plotting */
PROC GPLOT DATA=fun;
	PLOT pred*x prre*x / OVERLAY;  /* pred*x prre*x 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*x에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* prre*x에 사용하는 SYMBOL Setting */
RUN;
QUIT;  /* SAS Code 마지막에 PROC GPLOT을 하면 PLOT을 2번 함. 이를 방지하기 위해 QUIT 사용 */

저번 실습 마지막 프로그램과 유사하게 가상의 단순 회귀 모형을 만들어 회귀 분석을 진행하는 프로그램이다. 결국 두 프로그램은 동일한 Logic을 가진다. 이에 대해 설명하기 위하여 저번 시간에 작성한 프로그램 코드를 아래에 삽입했다.

DATA fun;
	seed=12;
	alpha=1.0;  /* 절편 설정(true value) */
	beta=0.1;  /* 기울기 설정(true value) */
	DO i=1 TO 200;
		x=i;  /* 독립 변수 값 */
		u=2*RANNOR(seed);  /* 오차항. Var(u)=4 */
		y=alpha+beta*x+u;  /* 종속 변수 값 */
		OUTPUT;  /* DO Loop 내부의 변수 값들을 모두 기억 */
	END;
RUN;
/* regression */
PROC REG DATA=fun;
	MODEL y=x;  /* regression model */
    /* regression 결과 Predicted(P)=pred로, Residual(R)=resid로
    표기하여 이 변수들을 DATA 바구니 out1에 담는다.*/
	OUTPUT OUT=out1 P=pred R=resid;
RUN;

DATA fun;
	/* out1 DATA 바구니를 fun DATA 테이블에 담는다. */
	SET out1;
RUN;
/* plotting */
PROC GPLOT DATA=fun;
	PLOT y*x;  /* 세로축을 y, 가로축을 x로 지정*/
	SYMBOL V=STAR I=JOIN C=BLUE;  /* GRAPH Option Setting */
RUN;
/* plotting */
PROC GPLOT DATA=fun;
	PLOT pred*x y*x/OVERLAY;  /* pred*x y*x 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*x에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* y*x에 사용하는 SYMBOL Setting */
RUN;
QUIT;  /* SAS Code 마지막에 PROC GPLOT을 하면 PLOT을 2번 함. 이를 방지하기 위해 QUIT 사용 */

두 프로그램에서 달라 보이는 부분은 GPLOT Procedure다. 이 부분 코드만 따로 살펴보자.

PROC GPLOT DATA=fun;
	/* prre=pred+resid */
	PLOT pred*x prre*x / OVERLAY;  /* pred*x prre*x 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*x에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* prre*x에 사용하는 SYMBOL Setting */
RUN;

이번 실습 때 작성한 GPLOT Procedure다. 위의 Procedure는 pred*x와 prre*x 변수의 관계를 나타내는 그래프를 겹쳐 그린다. 여기서 prre 변수는 Predicted value와 Residual 값을 합한 Observed value가 된다. 즉, prre 변수는 우리가 프로그램을 통해 추정한 단순 회귀 모형의 $y$ 값이다.

PROC GPLOT DATA=fun;
	PLOT pred*x y*x/OVERLAY;  /* pred*x y*x 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*x에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* y*x에 사용하는 SYMBOL Setting */
RUN;

저번 실습 때 작성한 GPLOT Procedure다. 위의 Procedure는 pred*x와 y*x 변수의 관계를 나타내는 그래프를 겹쳐 그린다. 즉, 위의 prre*x 그래프와 지금의 y*x 그래프는 동일한 그래프라고 할 수 있다. 결국 두 프로그램은 동일한 그래프를 그리고 있고, 출력 결과 역시 동일한 양상으로 나타날 것이다.

DATA fun;
	seed=1234;
	alpha=10;  /* 절편 설정(true value) */
	beta=1;  /* 기울기 설정(true value) */
	
	DO i=1 TO 200;
		x=i;  /* 독립 변수 값 */
		u=3*RANNOR(seed);  /* 오차항. Var(u)=9 */
		y=alpha+beta*x+u;  /* 종속 변수 값 */
		OUTPUT;  /* DO Loop 내부의 변수 값들을 모두 기억 */
	END;
RUN;
/* regression */
PROC REG DATA=fun;
	MODEL y=x;  /* regression model */
    /* regression 결과 Predicted(P)=pred로, Residual=resid로
    표기하여 이 변수들을 DATA 바구니 out1에 담는다.*/
	OUTPUT OUT=out1 P=pred R=resid;
RUN;

DATA fun;
	/* out1 DATA 바구니를 fun 테이블에 담는다. */
	SET out1;
    /* Observed value(prre)=Predicted value(pred)+Residual(resid) */
	prre=pred+resid;
RUN;
/* plotting */
PROC GPLOT DATA=fun;
	PLOT pred*x prre*x / OVERLAY;  /* pred*x prre*x 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*x에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* prre*x에 사용하는 SYMBOL Setting */
RUN;
QUIT;  /* SAS Code 마지막에 PROC GPLOT을 하면 PLOT을 2번 함. 이를 방지하기 위해 QUIT 사용 */

다시 이번 실습 때 작성한 프로그램으로 돌아와서, 해당 프로그램을 실행한 결과에 대해 분석해 보자.

프로그램 실행 화면 - 1

200개의 sample로 회귀 분석을 한 결과이다. R-Square 값은 0.9970으로, 이는 독립 변수가 종속 변수의 99.70% 정도를 설명한다고 할 수 있다. 우리가 설정한 True Intercept(실제 절편 값)은 $\alpha=10$이었고, 프로그램을 통해 계산된 OLS 추정치는 9.75118이다. 또한, 실제 기울기 값은 $\beta=1$이었고, 프로그램을 통해 계산된 OLS 추정치는 1.00542 임을 확인할 수 있다.
 
$\alpha$의 OLS 추정치 $a$의 t 값은 21.60으로 95% critical value인 1.98보다 훨씬 크다. 따라서 귀무가설 $\alpha_0=0$는 기각된다. 또한, $\beta$의 OLS 추정치 $b$의 t 값은 258.17, 즉, 95% critical value인 1.98보다 크므로 귀무가설 $\beta_0=0$ 역시 기각된다. $a, b$의 t Value를 계산하는 과정은 아래에 자세히 있으니 살펴보시면 됨.

  • $a$의 t Value$=\frac{a-\alpha_0}{\sqrt{s^2(\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n x_i^2})}}=\frac{Parameter Estimate}{Standard Error}=\frac{9.75118}{0.45138}=21.60>1.98$
  • $b$의 t Value$=\frac{b-\beta_0}{\sqrt{\frac{s^2}{\sum_{i=1}^n x_i^2}}}=\frac{Parameter Estimate}{Standard Error}=\frac{1.00542}{0.00389}=258.17>1.98$

t 값이 클수록 모수(parameter) 추정치가 유의하다고 한다. 즉, 귀무가설 $\alpha_0=0, \beta_0=0$가 강하게 기각된다, 해당 프로그램에서 분석된 t 값을 이용하여 위에서 $\alpha_0=0, \beta_0=0$ 귀무가설을 검증했다. 결과는 두 가설 모두 기각되어야 한다.
 
P 값으로도 해당 귀무가설을 검증할 수 있다. P > | t |는 귀무가설 $\alpha=0, \beta=0$일 때 계산된 t 값이 21.60(258.17) 보다 커질 확률이 0.0001일 확률을 나타낸다. 즉, 95%의 신뢰 수준에서 P <0.05이므로 해당 귀무가설은 기각되어야 한다. 즉, t 값과 P 값은 우리에게 같은 결론을 말해주지만, 동일한 결론에 도달하기 위한 기준이 다르다. t 값은 그 자체로 우리에게 귀무가설을 검증하게 해 주는 반면, P 값은 확률 값으로 환산되어 우리가 귀무가설을 검증하게 해 준다.

프로그램 실행 화면 - 2

위의 그래프는 잔차도(Residual plot)이다. 이는 Predicted value(예측값)과 residual(잔차) 사이의 데이터 패턴을 보여준다. 잔차 분석에 대해서는 추후에 다룰 예정이니 이 그래프 분석은 이쯤에서 마무리하려고 한다.

프로그램 실행 화면 - 3

위의 그림은 fit plot이다. 이는 예측된 회귀 직선과 산점도(scatter plot)을 합쳐 놓은 그림이다. 확대해서 봐야 보이지만, 파란색으로 색칠된 면적은 95% confidence limits이며 점선은 95% prediction limits이다.

프로그램 실행 화면 - 4

GPLOT Procedure를 통해 pred*x 관계를 도식화한 그림과 prre*x(y*x) 관계를 도식화한 그림을 겹쳐서 나타냈다.

[ 실제 데이터를 이용하여 단순 회귀 분석을 하는 프로그램 ]

/* 파일에서 실제 DATA를 추출하는 Logic은 이전에 설명한 부분이라 생략함. */
/* 기억 잘 안 나시는 분들은 https://bing-su-b.tistory.com/109 ㄱ-ㄱ */
DATA ip;
	INFILE '/home/u63329964/mydata/ip.prn';
	INPUT mon ip;
	logip=LOG(ip);
	ipg=DIF(logip)*1200;
    /* 1959년 이전 삭제, 코로나 초기 포함 */
    /* 코로나 기간을 제외하기 위해서는 mon>20200101로 코드 수정 필요 */
	IF mon<19590101 OR mon>20210201 THEN DELETE;
RUN;

DATA fyff;
	INFILE '/home/u63329964/mydata/fyff.prn';
	INPUT mon fyff;
	fyff4=LAG4(fyff);  /* 4달 전 값을 fyff4 변수로 정의 */
    /* 1959년 이전 삭제, 코로나 초기 포함 */
    /* 코로나 기간을 제외하기 위해서는 mon>20200101로 코드 수정 필요 */
	IF mon<19590101 OR mon>20210201 THEN DELETE;
RUN;

DATA all;
	MERGE ip fyff;
	BY mon;  /* mon 변수를 기준으로 병합 */
	num=_N_;
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT logip*num;  /* 세로축을 logip, 가로축을 num으로 지정 */
	SYMBOL VALUE=DOT INTERPOL=JOIN COLOR=BLUE;  /* GRAPH Option Setting */
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT ipg*num;  /* 세로축을 ipg, 가로축을 num으로 지정 */
	SYMBOL VALUE=DOT INTERPOL=JOIN COLOR=RED;  /* GRAPH Option Setting */
RUN;
/* regression */
PROC REG DATA=all;
	MODEL ipg=fyff;  /* regression model 1 */
	MODEL ipg=fyff4;  /* regression model 2 */
RUN;

해당 프로그램에서 사용할 데이터는 fyff(the U.S Federal funds interest rate, 미국 연방기금 금리에 대한 자료)와 ipg(monthly industrial production, 산업생산의 월별 성장률), 그리고 두 데이터를 날짜를 기준으로 merge한 데이터이다. 이 프로그램은 해당 데이터를 가지고 회귀 분석을 진행하는 프로그램이다. 분석에 대한 내용은 프로그램 실행 결과를 살펴보며 설명할 예정이다.

1. PLOT logip*num

DATA all;
	MERGE ip fyff;
	BY mon;  /* mon 변수를 기준으로 병합 */
	num=_N_;
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT logip*num;  /* 세로축을 logip, 가로축을 num으로 지정 */
	SYMBOL VALUE=DOT INTERPOL=JOIN COLOR=BLUE;  /* GRAPH Option Setting */
RUN;

PROC GPLOT 결과 - 1

num 변수는 mon 변수의 일련 번호라고도 할 수 있다. logip와 num 사이의 상관관계는 logip와 mon 사이의 상관관계라고도 할 수 있다. 즉, 해당 GPLOT Procedure는 logip 변수와 mon 변수의 상관관계를 그래프로 나타냈다고 할 수 있다. 요약하자면 해당 Procedure는 월간 산업생산의 지표(index of industrial production)를 표현한 것이라고 할 수 있다.
 
위의 그래프를 통해 분석할 수 있는 내용은 아래와 같다.

  • 1차 오일 쇼크 불경기 - 19750501: num=677 저점
  • 2차 오일 쇼크 불경기 - 19821201: num=768 저점
  • 글로벌 금융위기 불경기 - 20090601: num=1086 저점
  • 코로나 펜데믹 불경기 - 20200401: num=1216 저점

2. PLOT ipg*num

DATA all;
	MERGE ip fyff;
	BY mon;  /* mon 변수를 기준으로 병합 */
	num=_N_;
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT ipg*num;  /* 세로축을 ipg, 가로축을 num으로 지정 */
	SYMBOL VALUE=DOT INTERPOL=JOIN COLOR=RED;  /* GRAPH Option Setting */
RUN;

PROC GPLOT 결과 - 2

num 변수는 위에서 설명했듯이 mon 변수의 일련 번호라고 할 수 있다. 즉, ipg와 num 사이의 상관관계는 ipg와 mon 사이의 상관관계라고도 할 수 있다. 즉, 해당 GPLOT Procedure는 ipg 변수와 mon 변수의 상관관계를 그래프로 나타낸 것이다. 요약하자면 해당 Procedure는 월간 산업생산의 성장률을 표현한 것이라고 할 수 있다.
 
위의 그래프를 통해 분석할 수 있는 내용은 아래와 같다.

  • 1차 오일 쇼크 불경기 - 19750501: num=677 저점
  • 2차 오일 쇼크 불경기 - 19821201: num=768 저점
  • 글로벌 금융위기 불경기 - 20090601: num=1086 저점
  • 코로나 펜데믹 불경기 - 20200401: num=1216 저점

3. PROC REG 결과 1: MODEL ipg=fyff

fyff(미국 연방기금 금리)를 독립 변수로 사용하여 ipg(산업생산의 연간 생산률, 종속 변수)를 예측하기 위한 regression model을 설정했다. 이 모형에 대한 회귀 계수를 추정하는 Procedure의 실행 결과는 아래와 같다.

MODEL ipg=fyff의 Regression Procedure 결과

총 샘플 수는 19590101~20210201까지 총 746이다. 따라서 잔차의 sum of squares를 구할 때 자유도는 $746-2=744$임을 확인할 수 있고, 표본의 sum of squares를 구할 때 자유도는 $746-1=745$임을 확인할 수 있다.
 
해당 결과 표에서  'MODEL DF=1'임을 알 수 있다. 자유도 개념을 토대로 해당 결과를 해석해보자. 회귀 직선은 $Y_i=a+bX_i+e_i$로 표현된다. 이때, $Y_i$, $X_i$, $e_i$의 개별 원소 값을 토대로 $b$ 값을 추정할 수 있다. (추정 방법에 대해서는 이 링크의 OLS Method 관련 설명 보시면 자세히 써 있습니다.) b값을 구하면 a값은 $a=\bar{Y}-b\bar{X}$에 의해 자동 결정된다. 따라서 해당 선형 회귀 모델에서 자유롭게 변할 수 있는 원소는 개별 원소의 측정 값을 통해 결정되는 b값, 하나라고 할 수 있다. 따라서 'MODEL DF=1'이라는 결과를 얻을 수 있다.
 
이 외의 R-Square, t value, P value 등에 대한 해석은 위에서 다루었기 때문에 생략한다.

4. PROC REG 결과 1: MODEL ipg=fyff4

fyff(미국 연방기금 금리)에 따라 4달 후의 ipg(산업생산의 연간 생산률, 종속 변수)를 예측하기 위한 regression model을 설정했다. 이 모형에 대한 회귀 계수를 추정하는 Procedure의 실행 결과는 아래와 같다.

MODEL ipg=fyff4의 Regression Procedure 결과

총 샘플 수는 첫번째 모델과 동일하게 19590101~20210201까지 총 746이며, 잔차의 sum of squares를 구할 때 자유도는 $746-2=744$, 표본의 sum of squares를 구할 때 자유도는 $746-1=745$임을 확인할 수 있다.
 
회귀 직선 상의 a(intercept)에 대한 통계값은 아래와 같이 확인할 수 있다.

  • $a=Intercept=3.94413$
  • (Estimated) Standard Error of $a=0.74194$
  • t Value$=5.32=\frac{a}{s \sqrt{\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n x_i^2}}}$

회귀 직선 상의 b(fyff4. slope)에 대한 통계값은 아래와 같이 확인할 수 있다.

  • $b=fyff=-0.31869$
  • (Estimated) Standard Error of $b=0.12159$
  • t Value$=-2.62=\frac{b}{\frac{s}{\sqrt{\sum_{i=1}^n x_i^2}}}$

 
MODEL DF=1, R-Square, P value에 대한 해석은 위에서 충분히 다루었기 때문에 생략한다.

[ fyff를 내렸을 때 ipg가 가장 큰 영향을 받는 시기를 추정하는 프로그램 ]

/* 파일에서 실제 DATA를 추출하는 Logic은 이전에 설명한 부분이라 생략함. */
/* 기억 잘 안 나시는 분들은 https://bing-su-b.tistory.com/109 ㄱ-ㄱ */
DATA ip;
	INFILE '/home/u63329964/mydata/ip.prn';
	INPUT mon ip;
	logip=LOG(ip);
	ipg=DIF(logip)*1200;
    /* 1959년 이전 삭제. 코로나 기간 제외. */
	IF mon<19590101 OR mon>20200101 THEN DELETE;
RUN;

DATA fyff;
	INFILE '/home/u63329964/mydata/fyff.prn';
	INPUT mon fyff;
	fyff1=LAG(fyff);  /* 1달 전 값을 fyff1 변수로 정의 */
	fyff2=LAG2(fyff);  /* 2달 전 값을 fyff2 변수로 정의 */
	fyff3=LAG3(fyff);  /* 3달 전 값을 fyff3 변수로 정의 */
	fyff4=LAG4(fyff);  /* 4달 전 값을 fyff4 변수로 정의 */
	fyff5=LAG5(fyff);  /* 5달 전 값을 fyff5 변수로 정의 */
	fyff6=LAG6(fyff);  /* 6달 전 값을 fyff6 변수로 정의 */
    /* 1959년 이전 삭제. 코로나 기간 제외. */
	IF mon<19590101 OR mon>20200101 THEN DELETE;
RUN;

DATA all;
	MERGE ip fyff;
	BY mon;  /* mon 변수를 기준으로 병합 */
	n1=19590101;  /* 시작 날짜 */
	n2=20200101;  /* 종료 날짜 */
    /* 1959년 이전 삭제. 코로나 기간 제외. */
	IF mon<n1 OR mon>n2 THEN DELETE;
	num=_N_;
RUN;
/* regression */
PROC REG DATA=all;
	MODEL ipg=fyff;  /* regression model 0 */
	MODEL ipg=fyff1;  /* regression model 1 */
	MODEL ipg=fyff2;  /* regression model 2 */
	MODEL ipg=fyff3;  /* regression model 3 */
	MODEL ipg=fyff4;  /* regression model 4 */
	MODEL ipg=fyff5;  /* regression model 5 */
	MODEL ipg=fyff6;  /* regression model 6 */
RUN;

fyff(미국 연방기금 금리)를 내렸을 때 ipg(산업생산의 연간 성장률)가 몇 달 후에 가장 큰 영향을 받을지를 판단하는 프로그램이다. (단, 1달부터 6달 후까지만 고려) 이에 답하기 위해서는 최적의 LAG length를 찾아야한다. 따라서 해당 프로그램에서는 LAG length를 1씩 늘리며 여러 가지 회귀 모형을 관찰했다. 그 중에서 기울기의 t 값이 가장 큰(귀무가설 \beta_0=0과 가장 유의한 차이가 있는) LAG를 선택하면 될 것이다. 해당 프로그램을 실행한 결과 중 일부 DATA를 추출하여 표로 나타내면 아래와 같다.

Regressor Estimate t Value
fyff -0.14149 -1.45
fyff1 -0.23371 -2.40
fyff2 -0.29450 -3.04
fyff3 -0.34248 -3.54
fyff4 -0.36136 -3.74
fyff5 -0.35359 -3.66
fyff6 -0.33860 -4.37

fyff4에서 가장 큰 t 값을 가지므로, fyff(미국 연방기금 금리)를 내렸을 때, ipg(산업생산의 연간 성장률)는 4달 후에 가장 큰 영향을 받음을 확인할 수 있다.

[ 가장 좋은 regressor인 fyff4를 사용한 단순 회귀 모형을 분석하는 프로그램 ]

/* 파일에서 실제 DATA를 추출하는 Logic은 이전에 설명한 부분이라 생략함. */
/* 기억 잘 안 나시는 분들은 https://bing-su-b.tistory.com/109 ㄱ-ㄱ */
DATA ip;
	INFILE '/home/u63329964/mydata/ip.prn';
	INPUT mon ip;
	logip=LOG(ip);
	ipg=DIF(logip)*1200;
    /* 1959년 이전 삭제. 코로나 포함 X. */
	IF mon<19590101 OR mon>20200101 THEN DELETE;
RUN;

DATA fyff;
	INFILE '/home/u63329964/mydata/fyff.prn';
	INPUT mon fyff;
	fyff4=LAG4(fyff);
    /* 1959년 이전 삭제. 코로나 포함 X. */
	IF mon<19590101 OR mon>20200101 THEN DELETE;
RUN;

DATA all;
	MERGE ip fyff;
	BY mon;  /* mon 변수를 기준으로 병합 */
	num=_N_;
RUN;
/* regression */
PROC REG DATA=all;
	MODEL ipg=fyff4;  /* regression model */
    /* regression 결과 Predicted(P)=pred로, Residual(R)=resid로
    표기하여 이 변수들을 DATA 바구나 out1에 담는다. */
	OUTPUT OUT=out1 P=pred R=resid;
RUN;

DATA all;
	/* out1 DATA 바구니를 all 테이블에 담는다. */
	SET out1;
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT pred*fyff4;
	SYMBOL V=DOT I=JOIN C=BLUE;
RUN;
/* plotting */
PROC GPLOT DATA=all;
	PLOT resid*num;  /* 세로축을 resid, 가로축을 num으로 지정 */
	SYMBOL V=STAR I=JOIN C=VIOLET;  /* GRAPH Option Setting */
RUN;

PROC GPLOT DATA=all;
	PLOT pred*fyff4 ipg*fyff4/OVERLAY;  /* pred*fyff4 ipg*fyff4 그림을 겹쳐 그림(OVERLAY Option) */
	SYMBOL1 V=DOT I=JOIN C=BLUE;  /* pred*fyff4에 사용하는 SYMBOL Setting */
	SYMBOL2 V=STAR C=RED;  /* ipg*fyff4에 사용하는 SYMBOL Setting */
RUN;

바로 위 프로그램에서  fyff(미국 연방기금 금리)를 내렸을 때 ipg(산업생산의 연간 성장률)은 4달 후에 가장 큰 영향을 받음을 확인할 수 있었다. 해당 프로그램은 가장 좋은 regressor인 fyff4를 사용한 단순 회귀 모형을 분석하는 프로그램이다.
 
MODEL ipg=fyff4에 대한 분석은 이미 위에서 했기 때문에 생략한다. (동일한 방식으로 해석하면 되니, 잘 모르시는 분들은 위에 스크롤 올리고 관련 내용 읽고 오시면 됩니다.) 이번에는 해당 프로그램을 실행했을 때 plotting 결과를 살펴보려고 한다. 실행 결과는 아래 사진과 같다.

plotting - 1: pred*fyff4

fyff(미국 연방기금 금리)에 따른 4달 후의 ipg의 예측 값을 회귀 직선 위에 나타낸 산점도이다.

plotting - 2: resid*num

num 변수는 mon 변수의 일련 번호라고도 할 수 있다. 즉, 해당 GPLOT Procedure 결과는 시기에 따라 계산된 잔차를 나타낸 산점도라고 할 수 있다.

plotting - 3: pred*fyff4 &amp; ipg*fyff4 OVERLAY

위에서 관찰한 pred*fyff 산점도와 ipg*fyff4의 산점도를 겹쳐서 나타냈다. 즉, 예측 ipg 값을 나타낸 산점도에 실제 관측된 ipg 값을 나타낸 산점도를 겹쳐서 그린 것이다.

[ 정리 - OLS 추정치(a, b)의 통계량 ]

단순 회귀 분석을 마무리하며 OLS 추정치 a, b와 관련된 통계량을 최종적으로 정리해보려고 한다.
자세한 수식이 기억나지 않는 분들은 이전 게시글 보고 오시면 됩니다.

- OLS 추정치 a의 통계량

  • OLS estimator of $\alpha=a=\bar{Y}-b \bar{X}=(\alpha+\beta \bar{X}+\bar{u})-b \bar{X}$
  • $E(a)=E[(\alpha+\beta \bar{X}+\bar{u})-b \bar{X}]=[\alpha+\beta \bar{X}+E(\bar{u})-E(b)\bar{X}]=[\alpha+\beta \bar{X}+0-\beta \bar{X}]=\alpha$
    ($E(\bar{u})=0$, $E(b)=\beta$되는 자세한 과정은 이 게시물 참고)
  • (True) Variance of a$=(\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^nx_i^2})\sigma^2$
  • (True) Standard error of a$=\sqrt{a 분산}=\sqrt{(\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n x_i^2})\sigma^2}=\sigma \sqrt{\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n x_i^2}}$
  • (Estimated) Standard Error of a$=s \sqrt{\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n x_i^2}}$
  • t value of a$=\frac{a}{(Estimated) Standard Error}=\frac{a}{s \sqrt{\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n x_i^2}}}$

- OLS 추정치 b의 통계량

  • OLS estimator of $\beta=b=\frac{\sum_{i=1}^n x_iy_i}{\sum_{i=1}^n x_i^2}=\beta+\frac{\sum_{i=1}^n x_i u_i}{\sum_{i=1}^n x_i^2}$
  • E(b)$=\beta+\frac{\sum_{i=1}^n x_i E(u_i)}{\sum_{i=1}^n x_i^2}=\beta$
  • (True) Variance of b$=\frac{\sigma^2}{\sum_{i=1}^n x_i^2}$
  • (True) Standard error of b$=\sqrt{b 분산}=\sqrt{\frac{\sigma^2}{\sum_{i=1}^n x_i^2}}=\frac{\sigma}{\sqrt{\sum_{i=1}^n x_i^2}}$
  • (Estimated) Standard Error of b$=\frac{s}{\sqrt{\sum_{i=1}^n x_i^2}}$
  • t value of b$=\frac{b}{\frac{s}{\sqrt{\sum_{i=1}^n x_i^2}}}$

 

반응형
LIST

관련글 더보기

댓글 영역