6.4 모수적 방법을 이용한 생존함수의 추정과 비교
공학(시멘트의 양, 유리의 버티는 힘), 경영(고객 수), 교통(소방차 수) 모두 모수적 방법을 이용.
분포를 이루기 때문에 많은 분야에서 사용된다.
1. 데이터 입력
> setwd('c:/Rwork')
> lung=read.table('lung.txt',header=T)
> head(lung)
癤퓍ime status x1 x2 x3 trt type
1 411 1 70 64 5 1 1
2 126 1 60 63 9 1 1
3 118 1 70 65 11 1 1
4 82 1 40 69 10 1 1
5 8 1 40 63 58 1 1
6 25 0 70 48 9 1 1
> colnames(lung)<-c('time','status','x1','x2','x3','trt','type')
> head(lung)
time status x1 x2 x3 trt type
1 411 1 70 64 5 1 1
2 126 1 60 63 9 1 1
3 118 1 70 65 11 1 1
4 82 1 40 69 10 1 1
5 8 1 40 63 58 1 1
6 25 0 70 48 9 1 1
> attach(lung)
2. 모수적 모형에 근거한 생존시간 분석
> library(survival)
> weibull = survreg(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='weibull') #공변량factor(trt), 처리효과factor(type), 와이블 분포weibull, gaussian(정규분포), logistic(로지스틱)eh rksmd.
> summary(weibull)
Call:
survreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +
factor(type), data = lung, dist = "weibull")
Value Std. Error z p
(Intercept) 1.06044 1.35959 0.780 4.35e-01
x1 0.05420 0.00954 5.680 1.35e-08
x2 0.01168 0.01918 0.609 5.42e-01
x3 0.00379 0.01051 0.361 7.18e-01
factor(trt)2 0.28871 0.36899 0.782 4.34e-01
factor(type)2 -0.49964 0.45323 -1.102 2.70e-01
factor(type)3 -1.25968 0.49732 -2.533 1.13e-02
factor(type)4 -0.40243 0.38726 -1.039 2.99e-01
Log(scale) -0.13615 0.13146 -1.036 3.00e-01
Scale= 0.873
Weibull distribution
Loglik(model)= -203.4 Loglik(intercept only)= -219.7
Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05
Number of Newton-Raphson Iterations: 6
n= 40
> weibull
Call:
survreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +
factor(type), data = lung, dist = "weibull")
Coefficients:
(Intercept) x1 x2 x3 factor(trt)2
1.060436846 0.054195931 0.011681287 0.003792838 0.288708242
factor(type)2 factor(type)3 factor(type)4
-0.499640589 -1.259681146 -0.402431957
Scale= 0.8727099
Loglik(model)= -203.4 Loglik(intercept only)= -219.7
Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05
n= 40
logT = 1.060 + 0.054x1 + 0.012x2 + 0.004x3 + 0.289x4 - 0.500x5 -1.260x6 - 0.402x7 +0.873e
factor(trt)2 0.28871 0.36899 0.782 4.34e-01
-> x4는 처리그룹을 나타내는 가변수(dummy variable)로서 처리가 standard일 때 1, test 일 때 2의 값을 갖는다.
p-value가 0.434로 유의하지 않으므로 두 처리(standard, test)는 생존시간의 차이를 보이지 않는다.
type은 squamous,small,adeno,large, x5,6,7은 종양의 유형을 가르키는 가변수로서
x5 = 1, if 'small', = 0 o.w.
x6 = 1, if 'adeno', = 0 o.w.
x7 = 1, if 'large', = 0 o.w.
따라서 종양이 squamous인 경우에는 x5=x6=x7 = 0 된다.
factor(type)3 -1.25968 0.49732 -2.533 1.13e-02
-> adeno의 p-value만 0.0113으로 0.05보다 적으므로 유의, squamous와 차이를 보인다고 할 수 있다. adeno의 회귀계수가 -1.25968로 이 유형을 가진 사람들은 상대적으로 생존시간이 짧다는 것을 알 수 있다. type2 small, type4 large의 경우 p-value가 각각 0.270, 0,299로 유의하지 않다, 즉 squamous와 차이를 보이지 않는다.
Loglik(model)= -203.4 Loglik(intercept only)= -219.7
Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05
2(logL-logL0) = 32.59 > 14.067 모든 공변량의 회귀계수가 0이라는 귀무가설을 아주 강하게 기각한다.
통상적인 선형모형에서 모형에 대한 F-검정을 하는 것과 같이 여기서도 모든 공변량의 회귀계수가 0이라는 가설에 대하 우도비 검정(likelihood-ratio ttest)를 할 수 있다. 우리가 고려한 모형에서의 로그-우도(log-likelihood)를 logL, 공변량이 전혀 없는 귀무모형에서의 로그-우도를 logL0라고 했을 때 2(logL-logL0)이 귀무가설하에서 근사적으로 x2-분포를 따르게 되므로 이를 이용하여 검정할 수 있다.
이때 자유도는 귀무모형에서 제외되는 공변량의 수와 같게 되며, 위에서는 7이 된다.
LogL: 고려한 모형에서의 로그-우도 Log Likelihood for WEIBULL -203.4는 이 모형에 대한 로그-우도이다.
LogL0: 공변량이 전혀 없는 귀무모형에 대한 로그-우도를 구해보면 Log Likelihood for WEIULL = =219.7을
따라서 Chisq= 32.59는 아래와 같은 계산을 통해 얻을 수 있으며, 이 값이 x2-분포의 임계치인 x2 0.95(7) = 14.067보다 크므로
5%유의수준하에서 모든 공변량의 회귀계수가 0이라는 귀무사설을 기각하게 된다.
> 2*(-203.4-(-219.7))
[1] 32.6
3. 모수적 모형의 적합도 검토
고려한 모형이 타당한가를 검토하는 방법 중 하나는 로그-우도 비교. R을 이용하여 각 모형에 대하 로그-우도를 출력한 다음 이들을 비교하여 절대값이 가장 큰 모형을 택할 수 있다. 또는 AIC(Akaike information criterion) 값을 비교하여 이 값이 더 작은 것을 선택할 수 있다.
> library(flexsurv) #flexsurv library: Flexible parametric survival models
Warning message:
package ‘flexsurv’ was built under R version 3.2.5
> gengamma=flexsurvreg(formula=Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='gengamma')
> gengamma
Call:
flexsurvreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +
factor(type), data = lung, dist = "gengamma")
Estimates:
data mean est L95% U95% se exp(est)
mu NA 1.16226 -1.69652 4.02104 1.45859 NA
sigma NA 0.82367 0.47903 1.41628 0.22778 NA
Q NA 1.19119 -0.35449 2.73687 0.78863 NA
x1 56.50000 0.05426 0.03598 0.07254 0.00933 1.05576
x2 56.57500 0.01210 -0.02543 0.04963 0.01915 1.01217
x3 15.65000 0.00494 -0.01741 0.02729 0.01141 1.00495
factor(trt)2 0.50000 0.27185 -0.47815 1.02186 0.38266 1.31239
factor(type)2 0.27500 -0.57430 -1.63425 0.48566 0.54080 0.56310
factor(type)3 0.12500 -1.36101 -2.57769 -0.14434 0.62076 0.25640
factor(type)4 0.25000 -0.48503 -1.45929 0.48923 0.49708 0.61568
L95% U95%
mu NA NA
sigma NA NA
Q NA NA
x1 1.03664 1.07524
x2 0.97489 1.05088
x3 0.98274 1.02767
factor(trt)2 0.61993 2.77835
factor(type)2 0.19510 1.62524
factor(type)3 0.07595 0.86560
factor(type)4 0.23240 1.63107
N = 40, Events: 37, Censored: 3
Total time at risk: 5784
Log-likelihood = -203.4059, df = 10
AIC = 426.8117
> weibull=flexsurvreg(formula=Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='weibull')
> weibull
Call:
flexsurvreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +
factor(type), data = lung, dist = "weibull")
Estimates:
data mean est L95% U95% se exp(est)
shape NA 1.14586 0.88555 1.48269 0.15066 NA
scale NA 2.88763 0.20141 41.39933 3.92317 NA
x1 56.50000 0.05420 0.03551 0.07288 0.00953 1.05569
x2 56.57500 0.01168 -0.02587 0.04924 0.01916 1.01175
x3 15.65000 0.00379 -0.01679 0.02438 0.01050 1.00380
factor(trt)2 0.50000 0.28871 -0.43443 1.01185 0.36896 1.33470
factor(type)2 0.27500 -0.49964 -1.38783 0.38855 0.45317 0.60675
factor(type)3 0.12500 -1.25968 -2.23430 -0.28506 0.49726 0.28374
factor(type)4 0.25000 -0.40243 -1.16137 0.35651 0.38722 0.66869
L95% U95%
shape NA NA
scale NA NA
x1 1.03615 1.07560
x2 0.97446 1.05047
x3 0.98335 1.02468
factor(trt)2 0.64763 2.75068
factor(type)2 0.24962 1.47484
factor(type)3 0.10707 0.75197
factor(type)4 0.31306 1.42833
N = 40, Events: 37, Censored: 3
Total time at risk: 5784
Log-likelihood = -203.4363, df = 9
AIC = 424.8727
일반화감마분포 gengamma분포의 경우 로그-우도 값이 -203.4059, AIC는 426.8117
와이블분포 wibull 분포의 경우 로그-우도 값이 -203.4363, AIC는 424.8727로 큰 차이는 없다. 그래프로 그려보면 아래와 같이 큰 차이가 없음을 알 수 있다.
데이터 시각화
> plot(weibull,xlab='time',ylab='survival function',ci=F)
> lines(gengamma,col='blue',lty=2,ci=F)
> legend('topright',c('weibull','generalized gamma'),lty=1:2,col=c('red','blue'))
출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저)
'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글
| 제6.5장 준모수적 방법 (0) | 2017.01.26 |
|---|---|
| 제6.3장 비모수적 방법을 이용한 생존함수의 비교 (0) | 2017.01.26 |
| 제6.2장 비모수적 방법 - 2. 누적한계추정법 (0) | 2017.01.26 |
| 제6.2장 비모수적 방법 - 1. 생명표 방법 (0) | 2017.01.26 |
| 제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우 (2) | 2017.01.13 |
