R강의3
문건웅
2014년 7월 22일
그룹간의 비교를 한번에 : compareGroups 패키지의 이용
library(compareGroups)## Loading required package: Hmisc
## Loading required package: grid
## Loading required package: lattice
## Loading required package: survival
## Loading required package: splines
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, units
##
## Loading required package: xtable
##
## Attaching package: 'xtable'
##
## The following objects are masked from 'package:Hmisc':
##
## label, label<-
##
## Loading required package: gdata
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
##
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
##
## Attaching package: 'gdata'
##
## The following object is masked from 'package:Hmisc':
##
## combine
##
## The following object is masked from 'package:stats':
##
## nobs
##
## The following object is masked from 'package:utils':
##
## object.size
##
## Loading required package: SNPassoc
## Loading required package: haplo.stats
## Loading required package: mvtnorm
## Loading required package: parallel패키지의 디자인 : class와 method
compareGroups 패키지는 세개의 함수로 구성되어 있다.
- compareGroups 함수는 compareGroups 클래스의 object를 만든다. 이 object는
- 인쇄할수 있고 (print)
- 요약할수 있고 (summary)
- 그래프를 그릴 수 있고(plot)
갱신할 수 있다.(update)
- createTable 함수는 createTable 클래스의 object를 만든다.이 object는
- 인쇄할수 있고 (print)
요약할수 있다.(summary)
**export2latex, export2csv, export2html은 결과를 각각 CSV, LATEX, 또는 HTML형식으로 내보낸다.
그림 1은 패키지 설계도를 보여준다. 
사용할 데이타
이 패키지가 어떻게 동작하는지 보여주기 위해 PREDIMED 연구 (www.predimed.org)의 참여자로부터 85%의 데이타를 추출했다. PREDIMED연구는 스페인에서 시행한 다기관 연구로써 심혈관 위험이 높으나 연구에등록될 당시에는 심혈관 질환이 없는 사람들을 대상으로 세가지 식이요법 중 한가지를 무작위로 배정하였다. 세가지 식이요법은 extra-virgin의 올리브유가 보강된 지중해식식단(MedDiet+VOO), 혼합 너트류로 보강된 지중해식 식단(MedDiet+Nuts), 또는 대조식단(지방 성분을 줄이도록 조언한 식단)이다.참가자들은 분기마다 개인교육 및 집단교육을 받았으며 속한 실험군 에 따라 exra-virgin올리브유와 혼합넛트, 그리고 자그마한 선물을 제공받았다.연구의 일차종결점은주요심혈관사건(심근경색, 뇌졸중,심혈관원인으로 인한 사망)이었다.
먼저 PREDIMED 데이타를 읽어들인다.
data(predimed)이 데이타의 변수들과 라벨은 다음과 같다. 
분석에 앞서 먼저 생존분석을 하는데 필요한 Surv클래스의 변수를 만들어야 한다. 여기서 만들어진 변수 ** tmain** 은 time-to-death 와 time-to-cardiovascular event로 각각 사용된다.
predimed$tmain <- with(predimed, Surv(toevent, event == "Yes"))
label(predimed$tmain) <- "AMI, stroke, or CV Death"compareGroups
이것은 이패키의 주된 함수이다. 이 함수가 모든 계산을 수행하여 그결과를 object에 저장한다. 나중에** createTable**함수를 이 object에 적용하면 분석결과의 table 이 만들어진다. 예를 들어 predimed데이타에서 group을 반응변수로 다른 변수들을 설명변수로 하여 단변량분석을 수행하려면 다음과 같이 한다.
res=compareGroups(group ~ . - toevent - event, data = predimed)
res##
##
## -------- Summary of results by groups of 'Intervention group'---------
##
##
## var N p.value method
## 1 Sex 6324 <0.001** categorical
## 2 Age 6324 0.003** continuous normal
## 3 Smoking 6324 0.444 categorical
## 4 Body mass index 6324 <0.001** continuous normal
## 5 Waist circumference 6324 0.045** continuous normal
## 6 Waist-to-height ratio 6324 <0.001** continuous normal
## 7 Hypertension 6324 0.249 categorical
## 8 Type-2 diabetes 6324 0.017** categorical
## 9 Dyslipidemia 6324 0.423 categorical
## 10 Family history of premature CHD 6324 0.581 categorical
## 11 Hormone-replacement therapy 5661 0.850 categorical
## 12 MeDiet Adherence score 6324 <0.001** continuous normal
## 13 AMI, stroke, or CV Death 6324 0.011** Surv [Tmax=4.79]
## selection
## 1 ALL
## 2 ALL
## 3 ALL
## 4 ALL
## 5 ALL
## 6 ALL
## 7 ALL
## 8 ALL
## 9 ALL
## 10 ALL
## 11 ALL
## 12 ALL
## 13 ALL
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1createTable(res)##
## --------Summary descriptives table by 'Intervention group'---------
##
## ____________________________________________________________________________________
## Control MedDiet + Nuts MedDiet + VOO p.overall
## N=2042 N=2100 N=2182
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Sex: <0.001
## Male 812 (39.8%) 968 (46.1%) 899 (41.2%)
## Female 1230 (60.2%) 1132 (53.9%) 1283 (58.8%)
## Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
## Smoking: 0.444
## Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
## Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
## Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
## Body mass index 30.3 (3.96) 29.7 (3.77) 29.9 (3.71) <0.001
## Waist circumference 101 (10.8) 100 (10.6) 100 (10.4) 0.045
## Waist-to-height ratio 0.63 (0.07) 0.62 (0.06) 0.63 (0.06) <0.001
## Hypertension: 0.249
## No 331 (16.2%) 362 (17.2%) 396 (18.1%)
## Yes 1711 (83.8%) 1738 (82.8%) 1786 (81.9%)
## Type-2 diabetes: 0.017
## No 1072 (52.5%) 1150 (54.8%) 1100 (50.4%)
## Yes 970 (47.5%) 950 (45.2%) 1082 (49.6%)
## Dyslipidemia: 0.423
## No 563 (27.6%) 561 (26.7%) 622 (28.5%)
## Yes 1479 (72.4%) 1539 (73.3%) 1560 (71.5%)
## Family history of premature CHD: 0.581
## No 1580 (77.4%) 1640 (78.1%) 1675 (76.8%)
## Yes 462 (22.6%) 460 (21.9%) 507 (23.2%)
## Hormone-replacement therapy: 0.850
## No 1811 (98.3%) 1835 (98.4%) 1918 (98.2%)
## Yes 31 (1.68%) 30 (1.61%) 36 (1.84%)
## MeDiet Adherence score 8.44 (1.94) 8.81 (1.90) 8.77 (1.97) <0.001
## AMI, stroke, or CV Death 5.80% 3.58% 3.76% 0.011
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯1 반응변수고르기
“~”의 오른쪽에 “.”가 사용된다면 그 데이타프레임에 있는 모든 변수가 사용된다. 변수들중 * toevent와 event*변수를 분석에서 제외하고자 할때는 다음과 같이 한다.
compareGroups(group ~ . -toevent -event ,data=predimed)몇개의 설명변수(예를 들어 age,sex,그리고 waist)를 골라분석하고 그 결과를 compareGroup 클래스의 변수에 저장하려면 다음과 같이 한다.
res=compareGroups(group ~ age + sex + smoke + waist + hormo ,data=predimed)
res
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.003** continuous normal ALL
2 Sex 6324 <0.001** categorical ALL
3 Smoking 6324 0.444 categorical ALL
4 Waist circumference 6324 0.045** continuous normal ALL
5 Hormone-replacement therapy 5661 0.850 categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 전체 6324명의 데이타에서 나이, 성별, 허리둘레 는 모든 예에서 있지만 호르몬대체요법 중에 는 missing data가 존재한다.(아마도 남자 참여자일 것이다)
식이요법 군에 따라 흠연과 호르몬 대체요법에는 약간의 차이가 존재하지만 통계적으로 유의하지는 않았고(p값이 각각 0.714, 0.859) 나이와 성별, 그리고 허리둘레는명백히 다르다.
나이와 허리둘레는 정규분포를 하는 연속변수로 사용되었으며 성별, 흡연, 호르몬대체요법은 범주형 변수로 사용되었다. 필터가 적용되지 않았으므로 (예를 들어 치료받는 환자만 ) selection열에는 ALL로 표시되어 있다.
2. Subsetting
참여자 중 일부의 subset에 국한하여 분석을 하려면 다음과 같이 한다. (예를 들어 여자환자만)
compareGroups(group ~ age + sex + smoke + waist + hormo ,data=predimed,
subset=sex=="Female")Warning: Some levels of 'sex' are removed since no observation in that/those levels
Warning: the standard deviation is zero
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method
1 Age 3645 0.056* continuous normal
2 Sex 3645 . categorical
3 Smoking 3645 0.907 categorical
4 Waist circumference 3645 0.016** continuous normal
5 Hormone-replacement therapy 3459 0.898 categorical
selection
1 sex == "Female"
2 sex == "Female"
3 sex == "Female"
4 sex == "Female"
5 sex == "Female"
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 여자 대상자에 대한 결과만 나온 것을 알 수 있다.
Subset의 조합도 가능하다.예를 들면
compareGroups(group ~ age + smoke + waist + hormo, data=predimed,
selec = list(waist= !is.na(hormo)), subset = sex=="Female")
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method
1 Age 3645 0.056* continuous normal
2 Smoking 3645 0.907 categorical
3 Waist circumference 3459 0.007** continuous normal
4 Hormone-replacement therapy 3459 0.898 categorical
selection
1 sex == "Female"
2 sex == "Female"
3 (sex == "Female") & (!is.na(hormo))
4 sex == "Female"
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 변수가 두번나올수도 있다. 예를 들어
compareGroups(group ~ age + sex + bmi + bmi +bmi+ waist + hormo, data=predimed,
selec = list(bmi.1=!is.na(hormo),bmi.2=sex=="Female"))
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method
1 Age 6324 0.003** continuous normal
2 Sex 6324 <0.001** categorical
3 Body mass index 6324 <0.001** continuous normal
4 Body mass index 5661 <0.001** continuous normal
5 Body mass index 3645 0.002** continuous normal
6 Waist circumference 6324 0.045** continuous normal
7 Hormone-replacement therapy 5661 0.850 categorical
selection
1 ALL
2 ALL
3 ALL
4 !is.na(hormo)
5 sex == "Female"
6 ALL
7 ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 이 경우 bmi에 대한 결과가 먼저 모든 대상자(n = 6324)에 대한 결과가 나오고 또한 호르몬대체요법에 missing data가 아닌 (!is.na(hormo)) 사람만을 대상으로 결과가 나온다.
3 연속변수에 대한 Method
통계방법의 변경
디폴트 값으로 연속변수는 정규분포 하는 것으로 분석된다. 테이블이 만들어지는 경우(Section 4.2 createTable 함수 참조 ) 연속변수는 평균과 표준편차로 기술된다.디폴트 값을 바꾸려면 다음과 같이 한다(예를 들어허리둘레를 정규분포하지 않는 것으로 바꾸려면 다음과 같이 한다).
compareGroups(group~age+smoke+waist+hormo,data=predimed, method=c(waist=2))Warning: Cannot compute exact p-value with ties
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.003** continuous normal ALL
2 Smoking 6324 0.444 categorical ALL
3 Waist circumference 6324 0.085* continuous non-normal ALL
4 Hormone-replacement therapy 5661 0.850 categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 Waist circumference의 method 열을 보면 “continuous non-normal”로 표시되어 있다.
Method에는 다음과 같은 값을 사용할 수 있다.
- 1 : 정규분포하는 것으로 분석한다.
- 2 : 연속변수이나 정규분포하지 않는 것으로 분석한다.
- 3 : 범주형 변수로 분석한다.
- NA : Shapiro-Wilks test를 수행하여 정규분포인지 비정규분포인지 결정한다.
어떤 변수의 method에 = NA로 명시되어 있으면 정규분포 여부를 판단하기 위하여 Shapiro-Wilktest를 수행한다. 만일 signifiant threshold를 변경하고 싶으면 다음과 같이 한다.
compareGroups(group~age+smoke+waist+hormo,data=predimed,
method=c(waist=NA),alpha=0.01)Warning: Cannot compute exact p-value with ties
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.003** continuous normal ALL
2 Smoking 6324 0.444 categorical ALL
3 Waist circumference 6324 0.085* continuous non-normal ALL
4 Hormone-replacement therapy 5661 0.850 categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 Cut poinit가 0.01로 명시되어 있으므로 Shapiro-Wilks test 결과 waist circumference는 정규분포로부터 유의하게 떨어져 있으므로 “continuous non-normal”로 수행되었다.
5개 미만의 값을 갖는 연속형변수 : 범주형변수로 취급
Factor 가 아닌 모든 변수는 연속형 변수로 취급된다.단 어떤 변수가 5개 미만의 값을 갖는 경우(5개가 디폴트임 ) 예외로 취급된다. min.dis 값을 바꿈으로써 5개란 값을 바꿀 수 있다.
cuts<-"lo:55=1; 56:60=2; 61:65=3; 66:70=4; 71:75=5; 76:80=6; 81:hi=7"
library(car)
predimed$age7gr<-car::recode(predimed$age, cuts)
compareGroups(group ~ age7gr, data=predimed, method = c(age7gr=NA))Warning: Cannot compute exact p-value with ties
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.007** continuous non-normal ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 역자 주) 나이를 7개의 수준을 갖는 age7gr 변수로 바꾸는데 car패키지에 있는recode함수를 사용하였다 만일 에러가 나는 경우 car패키지를 설치하고 library(car)를 실행하라
compareGroups(group ~ age7gr, data=predimed, method = c(age7gr=NA), min.dis=8)Warning: variable 'age7gr' converted to factor since few different values
contained
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.009** categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 에러를 피하기 위해 반응변수의 최대 범위는 5개 이하로 되어 있다. 만약 어떤 벼누가 5개 이상의 값을 갖는다면 compareGroup은 에러메시지를 출력하게된다. 예를 들어
compareGroups(age7gr ~ sex+bmi+waist, data=predimed)이러한 default 값은 max.ylev 값을 변경하여 바꿀 수 있다.
compareGroups(age7gr ~ sex+bmi+waist, data=predimed,max.ylev=7)
-------- Summary of results by groups of 'Age'---------
var N p.value method selection
1 Sex 6324 <0.001** categorical ALL
2 Body mass index 6324 0.021** continuous normal ALL
3 Waist circumference 6324 0.034** continuous normal ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 마찬가지로 설명변수 level의 최대치도 정해져 있다. 만일 level이 default값을 넘게되면 경고메시지가 출력된다.
compareGroups(group ~ sex + age7gr, method= (age7gr=3), data=predimed, max.xlev=5)Warning: Variables 'age7gr' have been removed since some errors occurred
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Sex 6324 <0.001** categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 4 출력 꾸미기
이 장에서 설명되는 option들은 compareGroups 함수에 해당하지만, 이러한 옵션/세팅의 변화는 creatTable 함수로 표가 만들어지기 전에는 볼수 없다(나중에 설명됨)
Q1, Q3 변수에 대한 method가 “2”일때(즉, 연속변수이며 비정규분포로 분석할때; section 4.1.3 참조) createTable함수를 적용한 후에는 디폴트 값으로 median값과 1, 3사분위수(quartile)이최종결과에 보여진다.
resu1<-compareGroups(group ~ age + waist, data=predimed, method = c(waist=2))
createTable(resu1)
--------Summary descriptives table by 'Intervention group'---------
__________________________________________________________________________
Control MedDiet + Nuts MedDiet + VOO p.overall
N=2042 N=2100 N=2182
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
Waist circumference 101 [94.0;108] 100 [93.0;107] 100 [93.0;107] 0.085
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 허리둘레에 대한 25 및 75 퍼센타일이 계산되었다.대신 2.5와 97.5퍼센타일을 얻으려면
resu2<-compareGroups(group ~ age + smoke + waist + hormo, data=predimed,
method = c(waist=2), Q1=0.025, Q3=0.975)
createTable(resu2)
--------Summary descriptives table by 'Intervention group'---------
___________________________________________________________________________________
Control MedDiet + Nuts MedDiet + VOO p.overall
N=2042 N=2100 N=2182
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
Smoking: 0.444
Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
Waist circumference 101 [80.0;123] 100 [80.0;121] 100 [80.0;121] 0.085
Hormone-replacement therapy: 0.850
No 1811 (98.3%) 1835 (98.4%) 1918 (98.2%)
Yes 31 (1.68%) 30 (1.61%) 36 (1.84%)
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 2.5와 97.5 퍼센타일이 계산되었다. 최소값과 최대값을 얻으려면다음과 같이하면 된다.
resu3=compareGroups(group ~ age + smoke + waist + hormo, data=predimed,
method = c(waist=2), Q1=0, Q3=1)
createTable(resu3)
--------Summary descriptives table by 'Intervention group'---------
___________________________________________________________________________________
Control MedDiet + Nuts MedDiet + VOO p.overall
N=2042 N=2100 N=2182
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
Smoking: 0.444
Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
Waist circumference 101 [58.0;177] 100 [61.0;150] 100 [50.0;176] 0.085
Hormone-replacement therapy: 0.850
No 1811 (98.3%) 1835 (98.4%) 1918 (98.2%)
Yes 31 (1.68%) 30 (1.61%) 36 (1.84%)
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 5. 요약(summary)
createTable 클래스의 객체에 summary 함수를 적용하게 되면 보다 자세한 출력을 얻을수 있다.
res<-compareGroups(group ~ age + sex + smoke + waist + hormo, method = c(waist=2),
data=predimed)
summary(res[c(1,2,4)])
--- Descriptives of each row-variable by groups of 'Intervention group' ---
-------------------
row-variable: Age
N mean sd p.overall p.trend
[ALL] 6324 67.01 6.175
Control 2042 67.34 6.28 0.002666 0.1012
MedDiet + Nuts 2100 66.68 6.016
MedDiet + VOO 2182 67.02 6.213
p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.001672 0.206
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.1727
MedDiet + Nuts
MedDiet + VOO
-------------------
row-variable: Sex
Male Female Male (row%) Female (row%) p.overall p.trend
[ALL] 2679 3645 42.36 57.64
Control 812 1230 39.76 60.24 8.1e-05 0.3884
MedDiet + Nuts 968 1132 46.1 53.9
MedDiet + VOO 899 1283 41.2 58.8
p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.000133 0.3583
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.002076
MedDiet + Nuts
MedDiet + VOO
-------------------
row-variable: Waist circumference
N med Q1 Q3 p.overall p.trend
[ALL] 6324 100 93 107
Control 2042 101 94 108 0.0846 0.03956
MedDiet + Nuts 2100 100 93 107
MedDiet + VOO 2182 100 93 107
p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.1258 0.1106
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.7435
MedDiet + Nuts
MedDiet + VOO 지금은 변수들중 1,3,4번째 변수만 선택되었으므로 나이, 성별, 허리둘레에 대한 결과만 보여진다. 나이는 평균과 표준편차로 요약되며 성별은 빈도와 백분율로, 허리둘레는(method=2) 중앙값과 사분위수가 출력되었다.
6 그래프그리기(Plotting)
변수들은 분포를 보기 위해 plot을 할수 있다.그려지는 Plot은 연속형변수인지 범주형변수인지에 따라 달라진다. Plot은 화면으로 또는 pdf로 저장될수 있다.
plot(res[c(1,2)])
Plot을 할떄 grouping variable에 따라 그림을 그릴수도 있다.(bivar=TRUE)
plot(res[c(1,2)],bivar=TRUE)
7 업데이트(Updating)
compareGroups에서 나온 객체는 나중에 업데이트 될수 있다. 예를 들면
res<-compareGroups(group ~ age + sex + smoke + waist + hormo, data=predimed)
res
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method selection
1 Age 6324 0.003** continuous normal ALL
2 Sex 6324 <0.001** categorical ALL
3 Smoking 6324 0.444 categorical ALL
4 Waist circumference 6324 0.045** continuous normal ALL
5 Hormone-replacement therapy 5661 0.850 categorical ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 객체 res를 업데이트 하면 다음과 같다
res<-update(res, . ~. - sex + bmi + toevent, subset = sex== "Female" ,
method = c(waist=2, toevent=2), selec = list(bmi=!is.na(hormo)))Warning: Cannot compute exact p-value with ties
Warning: Cannot compute exact p-value with tiesres
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method
1 Age 3645 0.056* continuous normal
2 Smoking 3645 0.907 categorical
3 Waist circumference 3645 0.037** continuous non-normal
4 Hormone-replacement therapy 3459 0.898 categorical
5 Body mass index 3459 0.002** continuous normal
6 follow-up to main event (years) 3645 <0.001** continuous non-normal
selection
1 sex == "Female"
2 sex == "Female"
3 sex == "Female"
4 sex == "Female"
5 (sex == "Female") & (!is.na(hormo))
6 sex == "Female"
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 “Sex”는 설명변수목록 에서는 제거했지만 필터로 사용하여 여성 환자만 subset되어 있다. 변수 허리둘레는 “continuous non-normal” 로 바뀌었다. 두개의 새로운 변수가 추가되었는데 Body masss index와 Follow-up이다(continuous non-normal). Body mass index는 호르몬대체요법의 값이 missing vaue가 아닌 참가자의 데이타만 보여지고 있다.
8 Odds Ratio & Harzard Ratios
반응변수가 이원(binary) 데이타라면 최종 테이블에서 Odds ratio(OR) 를 인쇄할수 있다. 반응변수가 time-to-event라면 대신 Harzard ratio(HR)를 인쇄할 수 있다.
ref 이 문장을 통해 참조범주 (reference)를 바꿀수 있다.
res1=compareGroups(htn~age+sex+bmi+smoke,data=predimed,ref=1)
createTable(res1,show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
___________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
Sex: <0.001
Male 595 (54.6%) 2084 (39.8%) Ref. Ref.
Female 494 (45.4%) 3151 (60.2%) 1.82 [1.60;2.08] 0.000
Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
Smoking: <0.001
Never 536 (49.2%) 3356 (64.1%) Ref. Ref.
Current 233 (21.4%) 625 (11.9%) 0.43 [0.36;0.51] 0.000
Former 320 (29.4%) 1254 (24.0%) 0.63 [0.54;0.73] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 위의 예에서는 범주형 자료변수의 참조범주로 첫번쨰의 것이 사용되었다. 즉, 성별에서 남자를 reference로 흡연에서는 비흡연자를 reference로 사용하였다.
res2=compareGroups(htn~age+sex+bmi+smoke,data=predimed,ref=c(smoke=1,sex=2))
createTable(res2,show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
___________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
Sex: <0.001
Male 595 (54.6%) 2084 (39.8%) 0.55 [0.48;0.63] 0.000
Female 494 (45.4%) 3151 (60.2%) Ref. Ref.
Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
Smoking: <0.001
Never 536 (49.2%) 3356 (64.1%) Ref. Ref.
Current 233 (21.4%) 625 (11.9%) 0.43 [0.36;0.51] 0.000
Former 320 (29.4%) 1254 (24.0%) 0.63 [0.54;0.73] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 위에서는 흡연에 대하여는 비흡연자를 reference 로 사용했지만 성별에서는 여성을 reference로 사용했다.
ref.no “ref”의 경우와 마찬가지로 ref.no를 사용하면 yes/no로 되어있는 범주형 변수에서 “no”를 reference로 사용한다
res<-compareGroups(htn ~ age + sex + bmi + hormo + hyperchol, data=predimed,
ref.no="NO" )
createTable(res, show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
________________________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
Sex: <0.001
Male 595 (54.6%) 2084 (39.8%) Ref. Ref.
Female 494 (45.4%) 3151 (60.2%) 1.82 [1.60;2.08] 0.000
Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
Hormone-replacement therapy: 0.856
No 928 (98.4%) 4636 (98.3%) Ref. Ref.
Yes 15 (1.59%) 82 (1.74%) 1.08 [0.64;1.97] 0.773
Dyslipidemia: <0.001
No 409 (37.6%) 1337 (25.5%) Ref. Ref.
Yes 680 (62.4%) 3898 (74.5%) 1.75 [1.53;2.01] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 이떄 ‘no’, “NO”, ‘No’ 중 어떤 것을 사용해도 된다 이 함수는 대소문자를 구분하지 않는다.
fact.ratio 반응변수가 연속형변수인 경우 OR 또는 HR는 변수가 1 증가하는 것을 기준으로 계산된다.
res<-compareGroups(htn ~ age + bmi, data=predimed)
createTable(res, show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
__________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 이때 OR는 나이가 한살, BMI가 1 증가하는 것을 기준으로 계산된다. 나이를 10살 많아질때마다, 그리고 bmi는 2 증가할때마다로 계산하고 싶으면 다음과 같이 fact.ratio를 바꾸어 주면 된다.
res<-compareGroups(htn ~ age + bmi, data=predimed,
fact.ratio=c(age=10,bmi=2))
createTable(res, show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
__________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.43 [1.28;1.59] <0.001 <0.001
Body mass index 28.9 (3.69) 30.2 (3.80) 1.22 [1.17;1.26] <0.001 <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ref.y OR 또는 HR이 계산될때 디폴트 값으로는 반응변수의 첫번째 값이 참조치로 사용된다. 반응변수의 참조치를 바꾸고 싶을때 ref.y를 사용하면 된다.예를 들어
res<-compareGroups(htn ~ age + sex + bmi + hyperchol, data=predimed)
createTable(res, show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
___________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
Sex: <0.001
Male 595 (54.6%) 2084 (39.8%) Ref. Ref.
Female 494 (45.4%) 3151 (60.2%) 1.82 [1.60;2.08] 0.000
Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
Dyslipidemia: <0.001
No 409 (37.6%) 1337 (25.5%) Ref. Ref.
Yes 680 (62.4%) 3898 (74.5%) 1.75 [1.53;2.01] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 이때의 결과물에서 OR는 고혈압이 있는 것에 대한 OR이다. 즉 고혈압이 없는 것이 참조치이다. 하지만 고혈압이 없는 것에 대한 OR를 구하려면 다음과 같이 하면된다.
res<-compareGroups(htn ~ age + sex + bmi + hyperchol, data=predimed, ref.y=2)
createTable(res,show.ratio=TRUE)
--------Summary descriptives table by 'Hypertension'---------
___________________________________________________________________________
No Yes OR p.ratio p.overall
N=1089 N=5235
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Age 65.9 (6.19) 67.2 (6.15) 0.96 [0.98;0.95] <0.001 <0.001
Sex: <0.001
Male 595 (54.6%) 2084 (39.8%) Ref. Ref.
Female 494 (45.4%) 3151 (60.2%) 0.55 [0.48;0.63] 0.000
Body mass index 28.9 (3.69) 30.2 (3.80) 0.91 [0.92;0.89] <0.001 <0.001
Dyslipidemia: <0.001
No 409 (37.6%) 1337 (25.5%) Ref. Ref.
Yes 680 (62.4%) 3898 (74.5%) 0.57 [0.50;0.65] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 이때에는 “no hypertension”에 대한 OR을 구한것으로 고혈압이 있는 것을 참조치로 하여 계산되었다.
반응 변수가 생존분석에 사용되는 Surv클래스의 변수라면 bivariateplot을 그릴떄 설명변수가 범주형 변수라면 Kaplan-Meier Figure가 그려진다. 설명변수가 연속형 변수인 경우에는시간에 따른 그래프가 그려진다.
plot(compareGroups(tmain ~ sex,data=predimed), bivar=TRUE)
plot(compareGroups(tmain ~ age,data=predimed), bivar=TRUE) 
9 Time-to-event 설명변수
Surv클래스의 변수가 설명변수로 사용될 경우 주어진 시간 까지 Kaplan-Meier방법으로 계산된 event의 발생확률(probability) 가 계산된다.
res<-compareGroups(sex ~ age + tmain, timemax=c(tmain=3),
data=predimed)
res
-------- Summary of results by groups of 'Sex'---------
var N p.value method selection
1 Age 6324 <0.001** continuous normal ALL
2 AMI, stroke, or CV Death 6324 <0.001** Surv [Tmax=3.00] ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 이떄 tmain은 3년으로 계산되었다.
Surv클래스의 변수에 plot함수를 적용하면 Kaplan-Meier 그림이 그려진다. 이 그림은 grouping variable에 따라
계층화되어 그려질 수 있다.
plot(res[2])
plot(res[2],bivar=TRUE)