본문 바로가기
카테고리 없음

[STATA] 선형회귀분석 기초

by e-money2580 2023. 1. 10.
반응형


/*
y = a + bx + e

a와 b를 모수(parameter)라고 하고 e를 오차항이라고 한다.
모수 a와 b를 ㄱ가장 잘 추정하는 것이 회귀분석의 핵심 관심사이고 최소제곱법 OLS(ordinary least squares)는 가장 널릴 쓰이는 추정방법으로 잔차 의 제곱합을 최소화하는 a와 b를 선택하는 것
아래 조건이 만족되면 최소제곱법을 통해 얻어진 추정량이 최우수선형불편추정량(BLUE : best linear unbiased estimator)이다

(가정1) 실제 모형이 추정모형과 마찬가지로 선형이다
(가정2) 오차항의 조건부 기대값 E(e|x)이 0이다. 즉, cov(x,e)=0이다. 독립변수와 오차항의 공분산은 0이다.
(가정3) 오차항의 조건부 분산 var(ei|xi)은 모든 i에 대해 모두 동일하다(동분산성). 즉, var(ei|xi)=시그마제곱
(가정4) 두 개의 서로 다른 오차항 ei와 ej는 서로 독립이다. 즉, cov(ei,ej)=0이다. 
(가정5) 오차항은 정규분포를 따른다. ei|xi ~ N(0,시그마2)
(가정6) (상수항을 포함하여) 설명변수들 간에 완전한 선형관계가 없어야 한다.
*/


** 단순선형회귀모형 추정(독립변수 1개)
use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_1.dta"

reg food_exp income

/*
      Source |       SS           df       MS      Number of obs   =        40
-------------+----------------------------------   F(1, 38)        =     23.79
       Model |   190626.98         1   190626.98   Prob > F        =    0.0000
    Residual |  304505.173        38  8013.29403   R-squared       =    0.3850
-------------+----------------------------------   Adj R-squared   =    0.3688
       Total |  495132.153        39  12695.6962   Root MSE        =    89.517

------------------------------------------------------------------------------
    food_exp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      income |   10.20964   2.093263     4.88   0.000     5.972052    14.44723
       _cons |   83.41601   43.41016     1.92   0.062    -4.463272    171.2953
------------------------------------------------------------------------------
*/

corr food_exp income

/*
             | food_exp   income
-------------+------------------
    food_exp |   1.0000
      income |   0.6205   1.0000
*/

reg food_exp income, nocons

/*
      Source |       SS           df       MS      Number of obs   =        40
-------------+----------------------------------   F(1, 39)        =    394.28
       Model |  3377595.38         1  3377595.38   Prob > F        =    0.0000
    Residual |  334093.953        39  8566.51161   R-squared       =    0.9100
-------------+----------------------------------   Adj R-squared   =    0.9077
       Total |  3711689.33        40  92792.2333   Root MSE        =    92.555

------------------------------------------------------------------------------
    food_exp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      income |    14.0122   .7056746    19.86   0.000     12.58484    15.43956
------------------------------------------------------------------------------
*/

// 상수항이 없는 회귀결과 R2가 0.91이고 상수항이 있는 회귀의 R2인 0.37보다 높다고 설명력이 높다고 할 수 없다. R2가 서로 다른 방식으로 계산되기 때문에 서로 비교할 수 없다.



** 다중선형회귀모형 추정(독립변수 2개 이상)

use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_2.dta",clear

reg csat percent income

/*

      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =     84.94
       Model |  166809.926         2  83404.9632   Prob > F        =    0.0000
    Residual |  46151.4535        47   981.94582   R-squared       =    0.7833
-------------+----------------------------------   Adj R-squared   =    0.7741
       Total |   212961.38        49  4346.15061   Root MSE        =    31.336

------------------------------------------------------------------------------
        csat | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     percent |  -2.582867   .2347331   -11.00   0.000    -3.055089   -2.110645
      income |   .0021968   .0009402     2.34   0.024     .0003054    .0040882
       _cons |   962.2196   27.35218    35.18   0.000     907.1942    1017.245
------------------------------------------------------------------------------
*/
// percent 회귀추정치인 -2.58 : SAT 시험을 치르는 고등학교 졸업자 비율이 1%p 증가할 때 2.582점 떨어짐 
// 이때 다른 변수(income)는 변화가 없는 상태 가정
// 다중선형회귀모형에서 계수 추정치는 다른 변수가 고정되어 있을 때 해당 변수의 기울기를 의미
// income 변수는 금액이 수만달러이므로 1달러 변화 변화가 미미한 변화이므로, percent 변수의 계수가 훨씬 크다고 해서 percent 변수의 변화가 종속변수에 더 큰 영향을 반드시 미친다고 해석할 수 없음
// -> 측정단위를 표준화할 필요

reg csat percent income, beta

/*
      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =     84.94
       Model |  166809.926         2  83404.9632   Prob > F        =    0.0000
    Residual |  46151.4535        47   981.94582   R-squared       =    0.7833
-------------+----------------------------------   Adj R-squared   =    0.7741
       Total |   212961.38        49  4346.15061   Root MSE        =    31.336

------------------------------------------------------------------------------
        csat | Coefficient  Std. err.      t    P>|t|                     Beta
-------------+----------------------------------------------------------------
     percent |  -2.582867   .2347331   -11.00   0.000                -1.017304
      income |   .0021968   .0009402     2.34   0.024                 .2160285
       _cons |   962.2196   27.35218    35.18   0.000                        .
------------------------------------------------------------------------------
*/
// Beta : 모든 변수를 표준화 시킨 다음 선형회귀모형으로 추정했을 때의 계수


pwcorr csat percent income // 결측치가 존재할 경우 해당 개체(또는 시점)의 관측치 모두 삭제(짝별 삭제)
/*
             |     csat  percent   income
-------------+---------------------------
        csat |   1.0000 
     percent |  -0.8758   1.0000 
      income |  -0.4744   0.6786   1.0000 
*/

corr csat percent income // 결측치가 존재할 경우 사례별 삭제
/*

             |     csat  percent   income
-------------+---------------------------
        csat |   1.0000
     percent |  -0.8707   1.0000
      income |  -0.4744   0.6786   1.0000
*/


// 종속변수에 자연로그를 취한 후 회귀
gen ln_csat=ln(csat)
reg ln_csat percent income
/*
      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =     88.08
       Model |  .186137847         2  .093068924   Prob > F        =    0.0000
    Residual |  .049662574        47  .001056651   R-squared       =    0.7894
-------------+----------------------------------   Adj R-squared   =    0.7804
       Total |  .235800421        49  .004812253   Root MSE        =    .03251

------------------------------------------------------------------------------
     ln_csat | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     percent |  -.0027421   .0002435   -11.26   0.000    -.0032319   -.0022522
      income |   2.42e-06   9.75e-07     2.48   0.017     4.60e-07    4.38e-06
       _cons |   6.864061   .0283736   241.92   0.000     6.806981    6.921141
------------------------------------------------------------------------------
*/

// 종속변수와 독립변수에 자연로그 취한 후 회귀
gen ln_income=ln(income)
reg ln_csat percent income
/*
      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =     88.08
       Model |  .186137847         2  .093068924   Prob > F        =    0.0000
    Residual |  .049662574        47  .001056651   R-squared       =    0.7894
-------------+----------------------------------   Adj R-squared   =    0.7804
       Total |  .235800421        49  .004812253   Root MSE        =    .03251

------------------------------------------------------------------------------
     ln_csat | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     percent |  -.0027421   .0002435   -11.26   0.000    -.0032319   -.0022522
      income |   2.42e-06   9.75e-07     2.48   0.017     4.60e-07    4.38e-06
       _cons |   6.864061   .0283736   241.92   0.000     6.806981    6.921141
------------------------------------------------------------------------------
*/

// 독립변수에 자연로그 취한 후 회귀
reg csat percent ln_income
/*
      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =     84.00
       Model |  166407.418         2  83203.7092   Prob > F        =    0.0000
    Residual |  46553.9616        47  990.509821   R-squared       =    0.7814
-------------+----------------------------------   Adj R-squared   =    0.7721
       Total |   212961.38        49  4346.15061   Root MSE        =    31.472

------------------------------------------------------------------------------
        csat | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     percent |  -2.574919   .2376735   -10.83   0.000    -3.053056   -2.096781
   ln_income |   74.03974     33.091     2.24   0.030     7.469234    140.6102
       _cons |   265.3713   339.0011     0.78   0.438    -416.6109    947.3534
------------------------------------------------------------------------------
*/



** 회귀모형으로 추정된 종속변수의 적합값(fitted value) 및 잔차 구하기

reg csat percent income

predict yhat // fitted value 구하기 - yhat 변수에 fitted value 생성

predict ehat, resid // 잔차 구하기 = ehat 변수에 잔차 생성

list csat yhat ehat in 1/10

/*
     | csat       yhat        ehat |
     |-----------------------------|
  1. |  991   1001.965   -10.96527 |
  2. |  920   962.3282   -42.32817 |
  3. |  932   965.5681   -33.56812 |
  4. | 1005   1000.859     4.14096 |
  5. |  897   932.4681   -35.46806 |
     |-----------------------------|
  6. |  959    964.476   -5.475932 |
  7. |  897   859.8132    37.18686 |
  8. |  892   893.9464   -1.946328 |
  9. |  840          .           . |
 10. |  882      908.6   -26.60006 |
*/
// 종속변수(csat) - 적합값(fitted value, yhat) = 잔차(ehat)

reg csat percent income
gen e_sample=e(sample) // 관측치 중 결측치(0)를 e_sample 변수에 더미변수 형태로 기록 
list csat percent income e_sample in 5/12
/*
     | csat   percent   income   e_sample |
     |------------------------------------|
  5. |  897        47    41716          1 |
  6. |  959        29    35123          1 |
  7. |  897        81    48618          1 |
  8. |  892        61    40641          1 |
  9. |  840        71        .          0 |
     |------------------------------------|
 10. |  882        48    32027          1 |
 11. |  844        62    33819          1 |
 12. |  883        55    45248          1 |
     +------------------------------------+
*/


** 추정계수의 유의성 검정
// OLS 추정량이 좋은 추정량(BLUE)이기 위한 가정 하에서 모수 b(추정계수)의 유의성에 대한 검정

use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_3.dta", clear

reg lprice lnox ldist lproptax crime
/*
      Source |       SS           df       MS      Number of obs   =       506
-------------+----------------------------------   F(4, 501)       =     93.46
       Model |   36.143062         4   9.0357655   Prob > F        =    0.0000
    Residual |  48.4392089       501  .096685048   R-squared       =    0.4273
-------------+----------------------------------   Adj R-squared   =    0.4227
       Total |  84.5822709       505  .167489645   Root MSE        =    .31094

------------------------------------------------------------------------------
      lprice | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        lnox |   -.938093   .1459205    -6.43   0.000    -1.224784   -.6514015
       ldist |  -.2203829   .0514951    -4.28   0.000    -.3215559   -.1192099
    lproptax |  -.2456268   .0509896    -4.82   0.000    -.3458066    -.145447
       crime |  -.0158604   .0019787    -8.02   0.000     -.019748   -.0119728
       _cons |   13.30539    .349259    38.10   0.000      12.6192    13.99159
------------------------------------------------------------------------------
*/
test lnox

/*
 ( 1)  lnox = 0

       F(  1,   501) =   41.33
            Prob > F =    0.0000
*/
// lnox 변수의 추정계수 b=0이라는 귀무가설을 5% 유의수준에서 기각 = 추정계수는 유의함

test lnox=1

/*
 ( 1)  lnox = 1

       F(  1,   501) =  176.41
            Prob > F =    0.0000
*/
// lnox 변수의 추정계수 b=1이라는 귀무가설을 5% 유의수준에서 기각 = 추정계수는 1이 아님

test lnox ldist lproptax crime
/*
 ( 1)  lnox = 0
 ( 2)  ldist = 0
 ( 3)  lproptax = 0
 ( 4)  crime = 0

       F(  4,   501) =   93.46
            Prob > F =    0.0000
*/
// 회귀모형에서 추정한 모든 변수의 추정계수가 0인지 추정한 결과 귀무가설을 5% 유의수준에서 기각
// = 모든 변수의 추정계수는 유의함


** 그래프를 이용한 선형회귀모형 분석

use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_1.dta",clear

twoway (lfit food_exp income) (scatter food_exp income)
// lfit : 회귀직선, 적합선(linear fitting) 

twoway (lfitci food_exp income) (scatter food_exp income)
// lfitci : 회귀직선의 95% 신뢰구간 표시(linear fittting confidence interval)


** 관심변수 이외에 나머지 독립변수는 고정

use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_4.dta"

reg wage age ttl_exp tenure if union==1
// 노동조합에 가입한 근로자(union=1)의 근무연한(tenure)에 따른 임금(wage)에 관심

margins, at(tenure=(0(1)25)) atmeans noatlegend
estore(line1)
// margins, at(tenure=(0(1)25) : 근무연한(tenure, 독립변수)이 0~25일 때 임금(wage, 종속변수)의 예측값을 보여줌
// atmeans : 다른 변수는 평균값으로 고정  noatlegend : 공변량의 고정값을 표시하는 범례 없이  
/*

Adjusted predictions                                       Number of obs = 100
Model VCE: OLS

Expression: Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    7.45356   .6859693    10.87   0.000     6.091922    8.815199
          2  |   7.618791   .6237289    12.21   0.000     6.380699    8.856883
          3  |   7.784021    .564068    13.80   0.000     6.664355    8.903687
          4  |   7.949252   .5078964    15.65   0.000     6.941085    8.957418
          5  |   8.114482    .456504    17.78   0.000     7.208329    9.020635
          6  |   8.279713   .4116844    20.11   0.000     7.462526      9.0969
          7  |   8.444943   .3757968    22.47   0.000     7.698992    9.190894
          8  |   8.610174   .3515871    24.49   0.000     7.912279    9.308069
          9  |   8.775404   .3415477    25.69   0.000     8.097437    9.453371
         10  |   8.940634   .3469109    25.77   0.000     8.252022    9.629247
         11  |   9.105865   .3670022    24.81   0.000     8.377371    9.834359
         12  |   9.271095   .3996062    23.20   0.000     8.477883    10.06431
         13  |   9.436326   .4419622    21.35   0.000     8.559038    10.31361
         14  |   9.601556   .4915559    19.53   0.000     8.625826    10.57729
         15  |   9.766787     .54642    17.87   0.000     8.682152    10.85142
         16  |   9.932017   .6051227    16.41   0.000     8.730858    11.13318
         17  |   10.09725   .6666508    15.15   0.000     8.773956    11.42054
         18  |   10.26248   .7302904    14.05   0.000     8.812863    11.71209
         19  |   10.42771    .795535    13.11   0.000     8.848584    12.00683
         20  |   10.59294   .8620202    12.29   0.000     8.881843    12.30404
         21  |   10.75817   .9294798    11.57   0.000     8.913167    12.60317
         22  |    10.9234   .9977163    10.95   0.000     8.942949    12.90385
         23  |   11.08863    1.06658    10.40   0.000     8.971485    13.20578
         24  |   11.25386   1.135958     9.91   0.000     8.999002    13.50872
         25  |   11.41909   1.205761     9.47   0.000     9.025675    13.81251
         26  |   11.58432   1.275919     9.08   0.000     9.051643      14.117
------------------------------------------------------------------------------
*/
// _at 1은 tenure가 0, _at 26은 tenure가 25, margin은 각 tenure의 값에서 wage의 예측값


marginsplot, noci recast(line) addplot((scatter wage tenure if union==1)) legend(off)
// union=1인 값에 대해 wage를 y축, tenure를 x축으로 하는 산포도와 회귀선 그리기
// noci : 회귀을 중심으로 신뢰구간을 표시하지 않음, legend(off) : 범례를 표시하지 않음

marginsplot, recast(line) addplot((scatter wage tenure if union==1))


** 추정결과 표로 만들기

use "D:\STATA연습데이터\STATA기초통계와회귀분석\R_data8_5.dta"

reg lwg age wc hc inc if lwg>0
/*
      Source |       SS           df       MS      Number of obs   =       731
-------------+----------------------------------   F(4, 726)       =     34.52
       Model |     29.2399         4  7.30997501   Prob > F        =    0.0000
    Residual |  153.759787       726   .21179034   R-squared       =    0.1598
-------------+----------------------------------   Adj R-squared   =    0.1552
       Total |  182.999687       730  .250684503   Root MSE        =    .46021

------------------------------------------------------------------------------
         lwg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .0014282   .0021532     0.66   0.507     -.002799    .0056553
          wc |   .4286778   .0452609     9.47   0.000       .33982    .5175357
          hc |   .0271142   .0433457     0.63   0.532    -.0579837    .1122121
         inc |    .000249   .0015837     0.16   0.875    -.0028603    .0033583
       _cons |   .9531833   .0966086     9.87   0.000     .7635179    1.142849
------------------------------------------------------------------------------
*/

estimates table // 앞선 회귀결과의 변수와 추정계수를 표로 정리
/*
---------------------------
    Variable |   Active    
-------------+-------------
         age |  .00142817  
          wc |  .42867784  
          hc |  .02711422  
         inc |  .00024902  
       _cons |  .95318333  
---------------------------
*/

estimates table, b(%9.3f) stats(r2) stfmt(%9.2f) star(0.01 0.05 0.1)
// b(%9.3f) : 추정계수를 소수점 3자리까지 보여줌
// stats(r2) : ereturn list에 저장된 스케일러 중 r2 값을 불러옴
// stfmt(%9.2f) : stats()로 불러온 추정치의 포맷 지정
// star(0.01 0.05 0.1) : 유의수준 0.01 *** 0.05 ** 0.1 *
/*
-----------------------------
    Variable |    Active     
-------------+---------------
         age |     0.001     
          wc |     0.429***  
          hc |     0.027     
         inc |     0.000     
       _cons |     0.953***  
-------------+---------------
          r2 |      0.16     
-----------------------------
Legend: * p<.1; ** p<.05; *** p<.01
*/

reg lwg age wc hc inc if lwg>0
estimates store model1 // 회귀결과를 model1에 저장

reg lwg age wc hc if lwg>0
estimates store model2

// 저장된 회귀결과를 표로 불러옴  style(column) : 표의 각 열을 수직선으로 구분
estimates table model1 model2, stats(r2) stfmt(%9.2f) style(column)


reg lwg age wc hc inc if lwg>0

estout, cells(b(star fmt(%9.3f)) se(par fmt(%9.3f))) starlevels(* 0.1 ** 0.05 *** 0.01) legend  //se(par ) : 표준오차를 괄호안에 넣어서 표시
/*
----------------------------
                        .   
                     b/se   
----------------------------
age                 0.001   
                  (0.002)   
wc                  0.429***
                  (0.045)   
hc                  0.027   
                  (0.043)   
inc                 0.000   
                  (0.002)   
_cons               0.953***
                  (0.097)   
----------------------------
* p<0.1, ** p<0.05, *** p<0.01
*/


reg lwg age wc hc inc if lwg>0
estimates store reg1

reg lwg age wc hc if lwg>0
estimates store reg2

estout reg1 reg2 using D:\STATA연습데이터\STATA기초통계와회귀분석\table.xls , cells(b(star fmt(%9.3f)) se(par fmt(%9.3f))) starlevels(* 0.1 ** 0.05 *** 0.01) legend stats(r2, fmt(%9.2f))

 


[출처]  기초통계와 회귀분석(민인식, 최필선, 2012), 한국STATA학회 홈페이지(http://kastata.org/html/sub02-04.asp)

반응형

댓글