survival curve plotting

[cook]

일주일째 고민하다가 글 올려봅니다.

cox phm으로 survival curve plotting하려고 하는데요

한변수가 4개의 카테고리인 경우 곡선이 4개가 그려지는데요

(카테고리는 2,3,4,5 로 주었습니다.)

R을 이용해서 하려고 합니다...

km로 그린 survival curve와 비슷하게 나와야하는데 그렇지가 않아요

제가 R에서 coding 한 부분 좀 봐주시고 지적해주시면 좋겟어요

꼭 부탁드려요 아시는 분...

library(survival)

dset=read.table("생략",header=T)
attach(dset)

S=Surv(time,death)

phmcelltype=coxph(S~celltype)
d.phmcelltype=coxph.detail(phmcelltype)
times=c(0,d.phmcelltype$t)
h0=c(0,d.phmcelltype$hazard)
s0=exp(-cumsum(h0))
beta=phmcelltype$coef
celltype=as.numeric(celltype)
x1=c(1)-mean(celltype,na.rm=T)  #adeno
sx1=s0^exp(t(beta) %*% x1)
x2=c(2)-mean(celltype,na.rm=T) #large
sx2=s0^exp(t(beta) %*% x2)
x3=c(3)-mean(celltype,na.rm=T) #small
sx3=s0^exp(t(beta) %*% x3)
x4=c(4)-mean(celltype,na.rm=T)  #squamous
sx4=s0^exp(t(beta) %*% x4)

survfit1=survfit(S~celltype)
plot(survfit1,lty=2,col=1:4,xlab='t',ylab=expression‎(hat(S)(t)))
par(new=TRUE)
lines(times,sx1,type='s',col=1)
lines(times,sx2,type='s',col=2)
lines(times,sx3,type='s',col=3)
lines(times,sx4,type='s',col=4)

[안재형]

데이터를 올리시는게 나을것 같은데요...

[cook]

http://www.sph.emory.edu/~dkleinb/surv2datasets/vets.dat

원자료주소구요 ~여기서 첫번째칼럼이 5개의 digit으로 되어있는데 첫번째digit은 치료방법이구요 두버째 digit~다섯번째 digit까지가 종양크기 같아요 ~종양크기순으로 large,adeno,small,squamous가 한셀에 되어잇던걸 분리해서 celltype으로 만들었어요
두번째digit이 1이면 celltype=2로 ,세번째 digit이 1이면 celltype=3로,네번째 digit이 1이면 celltype=4로,다섯번째 digit이 1이면 celltype=5로 ...그리고 celltype 변수가 생존에 영향을 주는지 보구 plot그리려구요
근데 plot이 말씀드린대로 km curve에 적합되지 않아요

[지니]

제가 정확한 의도는 모르겠는데
위 coding시 plot이 안그려져서 coding일부 수정하였는데 plot이 나옵니다
원하는 내용이 맞는지는 모르겠네요

아래 부분을 수정하시면 plot이 나옵니다
plot(survfit,col=1:4,lty=2,xlab='t',ylab=expression!(hat(S)(t))) → plot(survfit,col=1:4,lty=2,xlab='t',ylab='expression',(hat(S)(t)))

[지니]

plot은 그려지는데 hat(S)(t) 부분에서 에러가 있네요
plot(survfit,col=1:4,lty=2,xlab='t',ylab='expression',(hat(s0)))
이렇게 바꾸니 에러는 없이 plot이 그려지는데, 맞는지는 한번 확인해보셔야 할것 같네요

[cook]

kmplot은 그려지는데요..
여기에 cox phm을 적합시키려구요..
두 그래프가 거의 일치하게 나와야 하는데
cox 예측 survival curve가 다르게 그려집니다.
말씀하신 부분은 카페글쓰기 오류인거 같아요 !기호가 없는데 생겨있네요^^

[안재형]

저도 "!"때문에 애 좀 먹었습니다. 안 없어지더라구요.

[cook]

그래요? 왜그렇죠? 소스에는 없는 문자에요
그 plot부분보다는 소스중에 meanx ~ sx5까지가 머가 잘못인거 같은데..모르겟어요
다른부분은 잘못된 곳이 없나요?

[cook]

위의것중에 phmcelltype=coxph(S~factor(celltype),data=dset) 이렇게 고쳐봤는데요..
저자님 책도 봤는데..birthwt 데이터의 저체중 분석에서 race가 3요인으로 나뉘잖아요..
summary결과엔 왜 race2,3만 나오는거죠?

[Au revoir]

제가 해보니 대충 비슷하게 나오는 것 같습니다. 'survfit' function을 이용하시면 직접 Cox model을 적합했을시의 Survival function estimate들을 구하실 수 있습니다. 본문에 올리신 맨 마지막 줄을 다음의 code로 바꾸어서 돌려보시기 바랍니다.
lines(survfit(phmcelltype, newdata=data.frame(celltype="squamous"),conf.type="none"),col=2)
lines(survfit(phmcelltype, newdata=data.frame(celltype="smallcell"),conf.type="none"),col=3)
lines(survfit(phmcelltype, newdata=data.frame(celltype="adeno"),conf.type="none"),col=4)
lines(survfit(phmcelltype, newdata=data.frame(celltype="large"),conf.type="none"),col=5)

[cook]

그래프가 의도대로 안되고 있어요..어디가 문제인지 지적좀 부탁드려요

[Au revoir]

일단 coxph를 이용해서 구한 object (phmcelltype)를 survfit function안에 넣으시면 직접 survival function estimates를 구할 수가 있습니다. 원문에 있는 code중에서 d.phmcelltype = ... 부터 sx4 = ... 까지는 굳이 하실 필요가 없습니다. 마지막 4줄을 위의 댓글에 있는 code로 바꾸셔도 돌려보시면 됩니다. 다음의 링크를 참조하시면 좋을 듯 합니다: http://stat.ethz.ch/R-manual/R-devel/library/survival/html/survfit.coxph.html

[cook]

suvfit 함수를 이용해서 하는게 아니구요..
coxph.detail함수를 이용해서 baseline survival function를 구해서 exp(bx)를 제곱하면 survival function이 구해지면 그거로 그래프를 그리려는게 목적이에요..저의 교재 내용이 그래요~^^