기본 콘텐츠로 건너뛰기

선형회귀 분석 (8)




 앞서 예제에서 선형 회귀 모델이 쉽게 적용되지 않는다는 것을 확인했습니다. 그러면 할 수 있는 일은 다른 방법을 모색하거나 비선형회귀 등 복잡한 모형을 도입해 보는 것입니다. 하지만 그 전에 해야 할 일은 어디서 잘못되었는지 알기 위해 데이터를 다시 뜯어보는 것입니다. 


 결국 많은 답이 데이터에 숨어 있으며 데이터에 대해서 잘 이해할 수록 정확한 결론을 도출해 낼 수 있습니다. 사실 이 경우에는 데이터 분포가 한쪽으로 치우친데다 다이아몬드 크기가 커질 수록 가격을 결정하는 다른 요소들이 많이 작용해서 단순한 예측 모델을 적용할 수 없는 것이 이유일 것입니다. 이를 확인하기 위해서 각 구간별 가격, 크기 분포를 확인해 보겠습니다. 


 다시 ggplot2 패키지를 읽어들인 후 이번에는 10개의 구간별로 carat 을 나눠 가격 분포를 다시 살펴보겠습니다. 앞서 예제에서 쓴 mygroup 함수를 다시 사용해 10개로 나누겠습니다. 

library("ggplot2")

mygroup<-function k="4){</span" y="">
   count=length(y)
   z=rank(y,ties.method = "min")
   return(floor((z-1)/(count/k))+1)
  }

diamonds$carat2<-mygroup carat="" diamonds="" span="">
table(diamonds$carat2)


 여기서 aggregate라는 매우 중요한 함수를 소개합니다. 이 함수를 이용해서 나만의 기능을 가진 다양한 연산 함수를 만들 수 있습니다. 


aggregate(price~carat2,data=diamonds,mean)

 이 코드는 새로 만든 carat2 그룹별로 평균을 구해 데이터프레임 형태로 보여달라는 이야기입니다. 

> aggregate(price~carat2,data=diamonds,mean)
   carat2      price
1       1   655.2940
2       2   750.7332
3       3   922.4984
4       4  1486.0546
5       5  2090.2473
6       6  3095.9929
7       7  5039.2091
8       8  5825.0127
9       9  7944.2637
10     10 12936.2783


 여기서 확장할 수 있는 것은 mean이라는 기능 대신 여러 개의 기능을 지닌 함수를 만들어 넣을 수 있다는 것입니다. myfx이라는 함수를 만들어 평균, 표준편차, 상한값, 하한값을 구해보겠습니다. 

myfx=function(x){
  result=c(mean(x),
           sd(x),
           min(x),
           max(x))
}

aggregate(price~carat2,data=diamonds,myfx)


 이렇게 하면 한 번에 여러가지 값을 구할 수 있으며 원하는 값만큼 더 포함시켜 구할 수 있습니다. 


> aggregate(price~carat2,data=diamonds,myfx)
   carat2    price.1    price.2    price.3    price.4
1       1   655.2940   173.7378   326.0000  2366.0000
2       2   750.7332   190.1553   345.0000  2346.0000
3       3   922.4984   215.3054   410.0000  1992.0000
4       4  1486.0546   430.2606   452.0000  4368.0000
5       5  2090.2473   621.2567   806.0000  6607.0000
6       6  3095.9929   814.0840   945.0000  9182.0000
7       7  5039.2091  1615.5915  1547.0000 16469.0000
8       8  5825.0127  1996.9402  1262.0000 18542.0000
9       9  7944.2637  2901.3972  2327.0000 18806.0000
10     10 12936.2783  3391.1394  3105.0000 18823.0000


 가격만 아니라 캐럿에 대해서도 같은 것을 구할 수 있습니다. 그리고 박스 플롯 등으로 시각화 할 수 있습니다. 


boxplot(price~carat2, outline=FALSE, data=diamonds)


> aggregate(carat~carat2,data=diamonds,myfx)
   carat2    carat.1    carat.2    carat.3    carat.4
1       1 0.29246125 0.02385143 0.20000000 0.31000000
2       2 0.33087712 0.01081135 0.32000000 0.35000000
3       3 0.39554695 0.01867721 0.36000000 0.42000000
4       4 0.49700353 0.03086971 0.43000000 0.53000000
5       5 0.62134302 0.06467172 0.54000000 0.70000000
6       6 0.78512112 0.07346134 0.71000000 0.90000000
7       7 0.98613037 0.03609363 0.91000000 1.01000000
8       8 1.06254319 0.03560108 1.02000000 1.13000000
9       9 1.31185889 0.13126791 1.14000000 1.51000000
10     10 1.85988565 0.30804080 1.52000000 5.01000000


 이렇게 보니 여러 가지를 확인할 수 있습니다. 1-3구간 까지는 선형으로 증가하던 가격이 4구간에서 한번 크게 증가하고 다시 7,9 구간에서 더 큰 폭으로 증가하는 것 같습니다. 이렇게 증가폭이 달라지는 만큼 하나의 선형 회귀 모델을 사용해서 회귀 계수를 얻는 방식이 잘 통하지 않는 것입니다. 쉽게 말해 기울기가 자꾸 변하는 그래프인 셈입니다. 


 여기서 한 가지 더 알아보기 위해 이번에는 10개의 그룹으로 균등하게 나누는 대신 특정값을 기준으로 나누는  cut breaks  함수를 사용해 보겠습니다. 0.5캐럿 이하, 0.5-1 캐럿, 1-2캐럿, 2캐럿 이상 이렇게 네 그룹으로 나누려면 어떻게 하는지 아래 함수에서 볼 수 있습니다. 5.01캐럿 다이아몬드까지 있으므로 중요한 점은 마지막 breaks 지점은 이것보다 크게 표시해야 한다는 것입니다. 그리고 각 그룹은 labels 함수로 이름을 정할 수 있습니다. 

diamonds$carat3=cut(diamonds$carat,breaks=c(0,0.5,1,2,6),labels = c(1,2,3,4))
table(diamonds$carat3)

 역시 carat3에 대해서도 같은 방식으로 알아보겠습니다. 


> aggregate(price~carat3,data=diamonds,myfx)
  carat3    price.1    price.2    price.3    price.4
1      1   839.7181   310.0636   326.0000  3378.0000
2      2  2811.3427  1302.8293   806.0000 16469.0000
3      3  7607.7241  3297.9883  1262.0000 18818.0000
4      4 14951.2504  2703.5945  5203.0000 18823.0000
> aggregate(carat~carat3,data=diamonds,myfx)
  carat3    carat.1    carat.2    carat.3    carat.4
1      1 0.35482569 0.06489307 0.20000000 0.50000000
2      2 0.72203473 0.15586754 0.51000000 1.00000000
3      3 1.25632870 0.24508433 1.01000000 2.00000000
4      4 2.15368449 0.23670424 2.01000000 5.01000000


이렇게 보면 몇 개의 그룹으로 나눠보면 선형 모델이 가능할 것 같습니다. 물론 꼭 이런 방법을 고집해야 하는 것은 아니지만, 한 번 해보겠습니다. subset 함수로 0.5 캐럿 이하 집단을 구해보겠습니다. 


diamonds1<-subset carat="" diamonds="" span="">
model=lm(price~carat, data=diamonds1)
summary(model)

par(mfrow=c(2,2))
plot(model)


> diamonds1<-subset carat="" diamonds="" span="">
> model=lm(price~carat, data=diamonds1)
> summary(model)

Call:
lm(formula = price ~ carat, data = diamonds1)

Residuals:
    Min      1Q  Median      3Q     Max 
-753.67 -145.77  -11.56  122.74 2040.33 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -377.353      8.721  -43.27   <2e-16 span="">
carat       3430.054     24.177  141.88   <2e-16 span="">
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 215.9 on 18930 degrees of freedom
Multiple R-squared:  0.5153, Adjusted R-squared:  0.5153 
F-statistic: 2.013e+04 on 1 and 18930 DF,  p-value: < 2.2e-16

> 



 이것도 플롯의 모양의 독특해보이지만, 이전 모델보다는 훨씬 나은 결과입니다. 참고로 줄무늬 처럼 나타나는 이유는 가격은 비교적 연속 분포지만, 캐럿은 사실 0.01단위로 끊어지는 변수이기 때문입니다. 이걸 보니 더 작게 그룹을 나누면 회귀 계수가 얼추 맞을 것 같습니다. 0.2-0.3캐럿 그룹을 구해봅시다. 


> diamonds1<-subset carat="" diamonds="" span="">
> model=lm(price~carat, data=diamonds1)
> summary(model)

Call:
lm(formula = price ~ carat, data = diamonds1)

Residuals:
    Min      1Q  Median      3Q     Max 
-336.71 -107.49   -4.88   97.90 1690.29 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -188.44      26.22  -7.186 7.87e-13 ***
carat        2880.53      92.29  31.213  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 149 on 4201 degrees of freedom
Multiple R-squared:  0.1883, Adjusted R-squared:  0.1881 
F-statistic: 974.3 on 1 and 4201 DF,  p-value: < 2.2e-16

> 
> par(mfrow=c(2,2))
> plot(model)


 이 정도면 몇 개 예측에서 크게 벗어나는 값 몇 개만 제외시키면 회귀 계수를 맞출 수 있을 것 같습니다. 회귀 계수는 이 구간에서 2880이므로 0.1캐럿 증가시 288달러 가격 증가로 볼 수 있습니다. 


 하지만 이것이 올바른 문제 해결일까요? 아닐 것입니다. 이렇게 잘게 잘라서 여러 개의 회귀 계수를 구하는 것이 과연 어떤 효용성이 있을까요? 오히려 10개 구간으로 잘라 가장 낮은 구간 대비 가장 높은 구간의 크기가 6배 크고 가격이 20배 비싸다는 쪽이 특별한 통계 분석이 없더라도 직관적으로 이해가 더 쉽습니다. 




 아무튼 이런 삽질(?)을 통해서 데이터의 특성을 더 잘 이해할 수 있게 되었고 처음 생각과는 다르게 단순하게 분석할 수 없다는 것도 알았습니다. 다음에는 샘플 숫자를 줄여 회귀 분석에 다른 예제들을 더 해보겠습니다. 부단한 삽질과 노력만이 실력이 느는 지름길일 것입니다. 

댓글

이 블로그의 인기 게시물

중등도 코로나 19 환자에서 제한적인 효과만을 확인한 렘데시비르

   현재 렘데시비르는 인공 호흡기 치료가 필요한 환자 등 중증 코로나 19 환자의 치료제로 사용되고 있습니다. 다만 공급이 부족해서 갑자가 중증 코로나 19 환자가 늘어날 경우 수급 문제가 불거질 수밖에 없습니다. 따라서 이를 중등도 코로나 19 환자 (경증이나 무증상은 당연히 치료가 필요없거나 대증 요법 정도) 사용한다는 것은 지금 당장에 급한 연구는 아니지만, 나중에 수급량이 많아지고 가격이 낮아질 때를 생각하면 의미 있는 연구 결과입니다.   렘데시비르 제조사인 길리어드 사이언스의 다아애나 브레이너드 박사 ( Dr. Diana Brainard of Gilead Sciences )는 병원에 입원이 필요했던 중등도 코로나 19 환자 584명을 렘데시비르 5일, 10일 투여 실험군과 위약을 투여한 대조군으로 나눠 임상 경과, 사망 등 중요한 합병증 발생 등에서 차이를 보이는지 조사했습니다. 연구 대상자의 40% 이상은 당뇨 같은 만성 질환이 있었으며 평균 연령은 57세였습니다.   연구 결과 대조군의 2%, 5일 투여군의 2%, 10일 투여군의 1%가 사망했으나 사망자 수가 4명, 3명, 2명 정도로 통계적 차이를 판단하기에는 사망자 수가 적어 결론을 내리기 어려웠습니다. 연구팀은 5일 투여군에서 임상 경과가 다소 호전되는 것을 판단했으나 10일 투여군에서는 통계적 차이를 검증할 수 없었습니다. 결론적으로 말해 중등도 코로나 19 환자에서 렘데시비르 투여는 확실한 임상적 이득이 없었습니다.     5일 투여량만 해도 가격이 3120달러라는 점을 감안할 때 중등도 환자에서 렘데시비르 투여는 설령 공급이 원활하다고 해도 적절한 치료법이 아닌 것으로 판단됩니다. 더구나 지금처럼 공급이 불안해 중증 환자 투여 물량 확보가 아쉬운 상황에서는 말할 것도 없습니다.  이번 연구는 렘데시비르보다 효과가 뛰어나면서 저렴한 차세대 항바이러스 치료제의 필요성을 다시 보여준 것으로 생각됩니다....

178Tb/s - 광섬유 통신 신기록 수립

  ( Lidia Galdino, lead researcher on the study. Credit: James Tye / UCL )  현재 가정에 들어오는 초고속 광통신은 10Gb/s 정도가 가장 빠른 편입니다. 그러나 사실 광섬유를 통해서 보낼 수 있는 정보량은 이것보다 훨씬 많습니다. 현재 개발 중인 기술 가운데서는 이미 1Tb/s 이상이 속도를 지닌 것들이 많아 미래에는 테라 인터넷 시대가 가능할 것으로 예상됩니다. 그런 가운데 영국 유니버시티 칼리지 런던 (UCL)의 리디아 갈디노( Lidia Galdino )가 이끄는 연구팀이 178Tb/s라는 신기록을 세웠다는 소식입니다.   광섬유 통신을 위한 주파수는 4.5THz까지 높일 수 있습니다. 최근에 개발된 신기술은 9THz 벽도 돌파했습니다. 그러나 연구팀은 16.8THz라는 새로운 영역에 도달했습니다. 주파수가 높을 수록 더 많은 정보를 전송할 수 있습니다. 하지만 16.8THz의 주파수로도 178Tb/s의 데이터 전송 속도는 사실 불가능합니다. 연구팀은  Geometric Shaping (GS) 방식으로 레이저 신호의 페이즈, 밝기, 편광 수준을 조절해 여러 신호가 서로 간섭 없이 전달되게 만들어 이런 신기록을 수립했습니다.   물론 새로운 기술이 항상 그렇듯이 바로 상용화가 가능한 것은 아닙니다. 특히 광섬유는 장거리 정보 전송이 가능해야 할 뿐 아니라 오류 없이 장시간 작동이 가능해야 합니다. 자꾸 에러가 생기거나 고장난다면 아무리 속도가 빨라도 무용지물입니다. 이 단계를 극복하고 경제성까지 갖춰 실용화 되는 경우는 사실 소수에 지나지 않습니다.   하지만 경쟁적으로 Tb 급 대역폭을 지닌 광통신 기술이 개발되고 있어 테라바이트 급 정보도 순식간에 전송할 수 있는 시대가 올 것은 분명합니다. 과거 전화선을 이용한 56K 모뎀에서 지금의 광섬유 기반 초고속 인터넷이 온 것처럼 속도는 계속 빨라질 수밖에 없기 때문입니다...

코로나 19 감염 위험성이 큰 동물은?

  ( A new genomic study ranks the potential of the SARS-CoV-2 spike protein to bind to the ACE2 receptor site in 410 vertebrate animals. Old World primates and great apes, which have identical amino acids at the binding site as humans, are predicted to have a very high propensity for binding ACE2 and are likely susceptible to SARS-CoV-2 infection. Credit: Matt Verdolivo/UC Davis )    코로나 19를 일으키는  SARS-CoV-2는 사스 및 메르스와 마찬가지로 박쥐를 자연 숙주로 삼는 베타 코로나바이러스 계통입니다. 따라서 본래 박쥐에 살던 바이러스가 중간 단계를 거쳐 사람으로 넘어온 것으로 추정됩니다. 그런데 이 이야기는 SARS-CoV-2 역시 사람에만 감염되는 것이 아니라 종을 넘어 다른 포유류나 조류, 혹은 척추동물에 감염될 수 있음을 시사합니다. 코로나 바이러스가 침투하는 경로인 ACE2 수용체는 다른 척추동물에서도 볼 수 있기 때문입니다. 다만 SARS-CoV-2와 ACE2 수용체의 결합 수준은 동물마다 큰 차이가 있어 감염 위험도는 모두 다를 수밖에 없습니다.   캘리포니아 대학 데이비스 캠퍼스 (UC Davis)의 조안나 다마스 (Joana Damas)와 그 동료들은 410종의 척추동물에서 ACE2 수용체 유전자를 조사해 코로나 19 감염 위험도를 조사했습니다. 연구팀에 따르면 ACE2 수용체 전체가 아니라 25개의 아미노산 배열이 특히 결합력에 큰 영향을 미쳐 감염 가능성을 결정합니다.   연구 결과 쉽게 예상할 수 있듯이 사람과 가까운 침팬치, 고릴라, 보노보...