2.1 표를 만드는 첫번째 방법 (table함수)

  1. one way table
data(mtcars)
result = table(mtcars$cyl) # 요약표
result
## 
##  4  6  8 
## 11  7 14
prop.table(result) # 비율표
## 
##       4       6       8 
## 0.34375 0.21875 0.43750
prop.table(result)*100 # 백분율표
## 
##      4      6      8 
## 34.375 21.875 43.750


  1. two way table
  • 표를 만들기는 쉽습니다. 그러나 보기 좋지는 않습니다. 새로운 범주형 변수를 만들어 am이 0이면 ‘automatic’, 1이면 ’manual’이라는 라벨을 붙입니다.
table(mtcars$cyl, mtcars$am) # 요약표
##    
##      0  1
##   4  3  8
##   6  4  3
##   8 12  2
mtcars$tm = factor(mtcars$am, labels = c('automatic', 'manual'))
table(mtcars$cyl, mtcars$tm)
##    
##     automatic manual
##   4         3      8
##   6         4      3
##   8        12      2
  • factor함수가 사용하기 어렵다면, ifelse를 사용할 수 있습니다. ifelse(조건, 조건이 참일 때의 값, 조건이 거짓일 때의 값)의 함수를 사용할 수 있습니다.
  • 표나 배열에서 합계를 붙이려면 addmargins() 함수를 쓸 수 있습니다.
mtcars$tm = ifelse(mtcars$am==0, 'automatic', 'manual')
result = table(mtcars$cyl, mtcars$tm)
result
##    
##     automatic manual
##   4         3      8
##   6         4      3
##   8        12      2
addmargins(result)
##      
##       automatic manual Sum
##   4           3      8  11
##   6           4      3   7
##   8          12      2  14
##   Sum        19     13  32


2.2 표를 만드는 두번째 방법 (xtabs함수)

  1. table() 함수와 달리 xtabs 함수는 규정식을 사용합니다. xtabs(~도수 가로 + 세로)의 형식을 사용합니다.
result1 = xtabs(~cyl+tm, data=mtcars)
result1
##    tm
## cyl automatic manual
##   4         3      8
##   6         4      3
##   8        12      2
  1. xtabs로 만든표는 변수 이름이 붙어 있기 때문에 그림을 그려도 변수의 이름이 붙어서 나옵니다.
plot(result1)

mosaicplot(result1,color=c("tan1","firebrick2")) # 색깔 지정


2.3 표를 만드는 세번째 방법 (prop.table함수)

  1. 비율로 된 표를 원할때는 prop.table 함수를 사용하면 됩니다.
prop.table(result1)
##    tm
## cyl automatic  manual
##   4   0.09375 0.25000
##   6   0.12500 0.09375
##   8   0.37500 0.06250
prop.table(result1)*100
##    tm
## cyl automatic manual
##   4     9.375 25.000
##   6    12.500  9.375
##   8    37.500  6.250
addmargins(prop.table(result1))*100
##      tm
## cyl   automatic  manual     Sum
##   4       9.375  25.000  34.375
##   6      12.500   9.375  21.875
##   8      37.500   6.250  43.750
##   Sum    59.375  40.625 100.000

2.4 카이제곱검정, 피셔검정

  1. 카이제곱검정은 chisq.test 함수를 사용합니다.
result
##    
##     automatic manual
##   4         3      8
##   6         4      3
##   8        12      2
chisq.test(result)
## Warning in chisq.test(result): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  result
## X-squared = 8.7407, df = 2, p-value = 0.01265
  1. 결과값을 보면 p-value가 0.01로 나옵니다. 단 table에서 도수가 2, 3인 칸이 있어서 카이제곱검정을 부정확할 수 있다는 경고가 나옵니다. 이런 경우에는 피셔검정을 시행하면 됩니다.
fisher.test(result)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  result
## p-value = 0.009105
## alternative hypothesis: two.sided
  1. 피셔검정 결과 p-value 0.009로 유의미하게 나옵니다.

2.5 Cochran-Armitage Trend test

  1. 두 개의 범주형 변수사이의 독립성 검정은 카이제곱검정 또는 피셔검정을 이용합니다. 하지만 한 변수가 순위가 있는 ordinal 변수인 경우 trend를 보려면 Cochran-Armitage Trend test를 시행합니다.

  2. acs 데이터의 고혈압(HBP)과 흡연(smoking)과의 관계를 알아봅시다. 함수 round(x,y)의 경우, x는 반올림하려는 숫자, y는 반올림하려는 소수점 자리수 입니다.

require(moonBook)
data(acs)
result = table(acs$HBP, acs$smoking)
result
##      
##       Ex-smoker Never Smoker
##   No         81    99    176
##   Yes       123   233    145
round(prop.table(result)*100,2)
##      
##       Ex-smoker Never Smoker
##   No       9.45 11.55  20.54
##   Yes     14.35 27.19  16.92


  1. acs 데이터에서 흡연력은 smoking 열에 있는데 “Ex-smoker”, “Never”, “Smoker” 세개 중 하나의 문자열 변수로 되어있습니다. 우리가 원하는 순서는 “Never”, “Ex-smoker”, “Smoker” 순입니다. 따라서 문자열 변수를 범주형 변수로 바꾸면서 순서를 정해줍시다.
acs$smoking=factor(acs$smoking, levels = c('Never','Ex-smoker','Smoker'))
str(acs$smoking)
##  Factor w/ 3 levels "Never","Ex-smoker",..: 3 1 1 1 3 1 2 2 1 3 ...
result=table(acs$HBP, acs$smoking)
result
##      
##       Never Ex-smoker Smoker
##   No     99        81    176
##   Yes   233       123    145


  1. prop.trend.test(x, n, score)를 사용합시다. x는 number of events, n은 number of trials로 되어있다. score는 일단 신경쓰지 맙시다. 흡연에 따른 고혈압 발생의 경향을 보려면 number of events는 고혈압이 있는 사람, 즉 결과에서 두번째 행이 됩니다. number of trials는 전체 대상 환자가 됩니다.
res = prop.trend.test(result[2,], colSums(result))
res
## 
##  Chi-squared Test for Trend in Proportions
## 
## data:  result[2, ] out of colSums(result) ,
##  using scores: 1 2 3
## X-squared = 41.967, df = 1, p-value = 9.282e-11
  • p-value가 유의미하게 나왔습니다. 이것을 그래프로 그리면 아래와 같습니다.
plot(t(result),col=c("deepskyblue","brown2"),main="Hypertension and Smoking",
     ylab="Hypertension",xlab="Smoking")


2.6 데이터의 요약

  1. acs는 857명의 자료인데 이 자료의 총콜레스테롤(TC)의 평균을 확인합니다. 결측치가 있어 평균을 구할 수 없으므로 NA를 제외하고 평균을 구합니다.
mean(acs$TC)
## [1] NA
mean(acs$TC, na.rm = T) # 결측치제거 후 계산
## [1] 185.2002


  1. 진단에 따른 총콜레스테롤의 평균을 알고싶으면 아래와 같이 시행합니다.
mytable(Dx~TC, data=acs)
## 
##           Descriptive Statistics by 'Dx'          
## ___________________________________________________ 
##        NSTEMI       STEMI     Unstable Angina   p  
##       (N=153)      (N=304)        (N=400)    
## --------------------------------------------------- 
##  TC 193.7 ± 53.6 183.2 ± 43.4  183.5 ± 48.3   0.057
## ---------------------------------------------------


  1. tapply(X, Index, FUN..) 자주 사용되는 데이터 요약함수 입니다. X는 보고자하는 값이고 Index는 자료를 나누고자 하는 범주형 변수, FUN은 적용시킬 함수입니다. 즉 전체자료 중 X를 Index 별로 나누어 함수 FUN을 적용시키는 것입니다.
tapply(acs$TC, acs$Dx, mean, na.rm=T)
##          NSTEMI           STEMI Unstable Angina 
##        193.7257        183.2020        183.4801


  1. aggregate(formula, data, FUN, … ,subset,na.action=na.omit) 또한 자주 사용되는 데이터 요약함수 입니다. aggregate 함수는 자체에서 NA가 있을 경우 제거합니다.
aggregate(TC~Dx, data=acs, mean)
##                Dx       TC
## 1          NSTEMI 193.7257
## 2           STEMI 183.2020
## 3 Unstable Angina 183.4801


  1. tapply와 aggregate가 하는 역할은 거의 유사합니다. 하지만 tapply 함수 결과는 벡터이고, aggregate 함수의 결과는 데이터프레임 입니다. 또한 평균과 표준편차를 같이 구하는 함수를 만들면 같이 구할 수 있습니다.
myfun = function(x){
  result=c(mean(x),sd(x))
}
aggregate(TC~Dx, data=acs, myfun)
##                Dx      TC.1      TC.2
## 1          NSTEMI 193.72568  53.61899
## 2           STEMI 183.20204  43.38024
## 3 Unstable Angina 183.48010  48.34923


  1. 사실 위의 경우와 같은 간단한 함수는 굳이 이름을 정하고 여러 줄로 쓸 필요가 없다. 함수의 인수로 함수를 쓰는 경우, 함수의 이름없이 간단히 함수를 호출할 수 있는데 이를 무기명함수라고 합니다.
aggregate(TC~Dx, data=acs, function(x) c(mean(x),sd(x)))
##                Dx      TC.1      TC.2
## 1          NSTEMI 193.72568  53.61899
## 2           STEMI 183.20204  43.38024
## 3 Unstable Angina 183.48010  48.34923


  • 무기명 함수를 쓰게되면 훨씬 간결하고 깔끔해 보입니다. 또한 여러가지 변수의 평균과 표준편차를 함께 볼수도 있습니다.
aggregate(cbind(TC,TG,HDLC,LDLC)~Dx, data=acs, function(x) c(mean(x),sd(x)))
##                Dx      TC.1      TC.2      TG.1      TG.2   HDLC.1   HDLC.2
## 1          NSTEMI 193.72568  53.61899 130.12162  88.54968 38.93919 11.87189
## 2           STEMI 183.42594  43.28396 106.65529  72.29184 38.45734 11.01411
## 3 Unstable Angina 183.35857  48.35118 136.80563 101.17377 37.81330 10.85066
##      LDLC.1    LDLC.2
## 1 126.09459  44.73373
## 2 116.92491  39.37014
## 3 112.87724  40.38533


  1. 하지만 moonBook 패키지의 mytable 함수를 쓰면 더 깔끔하게 볼 수 있습니다.
mytable(sex+Dx~TC+TG, data=acs)
## 
##                         Descriptive Statistics stratified by 'sex' and 'Dx'                        
## ____________________________________________________________________________________________________ 
##                            Male                                             Female                     
##     ------------------------------------------------ ----------------------------------------------- 
##        NSTEMI        STEMI     Unstable Angina   p      NSTEMI       STEMI     Unstable Angina   p  
##        (N=103)      (N=220)        (N=247)             (N=50)       (N=84)        (N=153)         
## ---------------------------------------------------------------------------------------------------- 
##  TC 192.6 ± 54.3  184.1 ± 42.6  178.7 ± 44.6   0.036 196.3 ± 52.7 180.7 ± 45.7  191.1 ± 53.1   0.192
##  TG 138.0 ± 100.2 104.3 ± 65.5  144.3 ± 114.2  0.000 112.5 ± 51.1 112.3 ± 87.2  126.3 ± 76.0   0.316
## ----------------------------------------------------------------------------------------------------