\(~\)

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 5
bb
## [1] 1 2 3 4 5
## Levels: 1 2 3 4 5
cc
## [1] "1" "2" "3" "4" "5"

차이가 보이시나요?

\(~\)

실제 분석에서 type은 중대한 영향을 끼칩니다. 이 벡터들에 1씩을 더해보도록 하겠습니다.

aa+1
## [1] 2 3 4 5 6
bb+1
## Warning in Ops.factor(bb, 1): '+' not meaningful for factors
## [1] NA NA NA NA NA
cc+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과 값들이 이상한 점이 없도록 고쳐주는 것이라고 생각하시면 됩니다.

그럼 본격적으로 데이터 정리를 해봅시다.

\(~\)

Importing the data and checking for duplicates

예제 데이터 example_data.csv를 컴퓨터에 저장한 후 RStudio를 켜기도 전에 해야할 일이 한가지 있습니다. 무엇일까요? (정답은 강의에서 공개하겠습니다!)

RStudio를 켠 후, 각자 데이터가 저장된 위치 경로를 setwd 안에 넣어서 working directory를 셋업합니다. 아래와 같이 역슬래시를 두개씩 붙이거나,

setwd("C:\\Users\\KNOU_stat\\Dropbox\\KNOU_2022_spring\\grad\\ClinicalResearchMethods\\R seminar")

아래와 같이 오른쪽으로 누운 슬래시를 한개씩 붙이셔도 됩니다.

setwd("C:/Users/KNOU_stat/Dropbox/KNOU_2022_spring/grad/ClinicalResearchMethods/R seminar")

이 강의에서는 분석과 자료정리에 dplyr 패키지를 적극 사용하려고 합니다. 처음에는 어려워보이지만 익숙해지면 훨씬 쉽고 빠르게 자료를 정리할 수 있습니다. 이 패키지는 보다 효율적인 데이터분석을 위해 Hadley Wickham이라는 통계학자가 만들고 있는 tidyverse의 일부입니다. 궁금하신 분들은 https://www.tidyverse.org/ 를 참고하세요. 필요한 패키지들을 설치하고 데이터를 read.csv() 함수를 이용해 읽어들입니다.

library(plyr)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.5
library(tableone)

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   volume_change        apt_1_mean       apt_1_95p    
##  Min.   :-6.717   Length:50          Min.   : 0.628   Min.   :1.267  
##  1st Qu.: 2.566   Class :character   1st Qu.: 1.831   1st Qu.:2.517  
##  Median :11.329   Mode  :character   Median : 2.510   Median :3.034  
##  Mean   :10.998                      Mean   : 3.129   Mean   :3.065  
##  3rd Qu.:15.638                      3rd Qu.: 3.100   3rd Qu.:3.650  
##  Max.   :63.617                      Max.   :27.473   Max.   :5.437  
##   adc_1_mean           adc_1_5p     
##  Length:50          Min.   : 53.86  
##  Class :character   1st Qu.: 75.78  
##  Mode  :character   Median : 85.10  
##                     Mean   : 85.03  
##                     3rd Qu.: 92.43  
##                     Max.   :116.59

이름이 너무 긴 변수가 있으면 코드가 길어져서 불편합니다. 먼저 이름을 짧게 고쳐줍니다.

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, volume_change, adc_1_mean을 확인해봅시다. 숫자가 아닌 문자가 섞여있어서 그렇다는 것을 알수 있죠. 이런 경우가 정말정말 많기 때문에 꼭 확인해야합니다. 그 문자들을 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))) %>%
                mutate_at(vars(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", 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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", 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   volume_change        apt_1_mean    
##  0: 3   Min.   :-2.474   Min.   :-6.717   Min.   :-42.2385   Min.   : 0.628  
##  1:16   1st Qu.: 8.507   1st Qu.: 2.566   1st Qu.:  0.8932   1st Qu.: 1.831  
##  2:15   Median :17.409   Median :11.329   Median : 23.9132   Median : 2.510  
##  3:12   Mean   :21.558   Mean   :10.998   Mean   : 33.8975   Mean   : 3.129  
##  4: 3   3rd Qu.:37.885   3rd Qu.:15.638   3rd Qu.: 65.4267   3rd Qu.: 3.100  
##  6: 1   Max.   :70.525   Max.   :63.617   Max.   :140.7366   Max.   :27.473  
##         NA's   :2                         NA's   :2                          
##    apt_1_95p       adc_1_mean        adc_1_5p     
##  Min.   :1.267   Min.   : 82.67   Min.   : 53.86  
##  1st Qu.:2.517   1st Qu.:111.85   1st Qu.: 75.78  
##  Median :3.034   Median :126.66   Median : 85.10  
##  Mean   :3.065   Mean   :124.33   Mean   : 85.03  
##  3rd Qu.:3.650   3rd Qu.:137.30   3rd Qu.: 92.43  
##  Max.   :5.437   Max.   :155.00   Max.   :116.59  
##                  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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", NA))))  %>%
                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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", NA))))  %>%
                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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", NA))))  %>%
                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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", NA))))  %>%
                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   volume_change        apt_1_mean    
##  0: 3   Min.   :-2.474   Min.   :-6.717   Min.   :-42.2385   Min.   : 0.628  
##  1:16   1st Qu.: 8.507   1st Qu.: 2.566   1st Qu.:  0.8932   1st Qu.: 1.831  
##  2:15   Median :17.409   Median :11.329   Median : 23.9132   Median : 2.510  
##  3:12   Mean   :21.558   Mean   :10.998   Mean   : 33.8975   Mean   : 3.129  
##  4: 3   3rd Qu.:37.885   3rd Qu.:15.638   3rd Qu.: 65.4267   3rd Qu.: 3.100  
##  6: 1   Max.   :70.525   Max.   :63.617   Max.   :140.7366   Max.   :27.473  
##         NA's   :2                         NA's   :2                          
##    apt_1_95p       adc_1_mean        adc_1_5p        log.CA19.9   
##  Min.   :1.267   Min.   : 82.67   Min.   : 53.86   Min.   :0.000  
##  1st Qu.:2.517   1st Qu.:111.85   1st Qu.: 75.78   1st Qu.:2.373  
##  Median :3.034   Median :126.66   Median : 85.10   Median :3.146  
##  Mean   :3.065   Mean   :124.33   Mean   : 85.03   Mean   :3.407  
##  3rd Qu.:3.650   3rd Qu.:137.30   3rd Qu.: 92.43   3rd Qu.:4.152  
##  Max.   :5.437   Max.   :155.00   Max.   :116.59   Max.   :8.489  
##                  NA's   :2                                        
##     log.CEA        CA19.9.group  TNM.b     CCI.b         RFS        
##  Min.   :-0.4620   0:30         1or2: 7   >=4 : 4   Min.   :   5.0  
##  1st Qu.: 0.3699   1:20         3   :12   0or1:19   1st Qu.: 221.5  
##  Median : 0.9163                4   :25   2   :15   Median : 631.0  
##  Mean   : 0.9885                5   : 6   3   :12   Mean   : 771.9  
##  3rd Qu.: 1.4229                                    3rd Qu.:1276.2  
##  Max.   : 4.7527                                    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(volume_change, adc_1_mean),
                          list(~as.double(replace(., .=="n/a" | .=="", NA))))  %>%
                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   volume_change        apt_1_mean    
##  0: 3   Min.   :-2.474   Min.   :-6.717   Min.   :-42.2385   Min.   : 0.628  
##  1:16   1st Qu.: 8.507   1st Qu.: 2.566   1st Qu.:  0.8932   1st Qu.: 1.831  
##  2:15   Median :17.409   Median :11.329   Median : 23.9132   Median : 2.510  
##  3:12   Mean   :21.558   Mean   :10.998   Mean   : 33.8975   Mean   : 3.129  
##  4: 3   3rd Qu.:37.885   3rd Qu.:15.638   3rd Qu.: 65.4267   3rd Qu.: 3.100  
##  6: 1   Max.   :70.525   Max.   :63.617   Max.   :140.7366   Max.   :27.473  
##         NA's   :2                         NA's   :2                          
##    apt_1_95p       adc_1_mean        adc_1_5p        log.CA19.9   
##  Min.   :1.267   Min.   : 82.67   Min.   : 53.86   Min.   :0.000  
##  1st Qu.:2.517   1st Qu.:111.85   1st Qu.: 75.78   1st Qu.:2.373  
##  Median :3.034   Median :126.66   Median : 85.10   Median :3.146  
##  Mean   :3.065   Mean   :124.33   Mean   : 85.03   Mean   :3.407  
##  3rd Qu.:3.650   3rd Qu.:137.30   3rd Qu.: 92.43   3rd Qu.:4.152  
##  Max.   :5.437   Max.   :155.00   Max.   :116.59   Max.   :8.489  
##                  NA's   :2                                        
##     log.CEA        CA19.9.group  TNM.b     CCI.b         RFS        
##  Min.   :-0.4620   0:30         1or2: 7   0or1:19   Min.   :   5.0  
##  1st Qu.: 0.3699   1:20         3   :12   2   :15   1st Qu.: 221.5  
##  Median : 0.9163                4   :25   3   :12   Median : 631.0  
##  Mean   : 0.9885                5   : 6   >=4 : 4   Mean   : 771.9  
##  3rd Qu.: 1.4229                                    3rd Qu.:1276.2  
##  Max.   : 4.7527                                    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", "volume_change", 
        "apt_1_mean", "apt_1_95p", "adc_1_mean", "adc_1_5p")

그리고 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.80 -0.3
## age            50    0      0   56 1e+01   58.5  49.2   66  25.0   79 -0.47 -0.3
## CA19.9         50    0      0  330 1e+03   23.2  10.7   64   1.0 4860  3.63 12.7
## log.CA19.9     50    0      0    3 2e+00    3.1   2.4    4   0.0    8  1.04  1.5
## CEA            50    0      0    6 2e+01    2.5   1.4    4   0.6  116  6.63 45.5
## log.CEA        50    0      0    1 9e-01    0.9   0.4    1  -0.5    5  1.81  5.8
## adc_decrease   50    2      4   22 2e+01   17.4   8.5   38  -2.5   71  0.50 -0.3
## adc_5_decrease 50    0      0   11 1e+01   11.3   2.6   16  -6.7   64  1.74  6.0
## volume_change  50    2      4   34 5e+01   23.9   0.9   65 -42.2  141  0.64 -0.4
## apt_1_mean     50    0      0    3 4e+00    2.5   1.8    3   0.6   27  6.06 40.2
## apt_1_95p      50    0      0    3 9e-01    3.0   2.5    4   1.3    5  0.40  0.5
## adc_1_mean     50    2      4  124 2e+01  126.7 111.9  137  82.7  155 -0.24 -0.4
## adc_1_5p       50    0      0   85 1e+01   85.1  75.8   92  53.9  117  0.05  0.5
## 
## =======================================================================================
## 
##      ### 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.803 -0.29
## age            50    0      0  56.44   12.60  58.50  49.25   65.8  25.00   79.0 -0.470 -0.35
## CA19.9         50    0      0 329.53 1029.52  23.25  10.72   63.6   1.00 4860.0  3.633 12.67
## log.CA19.9     50    0      0   3.41    1.88   3.15   2.37    4.2   0.00    8.5  1.045  1.53
## CEA            50    0      0   5.57   16.28   2.50   1.45    4.2   0.63  115.9  6.629 45.52
## log.CEA        50    0      0   0.99    0.89   0.92   0.37    1.4  -0.46    4.8  1.812  5.82
## adc_decrease   50    2      4  21.56   17.12  17.41   8.51   37.9  -2.47   70.5  0.502 -0.30
## adc_5_decrease 50    0      0  11.00   12.19  11.33   2.57   15.6  -6.72   63.6  1.737  6.04
## volume_change  50    2      4  33.90   45.45  23.91   0.89   65.4 -42.24  140.7  0.639 -0.40
## apt_1_mean     50    0      0   3.13    3.70   2.51   1.83    3.1   0.63   27.5  6.055 40.22
## apt_1_95p      50    0      0   3.07    0.88   3.03   2.52    3.6   1.27    5.4  0.404  0.51
## adc_1_mean     50    2      4 124.33   17.11 126.66 111.85  137.3  82.67  155.0 -0.236 -0.41
## adc_1_5p       50    0      0  85.03   11.48  85.10  75.78   92.4  53.86  116.6  0.049  0.52
## 
## =======================================================================================
## 
##      ### 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", "volume_change", 
        "apt_1_mean", "apt_1_95p", "adc_1_mean", "adc_1_5p")

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)  
##   volume_change (mean (SD))   33.90 (45.45)  
##   apt_1_mean (mean (SD))       3.13 (3.70)   
##   apt_1_95p (mean (SD))        3.07 (0.88)   
##   adc_1_mean (mean (SD))     124.33 (17.11)  
##   adc_1_5p (mean (SD))        85.03 (11.48)

몇가지를 고쳐보겠습니다. 성별과 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)       
##   volume_change (mean (SD))         33.90 (45.45)       
##   apt_1_mean (mean (SD))             3.13 (3.70)        
##   apt_1_95p (mean (SD))              3.07 (0.88)        
##   adc_1_mean (mean (SD))           124.33 (17.11)       
##   adc_1_5p (mean (SD))              85.03 (11.48)

어떤 경우에는 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        
##   volume_change (mean (SD))   42.39 (38.63)        20.94 (52.71)          0.111        
##   apt_1_mean (mean (SD))       3.29 (4.67)          2.89 (1.34)           0.713        
##   apt_1_95p (mean (SD))        3.02 (0.76)          3.13 (1.05)           0.670        
##   adc_1_mean (mean (SD))     127.61 (15.08)       119.74 (19.05)          0.117        
##   adc_1_5p (mean (SD))        87.01 (10.88)        82.05 (11.98)          0.136

두가지를 고쳐봅시다. 먼저, 두 그룹의 차이를 p-value 뿐만 아니라 SMD (standardized mean difference)로도 보여주려면, print() 함수 안에 smd=TRUE를 넣어주시면 됩니다. 그리고 TNMCCI 두 변수는 각 카테고리별 카운트가 너무 작으니 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
##   volume_change (mean (SD))   42.39 (38.63)        20.94 (52.71)          0.111          0.464
##   apt_1_mean (mean (SD))       3.29 (4.67)          2.89 (1.34)           0.713          0.116
##   apt_1_95p (mean (SD))        3.02 (0.76)          3.13 (1.05)           0.670          0.120
##   adc_1_mean (mean (SD))     127.61 (15.08)       119.74 (19.05)          0.117          0.458
##   adc_1_5p (mean (SD))        87.01 (10.88)        82.05 (11.98)          0.136          0.433

이 정도면 Table 1 이 잘 만들어졌습니다. 통계 분석에서 가장 시간이 오래걸리는 과정이 끝난 셈입니다.