Data Preparation
Seo Young Park, PhD
\(~\)
R을 이용해서 실제로 분석에서 데이터를 읽어들이고 분석 가능한 형태로 만드는 과정에 대해 살펴보겠습니다.
그런데, 실제 데이터 분석에서 가장 문제를 많이 일으키고 우리를 당황시키는 것이 R의 data type 입니다. 아주 간단하게 꼭 필요한 개념만 잠시 살펴보고 넘어가겠습니다.
\(~\)
Playing around with the R data type
값이 1, 2, 3, 4, 5 인 벡터 세가지를 만들어보겠습니다. aa는 numeric, bb는 factor, cc는 character 입니다.
aa<-c(1, 2, 3, 4, 5)
bb<-as.factor(aa)
cc<-as.character(aa)aa, bb, cc 세가지 벡터가 어떻게 저장되어있는지, 콘솔에 벡터이름을 쳐서 확인해봅시다.
aa## [1] 1 2 3 4 5bb## [1] 1 2 3 4 5
## Levels: 1 2 3 4 5cc## [1] "1" "2" "3" "4" "5"차이가 보이시나요?
\(~\)
실제 분석에서 type은 중대한 영향을 끼칩니다. 이 벡터들에 1씩을 더해보도록 하겠습니다.
aa+1## [1] 2 3 4 5 6bb+1## Warning in Ops.factor(bb, 1): '+' not meaningful for factors## [1] NA NA NA NA NAcc+1## Error in cc + 1: non-numeric argument to binary operator이러니, numeric type의 벡터가 들어가야하는 함수에 factor나 character를 넣으면 안되겠죠. 항상 함수가 input에 요구되는 data type을 잘 보고 꼭 맞게 넣어주어야 합니다. 그럼, 어떤 벡터의 type이 무엇인지는 어떻게 알 수 있을까요? typeof() 함수를 쓰면 됩니다.
typeof(aa)## [1] "double"typeof(bb)## [1] "integer"typeof(cc)## [1] "character"그런데 분석에 쓰이는 모든 변수를 이런식으로 일일히 typeof() 함수에 넣어주기는 번거로울 수 있어요. 일단 분석에 쓰이는 데이터셋처럼, 이 세개의 벡터를 data frame으로 묶어줍시다. 그리고 summary() 함수를 쓰시면, 한눈에 모든 변수들의 type을 짐작할 수 있습니다.
dd<-data.frame(aa, bb, cc)
summary(dd)## aa bb cc
## Min. :1 1:1 Length:5
## 1st Qu.:2 2:1 Class :character
## Median :3 3:1 Mode :character
## Mean :3 4:1
## 3rd Qu.:4 5:1
## Max. :5뿐만 아니라, numeric type 벡터(연속 변수)의 경우에는 최소값, 최대값 등등과, factor type 벡터(범주형 변수)의 경우에는 각 범주별 카운트까지 확인할 수 있습니다. 따라서 분석 전 데이터 정리의 목적은, 분석에 쓰일 data frame을 summary() 함수에 넣었을 때, type과 값들이 이상한 점이 없도록 고쳐주는 것이라고 생각하시면 됩니다.
\(~\)
dplyr package and pipe operator (%>%)
이 강의에서는 분석과 자료정리에 dplyr 패키지를 적극 사용하려고 합니다. 처음에는 어려워보이지만 익숙해지면 훨씬 쉽고 빠르게 자료를 정리할 수 있습니다. 이 패키지는 보다 효율적인 데이터분석을 위해 Hadley Wickham이라는 통계학자가 만들고 있는 tidyverse의 일부입니다. 궁금하신 분들은 https://www.tidyverse.org/ 를 참고하세요. 일단 필요한 패키지들을 설치합니다.
library(plyr)
library(dplyr)
library(tableone)dplyr 패키지에는 파이프 연산자(%>%) 라는 것이 있습니다. R에서 하는 많은 작업은 어떤 객체를, 어떤 함수 안에 입력값으로 넣어서 결과값을 얻는 식으로 이루어집니다. 이런 작업의 명령문은 함수(객체)의 형식을 취합니다. 이것을 파이프 연산자를 사용해서 다시 쓰면, 객체 %>% 함수() 형식이 됩니다. 예를 들어 1부터 10까지의 자연수를 더하기 위한 전통적인 R 코드는 아래와 같이 됩니다.
x<-1:10
sum(x)## [1] 55같은 작업을 파이프 연산자를 이용해서 쓰면 아래와 같이 됩니다.
x<-1:10
x %>% sum()## [1] 55여러 단계를 연속적으로 수행해야하는 경우를 생각해봅시다. 1부터 10까지 각 자연수에 자연로그를 취하고 소수점 첫째자리에서 반올림한 후 모두 더한 값을 구하려면 어떻게 해야할까요? 아래 두 줄의 코드는 정확히 같은 작업을 수행하는 프로그램입니다.
sum(round(log(x)))## [1] 15x %>% log() %>% round() %>% sum() ## [1] 15실제로 x에 적용되는 함수는 log, round, sum 의 순서인데, 전통적인 방법으로 코딩하면 명령문의 왼쪽에서부터 읽을 때 sum, round, log 순서로 함수가 등장합니다. 반면 파이프 연산자를 사용해서 코딩하면 실제로 적용되는 순서로 함수가 등장하니, 사람이 명령문을 읽고 이해하기 쉬워집니다.
자, 이제 파이프 연산자의 사용법을 알았으니, 본격적으로 데이터를 정리해봅시다.
\(~\)
Importing the data and checking for duplicates
예제 데이터 example_data.csv를 컴퓨터에 저장한 후 RStudio를 켜기도 전에 해야할 일이 한가지 있습니다. 무엇일까요? (정답은 강의에서 공개하겠습니다!)
RStudio를 켠 후, 각자 데이터가 저장된 위치 경로를 setwd 안에 넣어서 working directory를 셋업합니다. 아래와 같이 역슬래시를 두개씩 붙이거나,
setwd("C:\\Users\\SYP\\Dropbox\\Seminars\\R_data_prep")아래와 같이 오른쪽으로 누운 슬래시를 한개씩 붙이셔도 됩니다.
setwd("C:/Users/SYP/Dropbox/Seminars/R_data_prep")이제 데이터를 read.csv() 함수를 이용해 읽어들입니다.
dat0<-read.csv("example_dataset.csv", stringsAsFactors = F)데이터를 이렇게 읽어들인 이후 꼭 해야하는 일이 있습니다. 무엇일까요? 그 일을 하고 나면 아래와 같이 코드를 고쳐야 한다는 것을 알게 됩니다.
dat0<-read.csv("example_dataset.csv", stringsAsFactors = F, nrow=50)그리고 나서 꼭 해야하는 중요한 일이 있습니다. 각 행이 unique한지 확인하는 것이죠. 일단 50개의 환자 아이디가 unique한지 확인합니다.
length(unique(dat0$id))## [1] 50유니크한 환자 아이디는 50개 입니다. 즉 50명의 환자가 50개의 행으로 들어갔으니 중복된 아이디가 없다는 뜻입니다.
\(~\)
Correcting the types of the variables
자, 이제 summary()에 dat0를 넣어서 이상한 점이 있나 확인해봅시다.
summary(dat0)## id age sex Recur..1..Recur..0..Censored.
## Min. : 1.00 Min. :25.00 Min. :0.00 Min. :0.00
## 1st Qu.:13.25 1st Qu.:49.25 1st Qu.:0.00 1st Qu.:0.00
## Median :25.50 Median :58.50 Median :1.00 Median :0.00
## Mean :25.50 Mean :56.44 Mean :0.66 Mean :0.44
## 3rd Qu.:37.75 3rd Qu.:65.75 3rd Qu.:1.00 3rd Qu.:1.00
## Max. :50.00 Max. :79.00 Max. :1.00 Max. :1.00
## OP_date Recur_date local_6m CA19.9
## Length:50 Length:50 Min. :0.00 Length:50
## Class :character Class :character 1st Qu.:0.00 Class :character
## Mode :character Mode :character Median :0.00 Mode :character
## Mean :0.46
## 3rd Qu.:1.00
## Max. :1.00
## CEA TNM CCI adc_decrease
## Min. : 0.630 Min. :1.00 Min. :0 Length:50
## 1st Qu.: 1.450 1st Qu.:3.00 1st Qu.:1 Class :character
## Median : 2.500 Median :4.00 Median :2 Mode :character
## Mean : 5.574 Mean :3.56 Mean :2
## 3rd Qu.: 4.150 3rd Qu.:4.00 3rd Qu.:3
## Max. :115.900 Max. :5.00 Max. :6
## adc_5_decrease
## Min. :-6.717
## 1st Qu.: 2.566
## Median :11.329
## Mean :10.998
## 3rd Qu.:15.638
## Max. :63.617이름이 너무 긴 변수가 있으면 코드가 길어져서 불편합니다. 먼저 이름을 짧게 고쳐줍니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.)Environment 창에서 dat1를 눌러 이름이 잘 바뀌었나 확인해봅시다.
다음으로, character type으로 들어가있는 날짜 변수들을 Date type으로 고쳐줍니다. Date type에 대해서는 이 웹페이지에 잘 정리되어있습니다: https://www.statmethods.net/input/dates.html
먼저, Environment 창에서 dat1를 눌러 OP_date를 직접 확인하거나, 콘솔에서 call을 해서 날짜가 어떤 형식으로 들어가 있나 봅시다.
dat1$OP_date## [1] "2014-07-07" "2013-05-07" "2016-05-06" "2015-11-25" "2015-02-09"
## [6] "2015-06-22" "2017-02-10" "2012-09-20" "2016-05-20" "2017-05-19"
## [11] "2016-01-11" "2015-11-24" "2012-08-14" "2016-07-13" "2015-08-31"
## [16] "2015-12-16" "2015-02-04" "2015-07-23" "2015-10-23" "2011-12-29"
## [21] "2016-10-18" "2012-01-27" "2017-11-14" "2017-07-27" "2015-05-22"
## [26] "2015-02-26" "2017-05-17" "2016-03-14" "2017-03-14" "2016-09-20"
## [31] "2014-08-08" "2016-10-12" "2016-07-22" "2016-01-26" "2016-05-23"
## [36] "2017-06-29" "2016-09-28" "2017-09-05" "2016-02-29" "2013-12-20"
## [41] "2012-03-13" "2015-09-10" "2014-07-31" "2015-01-27" "2014-05-23"
## [46] "2017-08-21" "2014-05-09" "2015-07-21" "2015-02-03" "2015-09-07"그냥 character type이기 때문에, 덧뺄셈이 안됩니다.
dat1$OP_date + 1## Error in dat1$OP_date + 1: non-numeric argument to binary operator따라서 이것을 date type으로 바꾸어주여야 합니다. 지금 dat3$OP_date이 입력된 형식은 “네자리수 연도-두자리수 월 - 두자리수 일” 형식입니다. 이런 형식으로 들어가있으니 이것을 문자열이 아닌 날짜로 인식해라, 라는 명령어는 다음과 같습니다.
as.Date(dat1$OP_date, format="%Y-%m-%d")## [1] "2014-07-07" "2013-05-07" "2016-05-06" "2015-11-25" "2015-02-09"
## [6] "2015-06-22" "2017-02-10" "2012-09-20" "2016-05-20" "2017-05-19"
## [11] "2016-01-11" "2015-11-24" "2012-08-14" "2016-07-13" "2015-08-31"
## [16] "2015-12-16" "2015-02-04" "2015-07-23" "2015-10-23" "2011-12-29"
## [21] "2016-10-18" "2012-01-27" "2017-11-14" "2017-07-27" "2015-05-22"
## [26] "2015-02-26" "2017-05-17" "2016-03-14" "2017-03-14" "2016-09-20"
## [31] "2014-08-08" "2016-10-12" "2016-07-22" "2016-01-26" "2016-05-23"
## [36] "2017-06-29" "2016-09-28" "2017-09-05" "2016-02-29" "2013-12-20"
## [41] "2012-03-13" "2015-09-10" "2014-07-31" "2015-01-27" "2014-05-23"
## [46] "2017-08-21" "2014-05-09" "2015-07-21" "2015-02-03" "2015-09-07"이렇게하면 날짜에 덧뺄셈 등이 가능해집니다.
as.Date(dat1$OP_date, format="%Y-%m-%d") + 1## [1] "2014-07-08" "2013-05-08" "2016-05-07" "2015-11-26" "2015-02-10"
## [6] "2015-06-23" "2017-02-11" "2012-09-21" "2016-05-21" "2017-05-20"
## [11] "2016-01-12" "2015-11-25" "2012-08-15" "2016-07-14" "2015-09-01"
## [16] "2015-12-17" "2015-02-05" "2015-07-24" "2015-10-24" "2011-12-30"
## [21] "2016-10-19" "2012-01-28" "2017-11-15" "2017-07-28" "2015-05-23"
## [26] "2015-02-27" "2017-05-18" "2016-03-15" "2017-03-15" "2016-09-21"
## [31] "2014-08-09" "2016-10-13" "2016-07-23" "2016-01-27" "2016-05-24"
## [36] "2017-06-30" "2016-09-29" "2017-09-06" "2016-03-01" "2013-12-21"
## [41] "2012-03-14" "2015-09-11" "2014-08-01" "2015-01-28" "2014-05-24"
## [46] "2017-08-22" "2014-05-10" "2015-07-22" "2015-02-04" "2015-09-08"이렇게 date type으로 바꾼 벡터를 dat1$OP_date에 덮어쓰면 되겠죠. 그걸 “옛날” 방식으로 하려면 아래와 같이 하면 됩니다.
dat1$OP_date<-as.Date(dat1$OP_date, format="%Y-%m-%d")그렇지만 이걸 “요즘” 방식으로 dplyr를 써서 하려면 mutate()를 써서 다음과 같이 하면 됩니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate(OP_date = as.Date(OP_date, format="%Y-%m-%d"))마찬가지로 character type으로 저장된 날짜 변수 Recur_date에 대해서도 똑같이 해주면 아래와 같이 되겠죠.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate(OP_date = as.Date(OP_date, format="%Y-%m-%d"),
Recur_date = as.Date(Recur_date, format="%Y-%m-%d"))그런데 이렇게 바꿔줘야할 날짜 변수가 한 열개 쯤 있으면 어떻게 하는 게 효율적일까요? 옛날 방식으로 10줄의 코드를 실행시켜도 되긴 합니다. 하지만 dplyr의 mutate_at을 이용하면 훨씬 짧은 코드로 같은 일을 해낼 수 있습니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d")))\(~\)
이제, 연속변수여야 하는데 character type으로 저장된 CA19.9, adc_decrease를 확인해봅시다. 숫자가 아닌 문자가 섞여있어서 그렇다는 것을 알수 있죠. 이런 경우가 정말정말 많기 때문에 꼭 확인해야합니다. 그 문자들을 replace() 함수를 써서 숫자나 NA로 바꾸어줍니다. (R에서 NA는 결측을 의미합니다.) 그렇게 값을 바꾸어 준 후에도 여전히 character type으로 남아있기 때문에, as.numeric() 또는 as.double()을 써서 numeric type으로 바꾸어주는 것을 잊으면 안됩니다.
일단 “옛날” 방식으로 천천히 다음 코드를 이해해봅시다.
dat1$CA19.9## [1] "4490" "37.8" "36.6" "65" "38.9" "5.9" "23.3" "53.2" "19.8"
## [10] "84.7" "4860" "50.2" "15.9" "18.9" "<1.0" "59.4" "12.6" "81.4"
## [19] "2500" "6.6" "71.1" "4.3" "68.7" "11.4" "10.8" "10.8" "5.9"
## [28] "3.3" "14.4" "47.9" "25.6" "134.1" "10.7" "472" "160.6" "23.2"
## [37] "20.4" "3.4" "55.3" "9.2" "5.1" "12.1" "2647" "24.5" "17.9"
## [46] "33.8" "5.8" "5.6" "<1.0" "99.3"replace(dat1$CA19.9, dat1$CA19.9=="<1.0", 1)## [1] "4490" "37.8" "36.6" "65" "38.9" "5.9" "23.3" "53.2" "19.8"
## [10] "84.7" "4860" "50.2" "15.9" "18.9" "1" "59.4" "12.6" "81.4"
## [19] "2500" "6.6" "71.1" "4.3" "68.7" "11.4" "10.8" "10.8" "5.9"
## [28] "3.3" "14.4" "47.9" "25.6" "134.1" "10.7" "472" "160.6" "23.2"
## [37] "20.4" "3.4" "55.3" "9.2" "5.1" "12.1" "2647" "24.5" "17.9"
## [46] "33.8" "5.8" "5.6" "1" "99.3"as.double(replace(dat1$CA19.9, dat1$CA19.9=="<1.0", 1))## [1] 4490.0 37.8 36.6 65.0 38.9 5.9 23.3 53.2 19.8 84.7
## [11] 4860.0 50.2 15.9 18.9 1.0 59.4 12.6 81.4 2500.0 6.6
## [21] 71.1 4.3 68.7 11.4 10.8 10.8 5.9 3.3 14.4 47.9
## [31] 25.6 134.1 10.7 472.0 160.6 23.2 20.4 3.4 55.3 9.2
## [41] 5.1 12.1 2647.0 24.5 17.9 33.8 5.8 5.6 1.0 99.3위의 코드가 이해되시면, 아래의 코드도 이해가 되시겠죠?
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA))) 이제, 범주형 변수여야 하는데 연속변수로 들어가있는 경우를 찾아서 as.factor()를 이용하여 factor type으로 고쳐줍니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA))) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI),
list(~as.factor(.)))모든 type이 올바르게 들어가있는지 확인합니다.
summary(dat1)## id age sex Recur OP_date
## Min. : 1.00 Min. :25.00 0:17 0:28 Min. :2011-12-29
## 1st Qu.:13.25 1st Qu.:49.25 1:33 1:22 1st Qu.:2015-01-28
## Median :25.50 Median :58.50 Median :2015-11-08
## Mean :25.50 Mean :56.44 Mean :2015-08-11
## 3rd Qu.:37.75 3rd Qu.:65.75 3rd Qu.:2016-09-05
## Max. :50.00 Max. :79.00 Max. :2017-11-14
##
## Recur_date local_6m CA19.9 CEA TNM
## Min. :2012-04-09 0:27 Min. : 1.00 Min. : 0.630 1: 2
## 1st Qu.:2016-06-25 1:23 1st Qu.: 10.72 1st Qu.: 1.450 2: 5
## Median :2018-03-13 Median : 23.25 Median : 2.500 3:12
## Mean :2017-09-21 Mean : 329.53 Mean : 5.574 4:25
## 3rd Qu.:2019-10-12 3rd Qu.: 63.60 3rd Qu.: 4.150 5: 6
## Max. :2020-01-11 Max. :4860.00 Max. :115.900
##
## CCI adc_decrease adc_5_decrease
## 0: 3 Min. :-2.474 Min. :-6.717
## 1:16 1st Qu.: 8.507 1st Qu.: 2.566
## 2:15 Median :17.409 Median :11.329
## 3:12 Mean :21.558 Mean :10.998
## 4: 3 3rd Qu.:37.885 3rd Qu.:15.638
## 6: 1 Max. :70.525 Max. :63.617
## NA's :2\(~\)
Creating new variables necessary for analysis
이제 모든 변수가 올바른 type으로 들어가 있지만, 아직 분석을 할 준비가 덜 되었습니다. 연속변수의 경우 분포를 확인하고, 변환이 필요한 경우 변환을 취해줍니다.
hist(dat1$CA19.9)
hist(log(dat1$CA19.9))
hist(dat1$CEA)
hist(log(dat1$CEA))
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI),
list(~as.factor(.)))연속변수의 경우 cutoff를 적용해서 범주형 변수를 만들어 사용하기도 하죠. 이 예제에서는 CA 19-9이 37초과면 비정상, 37이하면 정상, 이렇게 두 그룹으로 나누어보겠습니다. 이런 것은 ifelse()함수를 사용하면 간단합니다. 단, factor type으로 바꾸는 것을 잊으면 안됩니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
CA19.9.group=ifelse(CA19.9>37, 1, 0)
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group),
list(~as.factor(.)))범주형 변수는 여러개의 카테고리를 하나로 합쳐야 할 경우가 있습니다. 그럴 때는 mapvalues() 함수를 사용하는 것이 편합니다. 마찬가지로, factor type으로 바꾸는 것을 잊으면 안됩니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
CA19.9.group=ifelse(CA19.9>37, 1, 0),
TNM.b=mapvalues(TNM, from=c(1, 2, 3, 4, 5), to=c("1or2", "1or2", "3", "4", "5")),
CCI.b=mapvalues(CCI, from=c(0, 1, 2, 3, 4, 6),
to=c("0or1", "0or1", "2", "3", ">=4", ">=4")),
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b, CCI.b),
list(~as.factor(.)))마지막으로, 생존분석을 할 경우 이벤트 날짜에서 time zero를 뺀, time to event를 계산해야 합니다. 아까 날짜를 모두 date type으로 바꾸어주었기 때문에 쉽게 계산할 수 있습니다. 단, 주의할 점은, 두 날짜의 차이를 numeric type으로 다시 바꾸어주어야 한다는 점입니다. as.numeric()이나 as.double()을 사용하면 됩니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
CA19.9.group=ifelse(CA19.9>37, 1, 0),
TNM.b=mapvalues(TNM, from=c(1, 2, 3, 4, 5), to=c("1or2", "1or2", "3", "4", "5")),
CCI.b=mapvalues(CCI, from=c(0, 1, 2, 3, 4, 6),
to=c("0or1", "0or1", "2", "3", ">=4", ">=4")),
RFS=as.double(Recur_date-OP_date)
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b, CCI.b),
list(~as.factor(.)))마지막으로, 분석에 필요한 모든 변수가 준비되었는지 확인합니다.
summary(dat1)## id age sex Recur OP_date
## Min. : 1.00 Min. :25.00 0:17 0:28 Min. :2011-12-29
## 1st Qu.:13.25 1st Qu.:49.25 1:33 1:22 1st Qu.:2015-01-28
## Median :25.50 Median :58.50 Median :2015-11-08
## Mean :25.50 Mean :56.44 Mean :2015-08-11
## 3rd Qu.:37.75 3rd Qu.:65.75 3rd Qu.:2016-09-05
## Max. :50.00 Max. :79.00 Max. :2017-11-14
##
## Recur_date local_6m CA19.9 CEA TNM
## Min. :2012-04-09 0:27 Min. : 1.00 Min. : 0.630 1: 2
## 1st Qu.:2016-06-25 1:23 1st Qu.: 10.72 1st Qu.: 1.450 2: 5
## Median :2018-03-13 Median : 23.25 Median : 2.500 3:12
## Mean :2017-09-21 Mean : 329.53 Mean : 5.574 4:25
## 3rd Qu.:2019-10-12 3rd Qu.: 63.60 3rd Qu.: 4.150 5: 6
## Max. :2020-01-11 Max. :4860.00 Max. :115.900
##
## CCI adc_decrease adc_5_decrease log.CA19.9 log.CEA
## 0: 3 Min. :-2.474 Min. :-6.717 Min. :0.000 Min. :-0.4620
## 1:16 1st Qu.: 8.507 1st Qu.: 2.566 1st Qu.:2.373 1st Qu.: 0.3699
## 2:15 Median :17.409 Median :11.329 Median :3.146 Median : 0.9163
## 3:12 Mean :21.558 Mean :10.998 Mean :3.407 Mean : 0.9885
## 4: 3 3rd Qu.:37.885 3rd Qu.:15.638 3rd Qu.:4.152 3rd Qu.: 1.4229
## 6: 1 Max. :70.525 Max. :63.617 Max. :8.489 Max. : 4.7527
## NA's :2
## CA19.9.group TNM.b CCI.b RFS
## 0:30 1or2: 7 >=4 : 4 Min. : 5.0
## 1:20 3 :12 0or1:19 1st Qu.: 221.5
## 4 :25 2 :15 Median : 631.0
## 5 : 6 3 :12 Mean : 771.9
## 3rd Qu.:1276.2
## Max. :2533.0
## 한가지 문제되는 점을 찾으셨나요? 그 점을 고치는 코드입니다.
dat1<- dat0 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(
adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
CA19.9.group=ifelse(CA19.9>37, 1, 0),
TNM.b=mapvalues(TNM, from=c(1, 2, 3, 4, 5), to=c("1or2", "1or2", "3", "4", "5")),
CCI.b=mapvalues(CCI, from=c(0, 1, 2, 3, 4, 6),
to=c("0or1", "0or1", "2", "3", ">=4", ">=4")),
CCI.b=factor(CCI.b, levels=c("0or1", "2", "3", ">=4")),
RFS=as.double(Recur_date-OP_date)
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b, CCI.b),
list(~as.factor(.)))분석에 필요한 모든 변수가 준비되었는지 다시한번 확인합니다.
summary(dat1)## id age sex Recur OP_date
## Min. : 1.00 Min. :25.00 0:17 0:28 Min. :2011-12-29
## 1st Qu.:13.25 1st Qu.:49.25 1:33 1:22 1st Qu.:2015-01-28
## Median :25.50 Median :58.50 Median :2015-11-08
## Mean :25.50 Mean :56.44 Mean :2015-08-11
## 3rd Qu.:37.75 3rd Qu.:65.75 3rd Qu.:2016-09-05
## Max. :50.00 Max. :79.00 Max. :2017-11-14
##
## Recur_date local_6m CA19.9 CEA TNM
## Min. :2012-04-09 0:27 Min. : 1.00 Min. : 0.630 1: 2
## 1st Qu.:2016-06-25 1:23 1st Qu.: 10.72 1st Qu.: 1.450 2: 5
## Median :2018-03-13 Median : 23.25 Median : 2.500 3:12
## Mean :2017-09-21 Mean : 329.53 Mean : 5.574 4:25
## 3rd Qu.:2019-10-12 3rd Qu.: 63.60 3rd Qu.: 4.150 5: 6
## Max. :2020-01-11 Max. :4860.00 Max. :115.900
##
## CCI adc_decrease adc_5_decrease log.CA19.9 log.CEA
## 0: 3 Min. :-2.474 Min. :-6.717 Min. :0.000 Min. :-0.4620
## 1:16 1st Qu.: 8.507 1st Qu.: 2.566 1st Qu.:2.373 1st Qu.: 0.3699
## 2:15 Median :17.409 Median :11.329 Median :3.146 Median : 0.9163
## 3:12 Mean :21.558 Mean :10.998 Mean :3.407 Mean : 0.9885
## 4: 3 3rd Qu.:37.885 3rd Qu.:15.638 3rd Qu.:4.152 3rd Qu.: 1.4229
## 6: 1 Max. :70.525 Max. :63.617 Max. :8.489 Max. : 4.7527
## NA's :2
## CA19.9.group TNM.b CCI.b RFS
## 0:30 1or2: 7 0or1:19 Min. : 5.0
## 1:20 3 :12 2 :15 1st Qu.: 221.5
## 4 :25 3 :12 Median : 631.0
## 5 : 6 >=4 : 4 Mean : 771.9
## 3rd Qu.:1276.2
## Max. :2533.0
## \(~\)
Generating summary statistics using tableone package
분석을 시작하기 전에 반드시 요약통계량을 뽑고 이상이 없는지 확인해야합니다. 일단, 분석에 필요한 변수 이름을 저장한 벡터를 만듭니다.
var1<-c("Recur", "RFS", "local_6m", "age", "sex", "CA19.9", "log.CA19.9", "CA19.9.group", "CEA", "log.CEA",
"TNM", "TNM.b", "CCI", "CCI.b", "adc_decrease", "adc_5_decrease")그리고 CreateTableOne() 함수안에 변수 이름과 data frame을 넣고 결과물을 저장한 후, 그 결과물을 summary() 함수에 넣어주면 요약통계량이 정리되어서 나옵니다. 이것을 보고 최소값, 최대값, 결측치 갯수 등을 확인하고 이상이 있으면 고쳐줍니다.
tobj<-CreateTableOne(vars = var1, data=dat1)
summary(tobj)##
## ### Summary of continuous variables ###
##
## strata: Overall
## n miss p.miss mean sd median p25 p75 min max skew kurt
## RFS 50 0 0 772 6e+02 631.0 221.5 1276 5.0 2533 0.8 -0.3
## age 50 0 0 56 1e+01 58.5 49.2 66 25.0 79 -0.5 -0.3
## CA19.9 50 0 0 330 1e+03 23.2 10.7 64 1.0 4860 3.6 12.7
## log.CA19.9 50 0 0 3 2e+00 3.1 2.4 4 0.0 8 1.0 1.5
## CEA 50 0 0 6 2e+01 2.5 1.4 4 0.6 116 6.6 45.5
## log.CEA 50 0 0 1 9e-01 0.9 0.4 1 -0.5 5 1.8 5.8
## adc_decrease 50 2 4 22 2e+01 17.4 8.5 38 -2.5 71 0.5 -0.3
## adc_5_decrease 50 0 0 11 1e+01 11.3 2.6 16 -6.7 64 1.7 6.0
##
## =======================================================================================
##
## ### Summary of categorical variables ###
##
## strata: Overall
## var n miss p.miss level freq percent cum.percent
## Recur 50 0 0.0 0 28 56.0 56.0
## 1 22 44.0 100.0
##
## local_6m 50 0 0.0 0 27 54.0 54.0
## 1 23 46.0 100.0
##
## sex 50 0 0.0 0 17 34.0 34.0
## 1 33 66.0 100.0
##
## CA19.9.group 50 0 0.0 0 30 60.0 60.0
## 1 20 40.0 100.0
##
## TNM 50 0 0.0 1 2 4.0 4.0
## 2 5 10.0 14.0
## 3 12 24.0 38.0
## 4 25 50.0 88.0
## 5 6 12.0 100.0
##
## TNM.b 50 0 0.0 1or2 7 14.0 14.0
## 3 12 24.0 38.0
## 4 25 50.0 88.0
## 5 6 12.0 100.0
##
## CCI 50 0 0.0 0 3 6.0 6.0
## 1 16 32.0 38.0
## 2 15 30.0 68.0
## 3 12 24.0 92.0
## 4 3 6.0 98.0
## 6 1 2.0 100.0
##
## CCI.b 50 0 0.0 0or1 19 38.0 38.0
## 2 15 30.0 68.0
## 3 12 24.0 92.0
## >=4 4 8.0 100.0
## 소숫점 갯수가 제각각이고 scientific notation으로 나온 값도 있어서 보기에 불편하죠. 이런 점을 해결하기 위해 digits=2라는 옵션을 넣어줍니다.
tobj<-CreateTableOne(vars = var1, data=dat1)
summary(tobj, digits=2)##
## ### Summary of continuous variables ###
##
## strata: Overall
## n miss p.miss mean sd median p25 p75 min max skew kurt
## RFS 50 0 0 771.94 637.45 631.00 221.50 1276.2 5.00 2533.0 0.80 -0.29
## age 50 0 0 56.44 12.60 58.50 49.25 65.8 25.00 79.0 -0.47 -0.35
## CA19.9 50 0 0 329.53 1029.52 23.25 10.72 63.6 1.00 4860.0 3.63 12.67
## log.CA19.9 50 0 0 3.41 1.88 3.15 2.37 4.2 0.00 8.5 1.04 1.53
## CEA 50 0 0 5.57 16.28 2.50 1.45 4.2 0.63 115.9 6.63 45.52
## log.CEA 50 0 0 0.99 0.89 0.92 0.37 1.4 -0.46 4.8 1.81 5.82
## adc_decrease 50 2 4 21.56 17.12 17.41 8.51 37.9 -2.47 70.5 0.50 -0.30
## adc_5_decrease 50 0 0 11.00 12.19 11.33 2.57 15.6 -6.72 63.6 1.74 6.04
##
## =======================================================================================
##
## ### Summary of categorical variables ###
##
## strata: Overall
## var n miss p.miss level freq percent cum.percent
## Recur 50 0 0.00 0 28 56.00 56.00
## 1 22 44.00 100.00
##
## local_6m 50 0 0.00 0 27 54.00 54.00
## 1 23 46.00 100.00
##
## sex 50 0 0.00 0 17 34.00 34.00
## 1 33 66.00 100.00
##
## CA19.9.group 50 0 0.00 0 30 60.00 60.00
## 1 20 40.00 100.00
##
## TNM 50 0 0.00 1 2 4.00 4.00
## 2 5 10.00 14.00
## 3 12 24.00 38.00
## 4 25 50.00 88.00
## 5 6 12.00 100.00
##
## TNM.b 50 0 0.00 1or2 7 14.00 14.00
## 3 12 24.00 38.00
## 4 25 50.00 88.00
## 5 6 12.00 100.00
##
## CCI 50 0 0.00 0 3 6.00 6.00
## 1 16 32.00 38.00
## 2 15 30.00 68.00
## 3 12 24.00 92.00
## 4 3 6.00 98.00
## 6 1 2.00 100.00
##
## CCI.b 50 0 0.00 0or1 19 38.00 38.00
## 2 15 30.00 68.00
## 3 12 24.00 92.00
## >=4 4 8.00 100.00
## 이상한 값이 없으면 논문에 넣을 Table 1을 만들어봅시다. Table 1에 들어갈 변수 이름의 벡터를 만들고, 위와 같은 방식으로 CreateTableOne() 함수 안에 넣어 결과물을 저장합니다. 그리고 그 결과물을 ‘print()’ 함수안에 넣으면 Table 1이 나옵니다.
var2<-c("age", "sex", "CA19.9", "CA19.9.group", "CEA", "TNM",
"CCI", "adc_decrease", "adc_5_decrease")
tobj2<-CreateTableOne(vars = var2, data=dat1)
print(tobj2)##
## Overall
## n 50
## age (mean (SD)) 56.44 (12.60)
## sex = 1 (%) 33 (66.0)
## CA19.9 (mean (SD)) 329.53 (1029.52)
## CA19.9.group = 1 (%) 20 (40.0)
## CEA (mean (SD)) 5.57 (16.28)
## TNM (%)
## 1 2 ( 4.0)
## 2 5 (10.0)
## 3 12 (24.0)
## 4 25 (50.0)
## 5 6 (12.0)
## CCI (%)
## 0 3 ( 6.0)
## 1 16 (32.0)
## 2 15 (30.0)
## 3 12 (24.0)
## 4 3 ( 6.0)
## 6 1 ( 2.0)
## adc_decrease (mean (SD)) 21.56 (17.12)
## adc_5_decrease (mean (SD)) 11.00 (12.19)몇가지를 고쳐보겠습니다. 성별과 CA19-9 그룹이 두개의 범주 밖에 없으니 둘 중 하나의 범주에 해당하는 값만 표에 들어갔는데요, 만약, 두 범주 다 나타내고 싶으면 ShowAllLevels = TRUE라는 옵션을 주면 됩니다. 그리고 CA 19-9과 CEA는 분포가 많이 skewed 되어있으니, 평균과 표준편차보다는 median[IQR]로 요약하는 것이 더 좋겠죠. 그렇게 나타내고 싶은 변수의 이름은 nonnormal 옵션으로 지정해주면 됩니다.
print(tobj2, showAllLevels = T, nonnormal = c("CA19.9", "CEA"))##
## level Overall
## n 50
## age (mean (SD)) 56.44 (12.60)
## sex (%) 0 17 (34.0)
## 1 33 (66.0)
## CA19.9 (median [IQR]) 23.25 [10.72, 63.60]
## CA19.9.group (%) 0 30 (60.0)
## 1 20 (40.0)
## CEA (median [IQR]) 2.50 [1.45, 4.15]
## TNM (%) 1 2 ( 4.0)
## 2 5 (10.0)
## 3 12 (24.0)
## 4 25 (50.0)
## 5 6 (12.0)
## CCI (%) 0 3 ( 6.0)
## 1 16 (32.0)
## 2 15 (30.0)
## 3 12 (24.0)
## 4 3 ( 6.0)
## 6 1 ( 2.0)
## adc_decrease (mean (SD)) 21.56 (17.12)
## adc_5_decrease (mean (SD)) 11.00 (12.19)어떤 경우에는 Table 1을 그룹별로 보여주고 싶을 때도 있습니다. 여기에선 CA19-9 그룹에 따라 보여주고 두 그룹을 비교하는 p-value를 계산하기로 하겠습니다. 이런 경우에는 CreateTableOne() 함수안에 strata 옵션으로 그룹을 지정해주면 됩니다. 그리고 그 결과물을 ‘print()’ 함수안에 넣으면 됩니다.
options(width = 300)
tobj3<-CreateTableOne(vars = var2, data=dat1, strata="CA19.9.group")
print(tobj3, nonnormal = c("CA19.9", "CEA"))## Stratified by CA19.9.group
## 0 1 p test
## n 30 20
## age (mean (SD)) 56.50 (12.19) 56.35 (13.50) 0.968
## sex = 1 (%) 16 (53.3) 17 ( 85.0) 0.044
## CA19.9 (median [IQR]) 11.10 [5.82, 19.58] 76.25 [54.77, 238.45] <0.001 nonnorm
## CA19.9.group = 1 (%) 0 ( 0.0) 20 (100.0) <0.001
## CEA (median [IQR]) 1.90 [1.33, 2.88] 3.80 [2.48, 5.23] 0.006 nonnorm
## TNM (%) 0.336
## 1 1 ( 3.3) 1 ( 5.0)
## 2 5 (16.7) 0 ( 0.0)
## 3 8 (26.7) 4 ( 20.0)
## 4 13 (43.3) 12 ( 60.0)
## 5 3 (10.0) 3 ( 15.0)
## CCI (%) 0.828
## 0 2 ( 6.7) 1 ( 5.0)
## 1 10 (33.3) 6 ( 30.0)
## 2 7 (23.3) 8 ( 40.0)
## 3 8 (26.7) 4 ( 20.0)
## 4 2 ( 6.7) 1 ( 5.0)
## 6 1 ( 3.3) 0 ( 0.0)
## adc_decrease (mean (SD)) 22.39 (18.21) 20.40 (15.86) 0.696
## adc_5_decrease (mean (SD)) 12.27 (13.55) 9.08 (9.82) 0.370두가지를 고쳐봅시다. 먼저, 두 그룹의 차이를 p-value 뿐만 아니라 SMD (standardized mean difference)로도 보여주려면, print() 함수 안에 smd=TRUE를 넣어주시면 됩니다. 그리고 TNM과 CCI 두 변수는 각 카테고리별 카운트가 너무 작으니 Chi-squared test 대신에 Fisher’s exact test를 쓰는 것이 더 적합합니다. Fisher’s exact test가 필요한 변수 이름을 exact 옵션으로 정해줍니다.
print(tobj3, nonnormal = c("CA19.9", "CEA"), smd=T, exact=c("TNM", "CCI"))## Stratified by CA19.9.group
## 0 1 p test SMD
## n 30 20
## age (mean (SD)) 56.50 (12.19) 56.35 (13.50) 0.968 0.012
## sex = 1 (%) 16 (53.3) 17 ( 85.0) 0.044 0.730
## CA19.9 (median [IQR]) 11.10 [5.82, 19.58] 76.25 [54.77, 238.45] <0.001 nonnorm 0.731
## CA19.9.group = 1 (%) 0 ( 0.0) 20 (100.0) <0.001 NaN
## CEA (median [IQR]) 1.90 [1.33, 2.88] 3.80 [2.48, 5.23] 0.006 nonnorm 0.397
## TNM (%) 0.308 exact 0.697
## 1 1 ( 3.3) 1 ( 5.0)
## 2 5 (16.7) 0 ( 0.0)
## 3 8 (26.7) 4 ( 20.0)
## 4 13 (43.3) 12 ( 60.0)
## 5 3 (10.0) 3 ( 15.0)
## CCI (%) 0.897 exact 0.442
## 0 2 ( 6.7) 1 ( 5.0)
## 1 10 (33.3) 6 ( 30.0)
## 2 7 (23.3) 8 ( 40.0)
## 3 8 (26.7) 4 ( 20.0)
## 4 2 ( 6.7) 1 ( 5.0)
## 6 1 ( 3.3) 0 ( 0.0)
## adc_decrease (mean (SD)) 22.39 (18.21) 20.40 (15.86) 0.696 0.116
## adc_5_decrease (mean (SD)) 12.27 (13.55) 9.08 (9.82) 0.370 0.270이 정도면 Table 1 이 잘 만들어졌습니다. 통계 분석에서 가장 시간이 오래걸리는 과정이 끝난 셈입니다.