책에 SAS 코드 추가하기

[안재형]

책 "R을 이용한 누구나 하는 통계분석"이 R책이어서 회사에서 좀 어려움을 겪다가 잘 해결되었었습니다.

몇주전에 매니저가 SAS legal department에서 연락이 왔다고 하더군요.

legal department에서 다음판에 책에 SAS 결과물을 추가하라네요.

정확히 얘기하면 다음판이나 책이 다른 나라로 번역될 경우라는데, 영어로 번역할 계획은 전혀 없으니 다음판이겠죠^^

그래서 몇주동안 틈틈히 돌려봤습니다. 

사분위수는 다른 방법을 사용해서 SAS와 R이 약간 차이가 있고, 

survival에서 Cox regression에서 conditional likelihood를 maximize하는 방법이 다른지 약간 다른 결과가 나오네요.

책에 R결과와 동일한 결과를 얻을수 있는 SAS 코드를 넣어 데이터와 내용을 자료실에 올려놓았습니다.

SAS-R.pdf에 있는 내용들을 책에 각 chapter 끝부분에 추가할 생각입니다. 

sample size는 proc power로 해야하는데 영 번거롭네요. 

이 부분은 SAS가 강점을 가지는 부분인듯합니다.

회사의 의도는 조금이라도 SAS의 사용을 늘리려는것이겠지만, R의 결과가 SAS와 동일하니,

연구비를 절약해야하는 우리나라에서는 제 책이 SAS -> R로 전환시키는 Rosetta stone이 되지나 않을까 합니다.

우리나라뿐 아니라 미국에서도 제약회사, 큰 금융기관등 안정성이 요구되고, 이미 SAS로 모든 게 set-up 되어 있는 곳 빼고는 점점 R이 SAS를 대체할 것 같다는 느낌이 듭니다.

약간의 Bug나 version up할때마다 가끔씩 생기는 문제는 그리 큰 문제가 될 것 같지는 않습니다.

다만 대용량 데이타의 경우 SAS가 큰 강점인데, 언젠가는 R도 따로오겠지요.

추가: Survival 전공한 Au revoir님 도움으로 Cox regression은 method="breslow"로 하니 SAS와 같은 결과가 나옵니다.

> coxph(Surv(time,status==1)~x,data=aml,method="breslow")

SAS에서는 ties=efron이라고 하면 책과 같은 결과가 나옵니다.

proc phreg data=aml;
  class x (ref="Maintained");
  model time*status(0) = x / ties=efron;
run;

SAS에서 proc sql로 중복을 찾아보니 2개가 있네요.

proc sql;

  select time, status, x, count(x)

  from aml

  group by time, status, x

  having count(x) ge 2

;quit;

                             time        status  x

                     -----------------------------------------------------

                                5             1  Nonmaintained           2

                                8             1  Nonmaintained           2

[cardiomoon]

열심히 하시는 모습이 보기 좋습니다.
저도 번역이 많이 되어가고 있긴 한데 올리려니 번거롭군요. 하지만 정말 공부할수록 재미 있는 것 같습니다.

[안재형]

저야 제 전공인데요^^ 전공도 아닌 분야를 열심히 공부하시는 문선생님이 더 대단하시죠.
이제 슬슬 출판사 알아보셔야겠는데요. 책에 그림 올리는게 상당히 번거로울겁니다.

[강성찬]

일하시면서 귀찮으실텐데 그래도 도움은 많이 되겠네요. 제가 소속된 곳의 의대 교수님은 SAS만 진정한 통계패키지라 생각하시는 분이라서요...ㅎㅎ

[안재형]

미국에서도 마찬가지입니다. R이 SAS하고 결과가 다르거나 없는 결과가 있으면 난리를 피우는 사람이 있더군요.

SAS로서는 그런 SAS에 대한 사대주의정신이 도움이 되겠죠^^

[Au revoir]

죄송한데요 COX REGRESSION의 경우에 차이가 얼마나 나는지요? 그냥 maximization방법의 차이에서오는 오차정도라고 생각할 정도인지 아니면 그 이상인지요?

[안재형]

이유는 모르겠는데 좀 차이가 나던데 오차라고 보기에는 좀 커... 그 데이터가 R에 내장되어 있는 데이터니까 한번 돌려봐.

[Au revoir]

옹... 신기하군요. 함 해봐야겠습니다.

[안재형]

beta coefficient가 R은 0.9155, SAS는 0.90422이니 꽤 다르지...

[Au revoir]

Tied observation때문인것 같아요. R과 SAS의 default가 다른 것으로 알고 있거든요 (R은 Efron, SAS는 Breslow). 혹, SAS의 PROC PHREG나 R의 COXPH사용시 따로 Tied observation을 처리하는 방법을 지정해주지 않았다면 이것때문일 것 같네요. alm.csv를 보니 일단 time=8, status=1, x='Nonmaintained'인 observation이 두개 있네요. 이따가 직접 돌려봐서 확인해보도록 하죠.

[안재형]

그렇군... 역시 전공자는 달라^^ 고마워~

[안재형]

> coxph(Surv(time,status==1)~x,data=aml,method="breslow") 라고 하니 SAS하고 똑같이 나오네^^

[Au revoir]

앗, 직접 해보셨군요. 문제는 해결이 된듯하네요!

[안재형]

고마와~