[Data Mining with R Programming] 2006 Birth Data 분석

이번 글에서는 "Data for R in a nutshell"이 제공하는 2006년도 출산 데이터와 다양한 시각화 도구를 이용하여 데이터 분석 및 마이닝을 연습해 보도록 한다.

"R in a nutshell"의 데이터를 불러오는 방법에 대한 자세한 내용은 패키지와 데이터 기본을 참고하도록 한다.

---

[준비할 사항]

데이터의 시각화를 위해 lattice 라이브러리와 nuteshell 라이브러리를 불러온다. lattice 라이브러리는 barchart, histogram, density plot 등 다양한 데이터 시각화 도구를 제공하며, nuteshell 라이브러리는 데이터 분석 연습을 위한 샘플을 제공한다.

> library(lattice)
> library(nutshell)

데이터는 2006년도 출산 샘플 데이터인 births2006.smpl 을 불러오고, 차후 보다 간편하게 다루기 위해 변수 births에 저장한다:

> data(births2006.smpl)
> births = births2006.smpl

이 샘플 데이터는 총 427,323개의 데이터를 가지고 있으며 총 13개의 변수로 되어 있다:

> nrow(births)
[1] 427323
> ncol(births)
[1] 13

또는 다음과 같이 차원을 출력할 수 있다:

> dim(births)
[1] 427323     13

또한 데이터 타입은 data.frame이다:

> class(births)
[1] "data.frame"

우선 변수 이름을 살펴보도록 하자:

> names(births)
 [1] "DOB_MM"    "DOB_WK"    "MAGER"     "TBO_REC"   "WTGAIN"    "SEX"       "APGAR5"    "DMEDUC"    "UPREVIS"  
[10] "ESTGEST"   "DMETH_REC" "DPLURAL"   "DBWT" 

각 변수 이름이 의미하는 바는 다음과 같다:

DOB_MM

출산한 월.

DOB_WK

출산한 요일.

MAGER

산모의 나이.

TBO_REC

전체 출산 순서. 만약 첫째 출산이면 1, 둘째 출산이면 2.

WTGAIN

산모의 몸무게 변화.

SEX

태아의 성별. 여아일 경우 F, 남아일 경우 M.

APGAR5

AGPAR 점수. 아프가 점수라고 하며, 신생아의 상태를 평가하는 점수로서, A(Appearance)는 태아의 피부색, P(Pulse)는 맥박수, G(Grimace)는 반사흥분도, A(Activity)는 활동성, R(Respiration)은 호흡을 의미한다.

DMEDUC

산모의 교육 수준.

UPREVIS

태아기 방문수.

ESTGEST

임신 경과 주.

DMETH_REC

출산 방법. 제왕절개 및 자연분만.

DPLURAL

쌍생아 수. 1: 싱글, 2: 쌍둥이, 3: 세 쌍둥이, 4: 네 쌍둥이, 5: 다섯 쌍둥이 이상.

DBWT

신생아의 몸무게(그램).

데이터의 일부(1에서 5행까지)를 출력하여 확인해 보자:

> births[1:5,]
        DOB_MM DOB_WK MAGER TBO_REC WTGAIN SEX APGAR5                 DMEDUC UPREVIS ESTGEST DMETH_REC
591430       9      1    25       2     NA   F     NA                   NULL      10      99   Vaginal
1827276      2      6    28       2     26   M      9     2 years of college      10      37   Vaginal
1705673      2      2    18       2     25   F      9                   NULL      14      38   Vaginal
3368269     10      5    21       2      6   M      9                   NULL      22      38   Vaginal
2990253      7      7    25       1     36   M     10 2 years of high school      15      40   Vaginal
         DPLURAL DBWT
591430  1 Single 3800
1827276 1 Single 3625
1705673 1 Single 3650
3368269 1 Single 3045
2990253 1 Single 3827

---

[출산방법에 대한 요일별 분석]

샘플 데이터에서의 출산방법은 제왕절개(C-section)과 자연분만(vaginal)이며, 출산방법이 알려지지 않은(unknown) 데이터도 있다.

출산방법에 대한 요일별 분석을 하기 전에, 요일별 출산 빈도수를 분석해 보자.

R 프로그래밍에서 data.frame의 변수명에 대한 데이터를 출력하려면 "$변수명"을 붙인다. 즉, 출산요일 변수인 DOB_WK에 대하여 데이터를 출력하려면:

> births$DOB_WK

과 같이 입력한다. 이제 요일별 빈도수를 출력하려면 table() 함수를 사용한다:

> births.dow = table(births$DOB_WK)
> births.dow

    1     2     3     4     5     6     7 
40274 62757 69775 70290 70164 68380 45683 

숫자 1은 일요일, 2는 월요일, 3은 화요일, 4는 수요일, 5는 목요일, 6은 금요일, 7은 토요일을 의미한다.

바차트(barchart)를 이용하여 이를 시각화 해보자:

> barchart(births.dow, xlab="Day of Week", col="red", horizontal = FALSE)

이제 본격적으로 출산 방법에 대한 요일별을 분석을 해보자. 우선 DOB_WK에 대한 DMETH_REC을 테이블로 출력해 본다:

> dob.dm.tbl = table(WK=births$DOB_WK, DM=births$DMETH_REC)
> dob.dm.tbl
   DM
WK  C-section Unknown Vaginal
  1      8836      90   31348
  2     20454     272   42031
  3     22921     247   46607
  4     23103     252   46935
  5     22825     258   47081
  6     23233     289   44858
  7     10696     109   34878

위의 테이블을 보면 각 요일(WK)에 따른 출산방법의 빈도수를 확인할 수 있다. 출산방법이 알려지지 않은(Unknown) 열을 제거하려면 제거하고자 하는 열 번호에 (-) 부호를 붙인다:

> dob.dm.tbl = dob.dm.tbl[,-2]
> dob.dm.tbl
   DM
WK  C-section Vaginal
  1      8836   31348
  2     20454   42031
  3     22921   46607
  4     23103   46935
  5     22825   47081
  6     23233   44858
  7     10696   34878

위의 테이블을 바차트로 시각화 해보자.

> barchart(dob.dm.tbl, xlab="Day of Week", horizontal=FALSE)

푸른색으로 표현된 것이 제왕절개(C-section), 분홍색으로 표현된 것이 자연분만(vaginal)이다. 산부인과가 대체적으로 주말에 쉬는 것을 고려하면 주중에 횟수가 더 많은 것은 당연한 결과이다. 두 가지 출산방법에 대해 각각 나누어 표현해 보자. 이를 위해서는 table() 함수에 옵션 group=FALSE를 추가한다:

> barchart(dob.dm.tbl, xlab="Day of Week", horizontal=FALSE, group=FALSE, col=c("red","blue"))

위의 그림에서 왼쪽은 제왕절개, 오른쪽은 자연분만의 요일별 빈도이다. 일단 우리가 알 수 있는 것은 자연분만 횟수가 제왕절개에 비해 빈도수가 높다는 것이며, 출산은 주중에 더 많이 일어난다는 것이다. 이제 각각의 방법에 대해 요일별 비율로 변환해 보자:

> dob.dm.tbl[,1] = dob.dm.tbl[,1] / sum(dob.dm.tbl[,1])
> dob.dm.tbl[,2] = dob.dm.tbl[,2] / sum(dob.dm.tbl[,2])
> dob.dm.tbl
   DM
WK   C-section   Vaginal
  1 0.06690493 0.1067210
  2 0.15487476 0.1430901
  3 0.17355453 0.1586686
  4 0.17493261 0.1597853
  5 0.17282763 0.1602823
  6 0.17591695 0.1527143
  7 0.08098858 0.1187385

비율은 단순히 각각의 출산방법에 대해 요일별 횟수를 모두 더한 후 각 요일별 해당 출산방법으로 나눈 것이다.

이제 바차트를 통해 계산된 비율을 시각화 해보자:

> barchart(dob.dm.tbl, ylab="Day of Week", group=FALSE, horizontal=FALSE, col=c("red","blue"))

위의 그림을 살펴보면 제왕절개 방법이 자연분만 보다 주중과 주말의 횟수의 비율 차이가 현저하다. 이는 자연분만 보다 제왕절개가 아무래도 큰 수술이다 보니 주말에 시행되지 않았을 가능성이 크다고 추측해 볼 수 있다.

---

[출산방법 및 쌍생아 수에 대한 신생아의 몸무게 분석]

lattice 라이브러리가 제공하는 히스토그램(histogram)을 통해 출산방법(DMETH_REC)에 대한 신생아 몸무게 분석(DBWT)을 해보자.

우선 births의 DMETH_REC에 대해 Unknown을 제외한 데이터를 births2로 저장하자:

> births2 = births[births$DMETH_REC!="Unknown",]

다음 명령은 신생아의 몸무게를 출산방법 별로 분류하여 히스토그램을 나타낸 것이다:

> histogram(~DBWT|DMETH_REC, data=births2, layout=c(1,2), col=c("red","blue"))

몸무게의 단위는 그램이다. 위의 그림을 보면 신생아의 몸무게는 대략 2000 - 4000 그램 사이에 집중되어 있는 것을 확인할 수 있다. 그러나, 출산방법의 차이에 따른 신생아의 몸무게 분포의 뚜렷한 차이는 없어 보인다. 즉, 출산방법과 신생아의 몸무게의 상관관계는 없는 것으로 결론 내릴 수 있다.

이번에는 쌍생아 수(DPLURAL)에 따른 신생아의 몸무게(DBWT) 분포를 확인해 보자. 이에 앞서 우리는 이미 쌍생아 수의 데이터 종류가 5개(Single / Twin / Triplet / Quadraplet / Quintuplet or Higher)임을 알고 있다. 그런데 만약 이를 모르고 있는 경우 알아내는 방법에 대해 생각해 보자.

split() 명령을 이용하여 DPLURAL에 따라 데이터를 분류한 후 이에 대한 길이를 알아내는 length() 명령을 이용하면 DPLURAL의 데이터 종류의 개수를 알아낼 수 있다:

> length(split(births, births$DPLURAL))
[1] 5

쌍생아 수(DPLURAL)에 따른 신생아의 몸무게(DBWT) 분포를 확인하는 명령은 다음과 같다:

> histogram(~DBWT|DPLURAL, data=births, layout=c(1,5), col=c("red","blue"))

아래로부터 각각 싱글, 두 쌍둥이, 세 쌍둥이, 네 쌍둥이, 다섯 쌍둥이 및 그 이상에 대해 신생아의 개별 몸무게 분포를 나타내고 있다. 쌍생아 수가 증가할 수록 개별 신생아의 몸무게는 대체적으로 낮아지고 있음을 확인할 수 있다.

lattice 라이브러리가 제공하는 densityplot을 통해 좀 더 세밀하게 데이터를 시각화해보자:

> densityplot(~DBWT|DPLURAL, data=births, layout=c(1,5), plot.points=FALSE)

네 쌍둥이까지는 비교적 스무스한 형태를 띠고 있으나 다섯 쌍둥이 이상은 그래프가 마치 노이즈처럼 요동치는 것을 볼 수 있다. 이는 다섯 쌍둥이가 통계적 자료로 사용하기에는 매우 드문 경우이기 때문이다.

이번에는 쌍생아의 몸무게 평균을 성별(SEX)로 나누어 분석해 보자.

우선 쌍생아 수에 대한 요인(factor)을 추출한다:

> fac = factor(births$DPLURAL)

다음으로 DBWT에 해당하는 데이터를 추출한다:

> res = births$DBWT

데이터 획득 시 언제나 모든 데이터가 완벽하게 획득되지는 않는다. 매우 자주 유실 데이터(missing values)가 존재하기 마련인데, 유실 데이터는 NA로 표현되며 is.na() 함수를 통해 다음과 같이 유실 데이터가 존재 유무를 확인할 수 있다:

> sum(is.na(res))
[1] 434

원리는 간단하다. is.na(res)는 NA인 경우 FALSE를, 이외의 경우 TRUE를 출력하는데 FALSE는 "0"을 TRUE는 "1"의 값을 갖기 때문에 합계를 계산하는 sum()이 1 이상의 값을 출력하면 NA, 즉 유실 데이터가 존재하며 0을 출력하면 유실 데이터가 없는 완전한 데이터 세트를 의미한다.

위에서 보듯이 총 434개의 유실 데이터가 존재한다. 데이터 분석 시 유실 데이터는 제외되어야 한다. 이에 대해서는 곧 다시 언급하기로 하고, R은 어떠한 반복 계산을 일괄처리는 함수들이 존재하는데 apply 함수들이다. 이에 대한 함수군은 lapply, apply, tapply, mapply이며 목적은 비슷하나 쓰임새가 약간씩 다르다. 각각에 대한 자세한 내용은 R의 헬프를 통해 알아보기 바란다. 일단 이들은 어떠한 그룹에 대해 각각 어떠한 함수로 일괄처리하는 역할을 한다.

다음 예를 살펴보자:

> x = list(a = 1:5, b = 6:10, c = 11:13)
> x
$a
[1] 1 2 3 4 5

$b
[1]  6  7  8  9 10

$c
[1] 11 12 13

> lapply(x, mean)
$a
[1] 3

$b
[1] 8

$c
[1] 12

변수 x는 벡터변수 a, b, c 세 개로 이루어진 list형이다. 위에서 보듯 lapply의 첫번째 인자는 x이며, 두번째 인자는 평균을 구하는 mean() 함수이다. 함수를 인자로 취하는 것에 익숙하지 않은 분들도 계시겠지만 이와 같이 함수를 인자로 받는 경우를 함수 핸들러라고 표현하기도 한다. 눈치챘겠지만 lapply에 의해 계산된 결과는 a, b, c에 대한 각각의 평균값이다. 결국 lapply는 각 그룹에 대한 어떠한 함수를 일괄적으로 처리하는 함수이며, tapply도 이와 유사하다.

다음 명령을 살펴보자:

> tmp = tapply(res, fac, mean, na.rm=TRUE)
> tmp
              1 Single                 2 Twin              3 Triplet           4 Quadruplet 
              3298.263               2327.478               1677.017               1196.105 
5 Quintuplet or higher 
              1142.800 

앞서 res는 births의 DBWT에 대한 그룹화를, fac는 DPLURAL에 대한 요인을 정의한 바 있다. 위의 명령은 DBWT에 대해 DPLURAL의 요인을 그룹화하여 각 그룹에 대한 평균(mean())을 구하되 유실 데이터는 제거(na.rm=TRUE)하고 구하도록 한 것이다.

Barchart를 통해 이를 시각화하려면 다음과 같이 입력한다:

> barchart(tmp, horizontal=FALSE, ylab="DPLURAL", col="cyan")

위의 그림에서 보듯이 싱글, 두 쌍둥이, 세 쌍둥이, 네 쌍둥이, 다섯 쌍둥이 및 그 이상으로 갈수록 신생아의 평균 몸무게가 낮아지는 것을 확인할 수 있다.

이번에는 성별까지 분류대상으로 포함시켜 보자:

> fac1 = births$DPLURAL
> fac2 = births$SEX
> tmp = tapply(res, INDEX=list(fac1, fac2), FUN=mean, na.rm=TRUE)
> tmp
                              F        M
1 Single               3242.302 3351.637
2 Twin                 2279.508 2373.819
3 Triplet              1697.822 1655.348
4 Quadruplet           1319.556 1085.000
5 Quintuplet or higher 1007.667 1345.500

성별에 대한 비교를 쉽게 할 수 있도록 Barchart 대신 Barplot을 사용하여 다음과 같이 시각화하였다:

> barplot(tmp, beside=TRUE, ylab="DBWT")

여아와 남아 모두 쌍생아의 크기가 커질수록 평균 몸무게는 낮아지는 결과를 보이고 있다. 다만 다섯 쌍둥이 및 그 이상은 매우 드문 경우라서 남아의 경우 일종의 데이터 노이즈와 같은 데이터를 보이고 있다. 충분한 샘플이 확보되면 아마도 네 쌍둥이와 비슷하거나 낮게 나올 것으로 생각된다.

---

[산모의 몸무게 변화와 신생아의 몸무게 관계 분석]

이번에는 산모가 임신 중 증가한 체중과 신생아의 몸무게 관계에 대해 분석하고자 한다. 단순하게 생각하면 산모의 임신 중 증가한 체중이 신생아의 몸무게와 상관관계가 있다고 판단할 수 있다. 과연 그러한지 분석을 통해 알아보자.

이번에는 lattice 라이브러리가 제공하는 xyplot을 이용하여 경향을 파악해 보기로 한다. xyplot은 scatter plot이라고도 하는데 해당 데이터의 값을 점(dot)으로 찍어 표시하는 것이다. xyplot은 점을 모두 찍어야하기 때문에 다른 플롯에 비해 시간이 많이 소요된다.

만약 각 요일별 신생아의 몸무게에 대한 scatter plot을 표현하려면 다음과 같이 입력한다:

> xyplot(DBWT~DOB_WK,data=births,col="black")

누구나 예상할 수 있겠지만 신생아가 태어난 요일과 신생아의 몸무게에는 별다른 상관관계가 없는 것으로 판단된다. 다만 화요일(DOB_WK = 2)에는 유독 아웃라이어(outlier)들이 많이 보인다.

신생아의 몸무게와 산모의 임신 중 몸무게 변화에 대한 xyplot을 만들어 보자:

> xyplot(DBWT~WTGAIN,data=births,col="red")

그래프를 살펴보면 신생아의 몸무게와 산모의 임신 중 증가한 몸무게 변화 사이에는 뚜렷한 경향성이 없어보인다. 상관관계를 보다 잘 판단하려면 xyplot 보다 smoothScatter가 더 도움이 될 수 있을 것 같다:

> smoothScatter(births$WTGAIN, births$DBWT, xlab="WTGAIN", ylab="DBWT")

smoothScatter은 상관관계를 보다 명확하게 판단할 수 있도록 할 뿐만아니라 그래프를 작성하는데 걸리는 시간도 xyplot 보다 훨씬 짧다.

진하게 표현된 부분이 거의 평행선에 가깝기 때문에 두 변수 사이의 상관관계는 없는 것으로 판단할 수 있다.

그러면 혹시라도 각 쌍생아 수에 대해 두 변수 사이의 상관관계가 있을까?

> xyplot(DBWT~WTGAIN|DPLURAL,data=births,layout=c(1,5), col="red")

마찬가지로 두 변수 사이의 뚜렷한 상관관계는 찾을 수 없다.

---

[신생아의 몸무게와 APGAR 관계 분석]

앞의 예제와 달리 신생아의 상태를 평가하는 AGPAR과 신생아의 몸무게 사이에는 강한 상관관계가 있을 것으로 예상된다. 왜냐하면 신생아의 몸무게가 낮으면 분명 AGPAR 점수는 낮기 때문이다.

이번 예제에서는 분석결과의 시각화를 위해 boxplot을 사용할 것인데 이에 대해 간단히 알아보도록 하겠다.

Boxplot은 샘플 미디언(median), 상위 및 하위 샘플 쿼타일(quartile), 최대 최소 데이터 값을 표현하기 위한 표현방법이다. 미디언 및 쿼타일에 대한 개념이 궁금하신 분들은 통계학 관련 사이트나 서적을 검색해 보기 바란다.

Boxplot을 통해 AGPAR 점수와 DBWT와의 관계를 시각화 해보자.

> boxplot(DBWT~APGAR5, data=births, ylab="DBWT", xlab="APGAR5", col="red")

그림에서 보듯이 왼쪽 부분의 한 두 개를 제외하고 APGAR과 DBWT는 양의 상관관계가 있음을 알 수 있다.

이번에는 DBWT와 DOB_WK 간의 관계가 있는지 분석해 보자. 물론 태어난 요일과 신생아의 몸무게 사이에는 어떠한 특별한 관계가 있을 것 같지 않다.

다음과 같이 입력해 보자:

> boxplot(DBWT~DOB_WK, data=births, ylab="DBWT", xlab="Day of Week", col="red")

예상대로 둘 사이에는 특별한 상관관계가 존재하지 않는다.

이번에는 lattice 라이브러리가 제공하는 bwplot을 활용하여 요인(factor)별로 결과를 시각화 해보자. "요인" 또는 "인자"라는 것은 주로 실험계획 및 분석법(Design of Experiments, DoE)에서 어떠한 요소가 시스템의 성능에 영향을 미치는지, 미친다면 양의 상관관계인지 또는 음의 상관관계인지를 파악하는데 사용된다. 그래서 실험 시 성능에 영향을 미칠 것 같은 요인과 이에 대한 수준(또는 레벨)을 정의한다. 이를 요인 분석(Factorial Design)이라고 한다.

가령 APGAR과 DBWT 사이의 관계를 분석할 때 성별(SEX)에 따라 다른 관계가 나타나는지를 파악한다고 가정해 보자. 이 때 성별이 요인이 된다.

다음과 같이 입력하여 요인 분석을 해보자:

> bwplot(DBWT~factor(APGAR5)|factor(SEX), data=births, xlab="APGAR5")

여아(F)와 남아(M)의 경우 모두 APGAR과 DBWT 사이에 대체적으로 양의 상관관계를 보이며, 성별 간에 특이점은 없어 보인다.

---

[산모의 체중 변화와 임신 기간 관계 분석]

분석에 앞서 table() 명령을 다시 한 번 이해하도록 하자. 이 명령은 해당 요인에 대한 수준(level)별 빈도수를 출력한다. 예를 들면 다음과 같다:

> table(births$DPLURAL)

              1 Single                 2 Twin              3 Triplet           4 Quadruplet 
                412979                  13658                    642                     39 
5 Quintuplet or higher 
                     5 

ESTGEST는 산모의 임신기간(Weeks)을 의미하며 이 인자에 대해 table() 명령을 통해 수준별 빈도수를 확인해 보자:

> table(births$ESTGEST)

    12     15     17     18     19     20     21     22     23     24     25     26     27     28     29 
     1      2     18     43     69    116    162    209    288    401    445    461    566    670    703 
    30     31     32     33     34     35     36     37     38     39     40     41     42     43     44 
  1000   1243   1975   2652   4840   7954  15874  33310  76794 109046  84890  23794   1931    133     32 
    45     46     47     48     51     99 
     6      5      5      2      1  57682 

가령 임신 12주는 1회, 15주는 2회 등등이다. 예상대로 임신 10개월 정도인 39주 또는 40주에 빈도수가 가장 많다.

Barplot을 통해 빈도수를 시각화 해보자:

> barplot(table(births$ESTGEST), ylab="Frequency", xlab="Estimated Gestation")

위의 그림을 보면 가장 오른쪽 빈도수가 뜬금없이 많이 있음을 알 수 있다. 99에 대한 데이터인데 사실 상식적으로 임신 99주라는 것이 도저히 납득이 안 된다. 이것은 사실 임신 기간이 알려지지 않은 데이터를 99에 표현한 것이다. 따라서 알려져 있는 의미있는 데이터만을 취급하려면 99주에 대한 데이터를 제거해야 한다:

> new_births = births[births$ESTGEST != 99, ]

위의 명령을 통해 임신시간이 알려져있지 않은(99) 데이터를 포함하는 행을 제거하였다. 제거가 잘 되었는지 확인해 보자:

> table(new_births$ESTGEST)

    12     15     17     18     19     20     21     22     23     24     25     26     27     28     29 
     1      2     18     43     69    116    162    209    288    401    445    461    566    670    703 
    30     31     32     33     34     35     36     37     38     39     40     41     42     43     44 
  1000   1243   1975   2652   4840   7954  15874  33310  76794 109046  84890  23794   1931    133     32 
    45     46     47     48     51 
     6      5      5      2      1