Data prep and creating table 1
Seo Young Park, PhD
2024-01-15
\(~\)
R을 이용해서 실제로 분석에서 데이터를 읽어들이고 분석 가능한 형태로 만드는 과정에 대해 살펴보겠습니다.
그런데, 실제 데이터 분석에서 가장 문제를 많이 일으키고 우리를 당황시키는 것이 R의 data type 입니다. 아주 간단하게 꼭 필요한 개념만 잠시 살펴보고 넘어가겠습니다.
\(~\)
R basic and R data type
R은 그냥 계산기처럼 사용할 수도 있지만 벡터나 데이터 프레임을 만들어 사용하면 많은 데이터를 쉽게 다룰 수 있습니다. 먼저 값이 1, 2, 3, 4, 5 인 벡터를 만들어보겠습니다.
aa<-c(1, 2, 3, 4, 5)aa는 다섯 개의 숫자가 저장된 벡터입니다. 대괄호를 이용하여 벡터의 값을 불러올 수 있습니다. 몇번째 원소를 불러오고 싶은지 직접 입력해도 되고, 조건문으로 입력해도 됩니다.
aa[1]## [1] 1aa[3]## [1] 3aa[aa>3]## [1] 4 5그런데 1, 2, 3, 4, 5라는 값을 가감승제가 가능한 숫자가 아니라, 범주(factor)나 문자(character)의 형식으로 저장할 수도 있습니다. 아래의 명령문을 실행하면 bb는 factor로, cc는 character로 저장됩니다.
bb<-as.factor(aa)
cc<-as.character(aa)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으로 묶어줍시다.
dd<-data.frame(aa, bb, cc)
dd## aa bb cc
## 1 1 1 1
## 2 2 2 2
## 3 3 3 3
## 4 4 4 4
## 5 5 5 5벡터처럼, 데이터 프레임의 원소도 골라내서 불러올 수 있습니다. 대괄호 안에 [행, 열] 형식으로 불러내고 싶은 원소의 위치를 입력하면 됩니다.
dd[1, ]## aa bb cc
## 1 1 1 1dd[2, ]## aa bb cc
## 2 2 2 2dd[, 2]## [1] 1 2 3 4 5
## Levels: 1 2 3 4 5dd[, 3]## [1] "1" "2" "3" "4" "5"dd[1, 3]## [1] "1"dd[aa>3, ]## aa bb cc
## 4 4 4 4
## 5 5 5 5summary() 함수를 쓰시면, 한눈에 모든 변수들의 type을
짐작할 수 있습니다.
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과 값들이
이상한 점이 없도록 고쳐주는 것이라고 생각하시면 됩니다.
그럼 본격적으로 데이터 정리를 해봅시다.
\(~\)
Importing the data and checking for duplicates
예제 데이터 example_data.csv를 컴퓨터에 저장한 후
RStudio를 켜기도 전에 해야할 일이 한가지 있습니다. 무엇일까요? (정답은
강의에서 공개하겠습니다!)
RStudio를 켠 후, 각자 데이터가 저장된 위치 경로를 setwd 안에 넣어서 working directory를 셋업합니다. 아래와 같이 역슬래시를 두개씩 붙이거나,
setwd("C:\\Users\\biostat81\\Dropbox\\Teaching\\AMC Radiology R 2024")아래와 같이 오른쪽으로 누운 슬래시를 한개씩 붙이셔도 됩니다.
setwd("C:/Users/biostat81/Dropbox/Teaching/AMC Radiology R 2024")필요한 패키지들을 설치하고 로드한 후 데이터를 read.csv()
함수를 이용해 읽어들입니다. 이번 시간엔 plyr, dplyr, tableone 패키지를
이용합니다. 아직 이 패키지를 설치하지 않으신 분들은 아래 명령어를
실행해서 패키지를 설치해주세요.
install.packages("plyr")
install.packages("dplyr")
install.packages("tableone")R에서 패키지를 쓰려면 ’설치’한 이후에 ’로드’를 해야합니다. RStudio를 껐다가 켜도 이미 설치한 패키지가 제거되진 않기 때문에 다시 설치할 필요는 없지만, 로드는 다시 해줘야 합니다. 아래 명령어를 이용해서 패키지를 로드합니다.
library(plyr)
library(dplyr)
library(tableone)이 강의에서는 분석과 자료정리에 dplyr 패키지를 적극
사용하려고 합니다. 처음에는 어려워보이지만 익숙해지면 훨씬 쉽고 빠르게
자료를 정리할 수 있습니다. 이 패키지는 보다 효율적인 데이터분석을 위해
Hadley Wickham이라는 통계학자가 만들고 있는 tidyverse의 일부입니다.
궁금하신 분들은 https://www.tidyverse.org/ 를 참고하세요.
이제 read.csv() 함수를 이용해 데이터를 읽어들입니다.
dat0<-read.csv("ex_data.csv")데이터를 이렇게 읽어들인 이후 꼭 해야하는 일이 있습니다. 무엇일까요? 그 일을 하고 나면 아래와 같이 코드를 고쳐야 한다는 것을 알게 됩니다.
dat0<-read.csv("ex_data.csv", nrow=52)그리고 나서 꼭 해야하는 중요한 일이 있습니다. 각 행이 unique한지
확인하는 것이죠. 일단 RStudio의 Environment 창에서 dat0 을
클릭해서 확인했을 때 31번 환자의 데이터가 중복된 것을 확인할 수
있습니다. 이렇게 완벽하게 똑같은 행이 여러개 들어간 경우는
distinct() 함수를 이용해서 자동으로 없앨 수 있습니다.
dat1<-dat0 %>% distinct()31번 환자의 중복된 데이터가 지워져서 이제 51행의 데이터 프레임이 되었습니다. 이제 이 51행의 환자 아이디가 unique 한지 확인합니다.
length(unique(dat1$id))## [1] 50유니크한 환자 아이디는 50개 입니다. 즉 50명의 환자가 51개의 행으로 들어갔으니 중복된 아이디가 있다는 뜻입니다. 어떤 아이디가 중복되었는지 찾아내서 고칩니다.
which(duplicated(dat1$id))## [1] 45dat1$id[which(duplicated(dat1$id))]## [1] 41dat1[dat1$id==41, ]## id age sex Recur..1..Recur..0..Censored. OP_date Recur_date local_6m
## 41 41 54 0 0 2012-03-13 2019-02-18 0
## 45 41 54 0 0 2012-03-13 2019-02-18 0
## CA19.9 CEA TNM CCI adc_decrease adc_5_decrease apt_1_mean apt_1_95p
## 41 5.1 0.93 2 3 70.52548981 63.61661 2.23 2.781
## 45 5.1 0.93 2 3 0 0.00000 0.00 2.781
## adc_1_5p
## 41 98.61853
## 45 98.6185341번 환자의 데이터가 두번 들어갔다는 걸 알 수 있죠. 31번 환자와는 달리, 환자 아이디는 같지만 다른 변수 값들이 다릅니다. 이런 경우, 어느 쪽이 틀린 정보인지 찾아내서 지워줍니다.
dat2<-dat1[-45, ]
length(unique(dat2$id))## [1] 50nrow(dat2)## [1] 50이제 유니크한 아이디 개수와 행의 수가 같습니다.
\(~\)
Correcting the types of the variables
자, 이제 summary()에 dat2를 넣어서 이상한
점이 있나 확인해봅시다.
summary(dat2)## 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 apt_1_mean apt_1_95p adc_1_5p
## Min. :-6.717 Min. : 0.628 Min. :1.267 Min. : 53.86
## 1st Qu.: 2.566 1st Qu.: 1.831 1st Qu.:2.517 1st Qu.: 75.78
## Median :11.329 Median : 2.510 Median :3.034 Median : 85.10
## Mean :10.998 Mean : 3.129 Mean :3.065 Mean : 85.03
## 3rd Qu.:15.638 3rd Qu.: 3.100 3rd Qu.:3.650 3rd Qu.: 92.43
## Max. :63.617 Max. :27.473 Max. :5.437 Max. :116.59이름이 너무 긴 변수가 있으면 코드가 길어져서 불편합니다. 먼저 이름을 짧게 고쳐줍니다.
dat3<- dat2 %>% rename(Recur=Recur..1..Recur..0..Censored.)Environment 창에서 dat3를 눌러 이름이 잘 바뀌었나 확인해봅시다.
다음으로, character type으로 들어가있는 날짜 변수들을 Date type으로 고쳐줍니다. Date type에 대해서는 이 웹페이지에 잘 정리되어있습니다: https://www.statmethods.net/input/dates.html
먼저, Environment 창에서 dat3를 눌러 OP_date를 직접 확인하거나, 콘솔에서 call을 해서 날짜가 어떤 형식으로 들어가 있나 봅시다.
dat3$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이기 때문에, 덧뺄셈이 안됩니다.
dat3$OP_date + 1## Error in dat3$OP_date + 1: non-numeric argument to binary operator따라서 이것을 date type으로 바꾸어주여야 합니다. 지금
dat3$OP_date이 입력된 형식은 “네자리수 연도-두자리수 월 -
두자리수 일” 형식입니다. 이런 형식으로 들어가있으니 이것을 문자열이 아닌
날짜로 인식해라, 라는 명령어는 다음과 같습니다.
as.Date(dat3$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(dat3$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으로 바꾼 벡터를 dat3$OP_date에 덮어쓰면
되겠죠. 그걸 “옛날” 방식으로 하려면 아래와 같이 하면 됩니다.
dat3$OP_date<-as.Date(dat3$OP_date, format="%Y-%m-%d")그렇지만 이걸 “요즘” 방식으로 dplyr를 써서 하려면
mutate()를 써서 다음과 같이 하면 됩니다.
dat3<- dat2 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate(OP_date = as.Date(OP_date, format="%Y-%m-%d"))마찬가지로 character type으로 저장된 날짜 변수 Recur_date에 대해서도 똑같이 해주면 아래와 같이 되겠죠.
dat3<- dat2 %>% 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을 이용하면 훨씬 짧은 코드로 같은 일을
해낼 수 있습니다.
dat3<- dat2 %>% 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으로 바꾸어주는 것을 잊으면
안됩니다.
일단 “옛날” 방식으로 천천히 다음 코드를 이해해봅시다.
dat2$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(dat2$CA19.9, dat2$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(dat2$CA19.9, dat2$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위의 코드가 이해되시면, 아래의 코드도 이해가 되시겠죠?
dat3<-dat2 %>% 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도 na와 .으로 입력된 값들을 NA로 바꾸어줍니다.
dat3<-dat2 %>% 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으로 고쳐줍니다.
dat3<-dat2 %>% 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(dat3)## 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 apt_1_mean apt_1_95p
## 0: 3 Min. :-2.474 Min. :-6.717 Min. : 0.628 Min. :1.267
## 1:16 1st Qu.: 8.507 1st Qu.: 2.566 1st Qu.: 1.831 1st Qu.:2.517
## 2:15 Median :17.409 Median :11.329 Median : 2.510 Median :3.034
## 3:12 Mean :21.558 Mean :10.998 Mean : 3.129 Mean :3.065
## 4: 3 3rd Qu.:37.885 3rd Qu.:15.638 3rd Qu.: 3.100 3rd Qu.:3.650
## 6: 1 Max. :70.525 Max. :63.617 Max. :27.473 Max. :5.437
## NA's :2
## adc_1_5p
## Min. : 53.86
## 1st Qu.: 75.78
## Median : 85.10
## Mean : 85.03
## 3rd Qu.: 92.43
## Max. :116.59
## \(~\)
Creating new variables necessary for analysis
이제 모든 변수가 올바른 type으로 들어가 있지만, 아직 분석을 할 준비가 덜 되었습니다. 연속변수의 경우 분포를 확인하고, 변환이 필요한 경우 변환을 취해줍니다.
hist(dat3$CA19.9)
hist(log(dat3$CA19.9))
hist(dat3$CEA)
hist(log(dat3$CEA))
CA19.9과 CEA에 각각 로그변환을 취한 변수 log.CA19.9와 log.CEA를 만들어줍니다.
dat3<-dat2 %>% 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으로 바꾸는 것을 잊으면
안됩니다.
dat3<-dat2 %>% 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으로 바꾸는 것을 잊으면 안됩니다.
dat3<-dat2 %>% 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"))
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b),
list(~as.factor(.)))마지막으로, 생존분석을 할 경우 이벤트 날짜에서 time zero를 뺀, time
to event를 계산해야 합니다. 아까 날짜를 모두 date type으로 바꾸어주었기
때문에 쉽게 계산할 수 있습니다. 단, 주의할 점은, 두 날짜의 차이를
numeric type으로 다시 바꾸어주어야 한다는 점입니다.
as.numeric()이나 as.double()을 사용하면
됩니다.
dat3<-dat2 %>% 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")),
RFS=as.double(Recur_date - OP_date)
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b),
list(~as.factor(.)))마지막으로, 분석에 필요한 모든 변수가 준비되었는지 확인합니다.
summary(dat3)## 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 apt_1_mean apt_1_95p
## 0: 3 Min. :-2.474 Min. :-6.717 Min. : 0.628 Min. :1.267
## 1:16 1st Qu.: 8.507 1st Qu.: 2.566 1st Qu.: 1.831 1st Qu.:2.517
## 2:15 Median :17.409 Median :11.329 Median : 2.510 Median :3.034
## 3:12 Mean :21.558 Mean :10.998 Mean : 3.129 Mean :3.065
## 4: 3 3rd Qu.:37.885 3rd Qu.:15.638 3rd Qu.: 3.100 3rd Qu.:3.650
## 6: 1 Max. :70.525 Max. :63.617 Max. :27.473 Max. :5.437
## NA's :2
## adc_1_5p log.CA19.9 log.CEA CA19.9.group TNM.b
## Min. : 53.86 Min. :0.000 Min. :-0.4620 0:30 1or2: 7
## 1st Qu.: 75.78 1st Qu.:2.373 1st Qu.: 0.3699 1:20 3 :12
## Median : 85.10 Median :3.146 Median : 0.9163 4 :25
## Mean : 85.03 Mean :3.407 Mean : 0.9885 5 : 6
## 3rd Qu.: 92.43 3rd Qu.:4.152 3rd Qu.: 1.4229
## Max. :116.59 Max. :8.489 Max. : 4.7527
##
## RFS
## Min. : 5.0
## 1st Qu.: 221.5
## Median : 631.0
## 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", "adc_decrease", "adc_5_decrease", "apt_1_mean")그리고 CreateTableOne() 함수안에 변수 이름과 data
frame을 넣고 결과물을 저장한 후, 그 결과물을 summary()
함수에 넣어주면 요약통계량이 정리되어서 나옵니다. 이것을 보고 최소값,
최대값, 결측치 갯수 등을 확인하고 이상이 있으면 고쳐줍니다.
tobj<-CreateTableOne(vars=var1, data=dat3)
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
## apt_1_mean 50 0 0 3 4e+00 2.5 1.8 3 0.6 27 6.1 40.2
##
## =======================================================================================
##
## ### 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
## 소숫점 갯수가 제각각이고 scientific notation으로 나온 값도 있어서
보기에 불편하죠. 이런 점을 해결하기 위해 digits=2라는
옵션을 넣어줍니다. (이 옵션이 일관적으로 잘 먹히지는 않네요.)
summary(tobj, digits=2)##
## ### Summary of continuous variables ###
##
## strata: Overall
## n miss p.miss mean sd median p25 p75 min max
## RFS 50 0 0 771.94 637.45 631.00 221.50 1276.2 5.00 2533.0
## age 50 0 0 56.44 12.60 58.50 49.25 65.8 25.00 79.0
## CA19.9 50 0 0 329.53 1029.52 23.25 10.72 63.6 1.00 4860.0
## log.CA19.9 50 0 0 3.41 1.88 3.15 2.37 4.2 0.00 8.5
## CEA 50 0 0 5.57 16.28 2.50 1.45 4.2 0.63 115.9
## log.CEA 50 0 0 0.99 0.89 0.92 0.37 1.4 -0.46 4.8
## adc_decrease 50 2 4 21.56 17.12 17.41 8.51 37.9 -2.47 70.5
## adc_5_decrease 50 0 0 11.00 12.19 11.33 2.57 15.6 -6.72 63.6
## apt_1_mean 50 0 0 3.13 3.70 2.51 1.83 3.1 0.63 27.5
## skew kurt
## RFS 0.80 -0.29
## age -0.47 -0.35
## CA19.9 3.63 12.67
## log.CA19.9 1.04 1.53
## CEA 6.63 45.52
## log.CEA 1.81 5.82
## adc_decrease 0.50 -0.30
## adc_5_decrease 1.74 6.04
## apt_1_mean 6.06 40.22
##
## =======================================================================================
##
## ### 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
## 이상한 값이 없으면 논문에 넣을 Table 1을 만들어봅시다. 위에서
CreateTableOne() 함수로 만든 결과물을 ‘print()’ 함수안에
넣으면 Table 1이 나옵니다.
print(tobj)##
## Overall
## n 50
## Recur = 1 (%) 22 (44.0)
## RFS (mean (SD)) 771.94 (637.45)
## local_6m = 1 (%) 23 (46.0)
## age (mean (SD)) 56.44 (12.60)
## sex = 1 (%) 33 (66.0)
## CA19.9 (mean (SD)) 329.53 (1029.52)
## log.CA19.9 (mean (SD)) 3.41 (1.88)
## CA19.9.group = 1 (%) 20 (40.0)
## CEA (mean (SD)) 5.57 (16.28)
## log.CEA (mean (SD)) 0.99 (0.89)
## TNM (%)
## 1 2 ( 4.0)
## 2 5 (10.0)
## 3 12 (24.0)
## 4 25 (50.0)
## 5 6 (12.0)
## TNM.b (%)
## 1or2 7 (14.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)
## apt_1_mean (mean (SD)) 3.13 (3.70)몇가지를 고쳐보겠습니다. 성별과 CA19-9 그룹이 두개의 범주 밖에 없으니
둘 중 하나의 범주에 해당하는 값만 표에 들어갔는데요, 만약, 두 범주 다
나타내고 싶으면 ShowAllLevels = TRUE라는 옵션을 주면
됩니다. 그리고 CA 19-9과 CEA는 분포가 많이 skewed 되어있으니, 평균과
표준편차보다는 median[IQR]로 요약하는 것이 더 좋겠죠. 그렇게 나타내고
싶은 변수의 이름은 nonnormal 옵션으로 지정해주면
됩니다.
print(tobj, showAllLevels = T, nonnormal = c("CA19.9", "CEA"))##
## level Overall
## n 50
## Recur (%) 0 28 (56.0)
## 1 22 (44.0)
## RFS (mean (SD)) 771.94 (637.45)
## local_6m (%) 0 27 (54.0)
## 1 23 (46.0)
## 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]
## log.CA19.9 (mean (SD)) 3.41 (1.88)
## CA19.9.group (%) 0 30 (60.0)
## 1 20 (40.0)
## CEA (median [IQR]) 2.50 [1.45, 4.15]
## log.CEA (mean (SD)) 0.99 (0.89)
## TNM (%) 1 2 ( 4.0)
## 2 5 (10.0)
## 3 12 (24.0)
## 4 25 (50.0)
## 5 6 (12.0)
## TNM.b (%) 1or2 7 (14.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)
## apt_1_mean (mean (SD)) 3.13 (3.70)어떤 경우에는 Table 1을 그룹별로 보여주고 싶을 때도 있습니다.
여기에선 CA19-9 그룹에 따라 보여주고 두 그룹을 비교하는 p-value를
계산하기로 하겠습니다. 이런 경우에는 CreateTableOne()
함수안에 strata 옵션으로 그룹을 지정해주면 됩니다. 그리고
그 결과물을 ‘print()’ 함수안에 넣으면 됩니다.
options(width = 300)
tobj<-CreateTableOne(vars=var1, data=dat3, strata="CA19.9.group")
print(tobj, showAllLevels = T, nonnormal = c("CA19.9", "CEA"))## Stratified by CA19.9.group
## level 0 1 p test
## n 30 20
## Recur (%) 0 26 ( 86.7) 2 ( 10.0) <0.001
## 1 4 ( 13.3) 18 ( 90.0)
## RFS (mean (SD)) 852.17 (688.56) 651.60 (546.68) 0.280
## local_6m (%) 0 18 ( 60.0) 9 ( 45.0) 0.451
## 1 12 ( 40.0) 11 ( 55.0)
## age (mean (SD)) 56.50 (12.19) 56.35 (13.50) 0.968
## sex (%) 0 14 ( 46.7) 3 ( 15.0) 0.044
## 1 16 ( 53.3) 17 ( 85.0)
## CA19.9 (median [IQR]) 11.10 [5.82, 19.58] 76.25 [54.77, 238.45] <0.001 nonnorm
## log.CA19.9 (mean (SD)) 2.28 (0.91) 5.10 (1.67) <0.001
## CA19.9.group (%) 0 30 (100.0) 0 ( 0.0) <0.001
## 1 0 ( 0.0) 20 (100.0)
## CEA (median [IQR]) 1.90 [1.33, 2.88] 3.80 [2.48, 5.23] 0.006 nonnorm
## log.CEA (mean (SD)) 0.72 (0.66) 1.39 (1.04) 0.007
## TNM (%) 1 1 ( 3.3) 1 ( 5.0) 0.336
## 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)
## TNM.b (%) 1or2 6 ( 20.0) 1 ( 5.0) 0.381
## 3 8 ( 26.7) 4 ( 20.0)
## 4 13 ( 43.3) 12 ( 60.0)
## 5 3 ( 10.0) 3 ( 15.0)
## CCI (%) 0 2 ( 6.7) 1 ( 5.0) 0.828
## 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
## apt_1_mean (mean (SD)) 3.29 (4.67) 2.89 (1.34) 0.713두가지를 고쳐봅시다. 먼저, 두 그룹의 차이를 p-value 뿐만 아니라 SMD
(standardized mean difference)로도 보여주려면, print() 함수
안에 smd=TRUE를 넣어주시면 됩니다. 그리고
TNM과 CCI 두 변수는 각 카테고리별 카운트가
너무 작으니 Chi-squared test 대신에 Fisher’s exact test를 쓰는 것이 더
적합합니다. Fisher’s exact test가 필요한 변수 이름을 exact
옵션으로 정해줍니다.
tobj<-CreateTableOne(vars=var1, data=dat3, strata="CA19.9.group")
print(tobj, showAllLevels = T, nonnormal = c("CA19.9", "CEA"), smd=T, exact=c("TNM", "CCI"))## Stratified by CA19.9.group
## level 0 1 p test SMD
## n 30 20
## Recur (%) 0 26 ( 86.7) 2 ( 10.0) <0.001 2.391
## 1 4 ( 13.3) 18 ( 90.0)
## RFS (mean (SD)) 852.17 (688.56) 651.60 (546.68) 0.280 0.323
## local_6m (%) 0 18 ( 60.0) 9 ( 45.0) 0.451 0.304
## 1 12 ( 40.0) 11 ( 55.0)
## age (mean (SD)) 56.50 (12.19) 56.35 (13.50) 0.968 0.012
## sex (%) 0 14 ( 46.7) 3 ( 15.0) 0.044 0.730
## 1 16 ( 53.3) 17 ( 85.0)
## CA19.9 (median [IQR]) 11.10 [5.82, 19.58] 76.25 [54.77, 238.45] <0.001 nonnorm 0.731
## log.CA19.9 (mean (SD)) 2.28 (0.91) 5.10 (1.67) <0.001 2.100
## CA19.9.group (%) 0 30 (100.0) 0 ( 0.0) <0.001 NaN
## 1 0 ( 0.0) 20 (100.0)
## CEA (median [IQR]) 1.90 [1.33, 2.88] 3.80 [2.48, 5.23] 0.006 nonnorm 0.397
## log.CEA (mean (SD)) 0.72 (0.66) 1.39 (1.04) 0.007 0.772
## TNM (%) 1 1 ( 3.3) 1 ( 5.0) 0.308 exact 0.697
## 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)
## TNM.b (%) 1or2 6 ( 20.0) 1 ( 5.0) 0.381 0.541
## 3 8 ( 26.7) 4 ( 20.0)
## 4 13 ( 43.3) 12 ( 60.0)
## 5 3 ( 10.0) 3 ( 15.0)
## CCI (%) 0 2 ( 6.7) 1 ( 5.0) 0.897 exact 0.442
## 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
## apt_1_mean (mean (SD)) 3.29 (4.67) 2.89 (1.34) 0.713 0.116이 정도면 Table 1 이 잘 만들어졌습니다. 통계 분석에서 가장 시간이 오래걸리는 과정이 끝난 셈입니다. 다음 시간에는 이렇게 정리된 자료를 가지고 univariate, multivariable logistic regression을 R에서 수행하는 방법에 대해서 알아보도록 하겠습니다.