상관분석, 회귀분석의 자동화(part 1)

[cardiomoon]

(1) 들어가는 말

임상에서 실제로 다루는 데이타는 상당히 방대한 경우가 많다. 일변량인 경우는 거의 없다고 봐야 하며 웬만한 경우에도 변수들이 수십개가 되는 경우도 허다하다.

R에서 상관분석, 회귀분석은 비교적 쉬운 편이지만 하나하나의 자료를 상관분석, 회귀곡선을 통해 회귀곡선을 구하고 멱승, 로그, 지수 변환에 의한 재표현이 필요한 경우가 많으며 또한 잔차분석등을 통하여 알맞은 모델을 찾으려는 노력이 필요하다.  

또한 두 변수의 상관관계가 단순한 직선이 아닌 경우도 많으며 조각별 회귀, 산점도 평활기 등을 이용해 적절한 회귀곡선을 찾아야 할 경우도 있을 것이다.

하지만 이번 장에서는 상관분석 및 단순회귀분석을 자동화하여 여러 변수가 있는 다변량 자료에서 서로 상관관계가 있는 이변량분석을 빨리 할 수 있게 자동화하는 것이 목표이다.

먼저 사용할 데이타는 R에 내장되어 있는 attitude자료를 사용한다.

attitude 자료는 7개의 변수가 모두 연속형자료이다.

> tail(attitude)

   rating complaints privileges learning raises critical advance

25     63         54         42       48     66       75      33

26     66         77         66       63     88       76      72

27     78         75         58       74     80       78      49

28     48         57         44       45     51       83      38

29     85         85         71       71     77       74      55

30     82         82         39       59     64       78      39

7개의 변수의 상관관계를 한눈에 보려면 pairs를 사용하면 된다.

> pairs(attitude)

(2) 전통적인 R의 방법

하지만 한 쌍의 변수들, 예를 들어 attitude 데이타에 있는 rating과 complaints의 상관관계 

r값 및 회귀식, 그래프를 그리려면 여러 줄에 걸친 수고를 해야 한다. 즉, 먼저 r값을 구하기 위해 cor.test()를 시행한다.

> attach(attitude)

The following object(s) are masked from 'attitude (position 3)':

    advance, complaints, critical, learning, privileges, raises, rating

> cor.test(rating,complaints)

Pearson's product-moment correlation

data:  rating and complaints 

t = 7.737, df = 28, p-value = 1.988e-08

alternative hypothesis: true correlation is not equal to 0 

95 percent confidence interval:

 0.6620128 0.9139139 

sample estimates:

      cor 

0.8254176 

여기서 우리가 r값은 0.825, p값은 1.988e-08로 매우 낮다는 것을 알수 있다. 또한 회귀식을 구하기 위해 단순회귀분석 및 그래프를 그리면

> out=lm(rating~complaints,data=attitude)

> summary(out)

Call:

lm(formula = rating ~ complaints, data = attitude)

Residuals:

     Min       1Q   Median       3Q      Max 

-12.8799  -5.9905   0.1783   6.2978   9.6294 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept) 14.37632    6.61999   2.172   0.0385 *  

complaints   0.75461    0.09753   7.737 1.99e-08 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.993 on 28 degrees of freedom

Multiple R-squared: 0.6813, Adjusted R-squared: 0.6699 

F-statistic: 59.86 on 1 and 28 DF,  p-value: 1.988e-08 

> plot(rating~complaints,data=attitude,col="blue")

이상에서 회귀곡선의 기울기는 0.75, y 절편은 14.38이라는 것을 알 수 있었고 그래프를 그리기 위해 다음의 코드를 입력하면

> abline(out,col="red")

다음과 같은 그래프를 얻을 수 있다.

이와 같은 과정을 한꺼번에 할수 있는 함수를 만들면 편할 것이다. 즉 한 줄만 입력하면 자동으로 r 값과 p값을 구해주고 그래프까지 그려주면 좋지 않을까? 그 함수의 이름을 임의로 my.lm2.test라고 이름 지어보자. 이 함수가 하는 일은 p값과 r값, 회귀식의 기울기 및 y절편을 구해주고 그림을 그린다. 서비스로 회귀식을 그래프의 제목으로 써 주고 p값과 r값을 그래프에 표시해준다.

> my.lm2.test("rating","complaints",attitude,drawgraph=1)

(3) 자동화 함수를 만들자

이제 본격적으로 pairs함수에 도전해보자. my.lm.test함수는 두개의 변수들의 배열을 받아들여 p값 및 r 값, 기울기, y절편의 배열의 list를 반환한다.그래프를 그리라고 하면 그래프도 그려준다.

예를 들어 아까 attitude 자료를 이용하여 my.lm.test를 실행하면 다음과 같다.

> x=names(attitude)

> y=x

> result=my.lm.test(y,x,attitude)

> result

$r

           rating complaints privileges learning raises critical advance

rating      1.000      0.825      0.426    0.624  0.590    0.156   0.155

complaints  0.825      1.000      0.558    0.597  0.669    0.188   0.225

privileges  0.426      0.558      1.000    0.493  0.445    0.147   0.343

learning    0.624      0.597      0.493    1.000  0.640    0.116   0.532

raises      0.590      0.669      0.445    0.640  1.000    0.377   0.574

critical    0.156      0.188      0.147    0.116  0.377    1.000   0.283

advance     0.155      0.225      0.343    0.532  0.574    0.283   1.000

$p

           rating complaints privileges learning raises critical advance

rating      0.000      0.000      0.019    0.000  0.001    0.409   0.413

complaints  0.000      0.000      0.001    0.001  0.000    0.321   0.233

privileges  0.019      0.001      0.000    0.006  0.014    0.438   0.063

learning    0.000      0.001      0.006    0.000  0.000    0.542   0.003

raises      0.001      0.000      0.014    0.000  0.000    0.040   0.001

critical    0.409      0.321      0.438    0.542  0.040    0.000   0.129

advance     0.413      0.233      0.063    0.003  0.001    0.129   0.000

$slope

           rating complaints privileges learning raises critical advance

rating       1.00       0.75       0.42     0.65   0.69     0.19    0.18

complaints   0.90       1.00       0.61     0.68   0.86     0.25    0.29

privileges   0.43       0.51       1.00     0.51   0.52     0.18    0.41

learning     0.60       0.53       0.47     1.00   0.72     0.14    0.61

raises       0.50       0.52       0.38     0.57   1.00     0.40    0.58

critical     0.13       0.14       0.12     0.10   0.36     1.00    0.27

advance      0.13       0.17       0.29     0.47   0.57     0.29    1.00

$intercept

           rating complaints privileges learning raises critical advance

rating       0.00      14.38      42.11    28.17  19.98    50.24   56.76

complaints   8.24       0.00      34.32    28.44  11.21    47.71   54.12

privileges  25.45      18.97       0.00    24.15  19.25    39.52   35.61

learning    17.50      21.33      31.22     0.00   9.65    46.08   30.33

raises      32.05      29.83      44.52    32.66   0.00    35.02   39.72

critical    66.55      65.48      68.44    69.26  51.58     0.00   63.07

advance     34.46      31.38      27.60    16.67   6.21    20.91    0.00

즉, 가로7개 * 세로 7개의 49쌍의 결과가 한꺼번에 나오게 된다. 이 결과를 이용해 보기 좋게 r값과 p값을 나타내주는 함수를 하나 만들었다.

> my.lm.print(result)

만일 회귀식까지 원한다면 다음과 같이 사용하면 된다.

> my.lm.print(result,mode=1)

이상으로 상관분석, 단순회귀분석의 자동화 함수를 만들었고 다음에는 다중회귀분석에 도전해 보기로 하겠다. 

[cardiomoon]

제 글이 도움이 되셨다면 댓글 부탁드립니다. 혹시 마음에 안드신다면 질책 부탁드립니다.

[안재형]

저두 가끔 자동화에대해서 질문을 받는데 그런분들에게 도움이 되리라 생각됩니다.

[남현수]

어차피 LaTex 또는 외부파일로 지접 출력되는것이 아니라면 지난번 카이스퀘어에 부분에서 보여주셨던 라인있는 결과물보다 가져다 쓰기 좋아 보입니다.

[cardiomoon]

외부화일로 출력하는 것은 언제든지 함수출력 전 sink("화일이름")하시고 함수출력한 후 sink()하시면 되는데요.
LaTex도 출력할 수 있지만 필요할까요? 필요하시다면 할수 있읍니다.

[남현수]

제가 외부화일이라한것은 PDF 나 이미지 화일처럼 완성된 형태로 그냥 보고서에 붙여쓸수 있는것을 말씀드린겁니다. LaTex 같은 경우도 논문전체를 Latex로 쓰지 않는 이상 현실적이지는 않다고 생각되어집니다. 어차피 그냥 택스트상태로 편집해서 쓸거라면 라인이 없는 편이 편집하기가 좋다고 말씀드린 겁니다.

[우창윤]

좋은 자료 감사드립니다^^