데이터마이너를 꿈꾸며

6.5 준모수적 방법

Cox의 비례 모형

모수적 모형은 가정이 타당할 때는 상당히 효율적이지만 어떤 모형이 적당한가에 대한 지식이 없다면 함부로 사용하기가 곤란하다. 이에 반해 Cox(1972)의 비례 위험 모형(proportional hazards model)은 준모수적(semi-parametric)방법으로서 생존 시간의 분포에 대한 가정을 필요로 하지 않는다. 또한 시간에 따라 바뀌는 공변량(time-dependent variable)의 경우에도 분석할 수 있다는 장점이 있어, 생존자료의 분석에 매우 자주 사용된다.

[예] Prenctice(1973)에 소개된 것으로 40명의 폐암 환자의 생존시간을 조사한 것이다. 40명의 환자 중 21명은 기존 치료방법인 처리1에, 나머지 19명은 새로운 치료방법은 처리2에 할당되었으며, 생존시간에 영향을 미칠 것으로 생각되는 공변량은 다음과 같다.

X1: 진단시의 환자상태(Perfermance Status : 0~100 점)

X2: 환자의 나이(단위 : 년)

X3: 진단 후 연구 참여시까지의 시간(단위 : 월)

Trt: 치료 방법 (1 – 기존 치료, 2 – 새로운 치료)

TYPE: 종양의 유형 (squamous, small, adeno, large)

1. 데이터 입력

> setwd('c:/Rwork/')

> lung<-read.table('lung.txt',header=T)

> 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

2. Cox비례위험모형에 근거한 생존시간 분석

> library(survival)

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> summary(coxfit1)

Call:

coxph(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung)

  n= 40, number of events= 37

                   coef exp(coef)  se(coef)      z Pr(>|z|)   

x1            -0.060281  0.941500  0.013777 -4.375 1.21e-05 ***

x2            -0.015086  0.985027  0.022340 -0.675   0.4995   

x3             0.001201  1.001201  0.011886  0.101   0.9195   

factor(trt)2  -0.448171  0.638795  0.431302 -1.039   0.2988   

factor(type)2  0.279682  1.322709  0.547259  0.511   0.6093   

factor(type)3  1.418190  4.129638  0.625283  2.268   0.0233 * 

factor(type)4  0.361145  1.434971  0.479210  0.754   0.4511   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95

x1               0.9415     1.0621    0.9164    0.9673

x2               0.9850     1.0152    0.9428    1.0291

x3               1.0012     0.9988    0.9781    1.0248

factor(trt)2     0.6388     1.5654    0.2743    1.4876

factor(type)2    1.3227     0.7560    0.4525    3.8663

factor(type)3    4.1296     0.2422    1.2125   14.0655

factor(type)4    1.4350     0.6969    0.5610    3.6707

Concordance= 0.764  (se = 0.058 )

Rsquare= 0.524   (max possible= 0.994 )

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

- coxph: cox의 비례위험 모형을 ㅈ거합함.

-> 결과 해석:

- 진단시의 환자상태(x1)에 대한 회귀계수가 -0.060281로 유의하다(p-value < 0.001). 위험률(exp(coef))를 보면 진단시의 환자상태에 대한 점수가 1점 더 높은 환자는 다른 조건이 동일한 환자와 비교할 때 0.9415(exp(-0.060281))배 위험하다는 것을 알 수 있다.

exp(-coef) = 1.0621은 exp(coef)의 역수인 1/0.9415이다.

- 나이(x2), 진단 후 연구 참여시까지의 시간(x3), 처리(trt)는 모두 p-value가 0.05보다 커서 유의한 영향을 미치지 않는다.

- TYPE: 종양의 유형 (squamous, small, adeno, large)을 보면 factor(type)3인 adeno의 경우 squamous와 차이를 보인다.(p-value: 0.0233)

회귀계수가 0이라는 가설에 대한 검정통계량. 차례로 우도비 검정(-2 LOG L), Wald 검정통계량, Score 검정이 있다. 세 통계랑 모두 가설을 기각하게 되므로 (p-value < 0.05), 적어도 하나의 공변량은 의미가 있다는 것(= 적어도 모든 변수가 0이라는 귀무가설을 기각)을 알 수 있다.

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

3. 생존함수추정치

> fit4 = survfit(coxfit1)

> summary(fit4)

Call: survfit(formula = coxfit1)

 time n.risk n.event survival  std.err lower 95% CI upper 95% CI

    1     40       1 9.90e-01 1.10e-02     9.68e-01        1.000

    2     39       1 9.79e-01 1.67e-02     9.47e-01        1.000

    8     38       2 9.54e-01 2.71e-02     9.03e-01        1.000

   10     36       1 9.37e-01 3.32e-02     8.74e-01        1.000

   11     35       1 9.20e-01 3.87e-02     8.47e-01        0.999

   12     34       2 8.83e-01 4.84e-02     7.93e-01        0.984

   15     32       1 8.64e-01 5.33e-02     7.65e-01        0.975

   16     31       1 8.44e-01 5.79e-02     7.38e-01        0.965

   18     30       1 8.22e-01 6.27e-02     7.07e-01        0.954

   19     29       1 7.96e-01 6.77e-02     6.74e-01        0.941

   20     28       1 7.67e-01 7.27e-02     6.37e-01        0.924

   21     27       1 7.34e-01 7.73e-02     5.97e-01        0.903

   43     25       1 6.97e-01 8.22e-02     5.54e-01        0.879

   44     24       1 6.61e-01 8.62e-02     5.12e-01        0.854

   51     23       1 6.26e-01 8.94e-02     4.73e-01        0.828

   54     22       1 5.84e-01 9.23e-02     4.28e-01        0.796

   56     21       1 5.44e-01 9.39e-02     3.87e-01        0.763

   82     20       1 5.06e-01 9.46e-02     3.50e-01        0.730

   84     19       1 4.64e-01 9.54e-02     3.11e-01        0.694

   90     18       1 4.25e-01 9.51e-02     2.74e-01        0.659

  100     17       1 3.81e-01 9.46e-02     2.34e-01        0.620

  118     15       1 3.35e-01 9.36e-02     1.94e-01        0.580

  126     14       1 2.93e-01 9.12e-02     1.59e-01        0.540

  153     13       1 2.53e-01 8.79e-02     1.28e-01        0.500

  164     12       1 2.14e-01 8.37e-02     9.97e-02        0.461

  177     11       1 1.80e-01 7.81e-02     7.66e-02        0.421

  200     10       1 1.40e-01 7.08e-02     5.20e-02        0.377

  201      9       1 1.07e-01 6.19e-02     3.41e-02        0.333

  231      8       1 8.00e-02 5.25e-02     2.21e-02        0.289

  250      6       1 5.18e-02 4.19e-02     1.06e-02        0.253

  287      5       1 2.88e-02 2.98e-02     3.81e-03        0.218

  340      4       1 9.18e-03 1.50e-02     3.74e-04        0.225

  411      3       1 2.19e-03 5.10e-03     2.30e-05        0.210

  991      2       1 1.28e-04 5.20e-04     4.54e-08        0.362

  999      1       1 3.36e-10 5.39e-09     7.43e-24        1.000

 > names(fit4)

 [1] "n"         "time"      "n.risk"    "n.event"   "n.censor"  "surv"     

 [7] "type"      "cumhaz"    "std.err"   "upper"     "lower"     "conf.type"

[13] "conf.int"  "call"     

> fit4$surv

 [1] 9.895862e-01 9.787704e-01 9.543895e-01 9.372414e-01 9.195316e-01

 [6] 8.834541e-01 8.637107e-01 8.439345e-01 8.215162e-01 7.962148e-01

[11] 7.670584e-01 7.344065e-01 7.344065e-01 6.974181e-01 6.609451e-01

[16] 6.256915e-01 5.835925e-01 5.436422e-01 5.055202e-01 4.643798e-01

[21] 4.249339e-01 3.809762e-01 3.809762e-01 3.352450e-01 2.932721e-01

[26] 2.533265e-01 2.142532e-01 1.795911e-01 1.401508e-01 1.065182e-01

[31] 8.001270e-02 5.176743e-02 2.881620e-02 9.176925e-03 2.194374e-03

[36] 1.281990e-04 3.359491e-10

4. 생존함수그래프

> plot(survfit(coxfit1),xlab='time',ylab='Survival function',xlim=c(0,998.9))

> legend(500,1.0,c('누적한계추정치','95%신뢰구간'),lty=c(1,2))

5. 누적함수그래프

누적위험함수의 추정치 그래프를 통해 위험함수의 형태를 짐작할 수 있다. 예를 들어 단조 증가하는 직선형태는 위험함수가 시간에 대해 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 위험함수가 증가하고, 아래 방향으로 휘면 감소한다는 것을 의미한다.

> H.hat = -log(fit4$surv)

> H.hat = c(H.hat,tail(H.hat,1)) 

> plot(c(fit4$time,1100),H.hat,xlab='time',ylab='comulative hazard function',type='s')

 - tail: 벡터, 매트릭스 데이터에서 마지막 n개의 행들을 선택함. 예를 들어 tail(H.hat,1)라고 입력하면 H.hat 벡터의 마지막 1개의 성분을 선택. 해당 처리를 하는 이유는 위험함수 곡선은 시간에 대한 상승곡선이므로, 마지막 함수값을 추가함을써 곡선의 모양을 자연스럽게 하기 위한 처리임.

6. 비례성 검토를 위한 로그-로그 그림

비례성의 가정이 타당한 것인지 검토하는 방법에 대해 알아보자. 폐암자료에서 두 처리그룹별로 시간에 따른 log(-log S^ (t)) 의 그래프를 그렸을 때 평행하게 되는 가를 볼 수 있다.

> coxfit2 = coxph(Surv(time,status)~x1+x2+x3+strata(trt)+factor(type),data=lung) 

> plot(survfit(coxfit2),fun='cloglog',lty=1:2,col=c('red','blue'))#fun='cloglog', 두 그래프가 대체적으로 평행하므로 비례성이 타당하다고 할 수 있다.

> legend('topleft',c('처리1','처리2'),lty=1:2,col=c('red','blue'))

- strata(trt): 처리를 층으로 입력해주면, 처리를 층으로 입력하여, 처리1과 처리2로 나누어 Cox 비례 위험 모형을 추정함.

- fun='cloglog'를 입력하면 로그-로그 그림을 그릴 수 있다.

- 비례성검토를 위한 로그-로그 그림(log-log plot)은 그래프가 평행하므로 비례성의 가정이 타당하다 할 수 있다.

7. 로그 – 랭크 테스트 log – rank test

1) 종양 유형 별 생존함수의 차이

> survdiff(Surv(time,status)~factor(type),data=lung)

Call:

survdiff(formula = Surv(time, status) ~ factor(type), data = lung)

                N Observed Expected (O-E)^2/E (O-E)^2/V

factor(type)=1 14       12    16.72     1.330     2.825

factor(type)=2 11       10     6.93     1.356     1.771

factor(type)=3  5        5     2.18     3.651     4.099

factor(type)=4 10       10    11.17     0.123     0.187

 Chisq= 7.4  on 3 degrees of freedom, p= 0.0614

X²(df=3)=7.4이고, p-value가 0.0614로 0.05보다는 약간 크지만 상당히 의미 있음을 알 수 있다.

2) 치료법 생존함수의 차이

> survdiff(Surv(time,status)~x1,data=lung)

Call:

survdiff(formula = Surv(time, status) ~ x1, data = lung)

      N Observed Expected (O-E)^2/E (O-E)^2/V

x1=20 2        2    0.128   27.3119   28.6351

x1=30 4        4    1.565    3.7893    4.2568

x1=40 7        7    1.932   13.2948   15.1966

x1=50 4        3    3.298    0.0270    0.0305

x1=60 7        7    6.642    0.0193    0.0247

x1=70 9        7   12.420    2.3650    3.7914

x1=80 6        6    6.982    0.1382    0.1761

x1=90 1        1    4.033    2.2811    3.7772

 Chisq= 58.8  on 7 degrees of freedom, p= 2.65e-10

X²(df=7)=58.8이고, p-value가 0.05보다 매우 작기 때문에 의미 있음을 알 수 있다.

8. anova를 이용하여 모형 비교

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> coxfit2 = coxph(Surv(time,status)~x1+factor(type),data=lung)

> anova(coxfit2,coxfit1)

Analysis of Deviance Table

 Cox model: response is  Surv(time, status)

 Model 1: ~ x1 + factor(type)

 Model 2: ~ x1 + x2 + x3 + factor(trt) + factor(type)

   loglik Chisq Df P(>|Chi|)

1 -88.143                  

2 -87.516 1.254  3    0.7401

anova()는 Cox Regression의 모형을 비교할 때 LRT(Likelihood ratio test)를 사용한다. 환자의 나이(x2), 진단 후 연구 참여시까지의 시간( x3)의 p-value는 0.7401으로 유의하지 않다.

출처: 보건정보데이터 분석 (이태림, 이재원, 김주한, 장대흥 공저), R을 이용한 누구나 하는 통계분석

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'))

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저)

6.3 비모수적 방법을 이용한 생존함수의 비교

[예] 흑색종(melanoma) 환자들에 대한 BCG와 CP(coryne-bacterium parvum)의 생존지속 효과를 비교하기 위한 연구에서 30명의 흑색종 환자 중 11명은 BCG 처리를 받고 나머지 19명은 CP처리를 받았다고 한다. 중도절단이 포함되어 있는 경우에 이 두 그룹의 생존분포를 비교하기 위한 방법을 알아보자.

1. 데이터 입력

> library(survival)

> setwd('c:/Rwork')

> melanoma = read.table('melanoma.txt',header=T)

> head(melanoma,3)

  癤퓍ime status   x

1    33.7      0 BCG

2     3.9      1 BCG

3    10.5      1 BCG

> colnames(melanoma)<-c('time','status','x')

> head(melanoma)

  time status   x

1 33.7      0 BCG

2  3.9      1 BCG

3 10.5      1 BCG

4  5.4      1 BCG

5 19.5      1 BCG

6 23.8      0 BCG

> attach(melanoma)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit2 = survfit(Surv(time,status)~x,data=melanoma)

> summary(fit2)

Call: survfit(formula = Surv(time, status) ~ x, data = melanoma)

                x=BCG

 time n.risk n.event survival std.err lower 95% CI upper 95% CI

  3.9     11       1    0.909  0.0867        0.754        1.000

  5.4     10       1    0.818  0.1163        0.619        1.000

  7.9      9       1    0.727  0.1343        0.506        1.000

 10.5      8       1    0.636  0.1450        0.407        0.995

 19.5      4       1    0.477  0.1755        0.232        0.981

                x=CP

 time n.risk n.event survival std.err lower 95% CI upper 95% CI

  6.9     19       1    0.947  0.0512        0.852        1.000

  7.7     18       1    0.895  0.0704        0.767        1.000

  8.0     16       1    0.839  0.0854        0.687        1.000

  8.3     13       1    0.774  0.1003        0.601        0.998

 24.4      3       1    0.516  0.2211        0.223        1.000

> fit2

Call: survfit(formula = Surv(time, status) ~ x, data = melanoma)

       n events median 0.95LCL 0.95UCL

x=BCG 11      5   19.5    10.5      NA

x=CP  19      5     NA    24.4      NA

3. 사망시점의 사분위수 추정치와 그의 신뢰구간

> quantile(fit2,probs=c(0.25,0.5,0.75),conf.int=T)

$quantile

        25   50 75

x=BCG  7.9 19.5 NA

x=CP  24.4   NA NA

$lower

       25   50   75

x=BCG 5.4 10.5 19.5

x=CP  8.0 24.4 24.4

$upper

      25 50 75

x=BCG NA NA NA

x=CP  NA NA NA

4. 데이터 시각화

> plot(fit2,xlab='time',ylab='survival function',lty=c(1,2),col=c(1,2))

> legend(5,0.2,c('cp 처리 그룹','BCG 처리 그룹'),lty=c(2,1),col=c(2,1))

> abline(h=0.5)

> abline(v=c(10.5,24.4))

5. 로그 순위 검정법(log-rank test)과 Gehan-Wilcoxon 검정법 비교

1) 로그 순위 검정법(log-rank test)

> survdiff(Surv(time,status)~x,data=melanoma)

Call:

survdiff(formula = Surv(time, status) ~ x, data = melanoma)

       N Observed Expected (O-E)^2/E (O-E)^2/V

x=BCG 11        5     3.68     0.469     0.747

x=CP  19        5     6.32     0.274     0.747

 Chisq= 0.7  on 1 degrees of freedom, p= 0.387

2) Gehan-Wilcoxon 검정법

> survdiff(Surv(time,status)~x,rho=1,data=melanoma)

Call:

survdiff(formula = Surv(time, status) ~ x, data = melanoma, rho = 1)

       N Observed Expected (O-E)^2/E (O-E)^2/V

x=BCG 11     4.31     3.07     0.500     0.929

x=CP  19     4.10     5.34     0.288     0.929

 Chisq= 0.9  on 1 degrees of freedom, p= 0.335  

두 검정법 모두 p value가 0.05보다 크기 떄문에 유의하지 않다. 즉 BCG, CP 두 그룹 간의 생존함수는 유의한 차이를 보이지 않는다.

이상 그래프에서 보듯이 두 그룹의 누적한계추정치의 그래프도 교차되지 않고 나란한 형태를 보이므로 로그-순위 검정법이 타당한 것이었음을 알 수 있다.

출처: 보건정보데이터 분석 (이태림, 이재원, 김주한, 장대흥 공저) 

비모수적 방법 non-parametric method

6.2.2 누적한계추정법 product-limit method

생명표 방법: 기간을 1년, 6개월 등 특정 단위로 나눠 구분

누적한계추정법: 매 사건이 발생할 때 마다 해당 시점을 표기

누적한계추정법(product-limit method)은 생존함수를 추정하는 대표적인 방법 중 하나로 연구자의 이름을 따서 Kaplan-Meier추정법이라고도 한다. 모든 환자의 생존시간 또는 중도절단 시간이 각각 관찰되었다고 하자. 모든 자료의 생존시간 도는 중도절단 x1,…,투을 순서대로 배열한 것을 t1<t2<…<tn이라 하고, δi = 0, 그렇지 않은 경우 δi =1로 정의한다. 즉 중도 절단 되지 않고 사건이 발생한 경우 status = 1, 중도절단된 경우 stats = 0 로 표기한다.

예] 신장이식수술을 받은 15명 환자들의 호전기간(remission duration)이 아래와 같다고 하자. 여기서 +로 표기된 것은 중도절단된 자료를 가르킨다. 각 시점에서 생존함수를 누적한계추정법을 이용하여 추정해보자

표. 신장이식 환자들의 호전기간 (단위 : 일)

1. 데이터 불러오기

> library(survival)

> setwd('c:/Rwork')

> kidney<-read.table('kidney.txt',header=T)

> kidney #status =1 사건, 0 = 절단, 절단 표시 매우 중요

time status

1   3.0      1

2   4.0      0

3   4.5      1

4   4.5      1

5   5.5      1

6   6.0      1

7   6.4      1

8   6.5      1

9   7.0      1

10  7.5      1

11  8.4      0

12 10.0      0

13 10.0      1

14 12.0      1

15 15.0      1

> attach(kidney)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

> summary(fit1)

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

time n.risk n.event survival std.err lower 95% CI upper 95% CI

3.0     15       1    0.933  0.0644       0.8153        1.000

4.5     13       2    0.790  0.1081       0.6039        1.000

5.5     11       1    0.718  0.1198       0.5177        0.996

6.0     10       1    0.646  0.1275       0.4389        0.951

6.4      9       1    0.574  0.1320       0.3660        0.901

6.5      8       1    0.503  0.1336       0.2984        0.846

7.0      7       1    0.431  0.1324       0.2358        0.787

7.5      6       1    0.359  0.1283       0.1781        0.723

10.0      4       1    0.269  0.1237       0.1094        0.663

12.0      2       1    0.135  0.1135       0.0258        0.703

15.0      1       1    0.000     NaN           NA           NA

> fit1

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

n  events  median 0.95LCL 0.95UCL

15      12       7       6      NA

Error: invalid multibyte character in parser at line 1

-> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

6.5일 에서의 생존확률은 50.3%, 7.0일 이 지나면 생존율이 43.1%로 50%이하가 된다.

3. 사망시점의 사분위수 추정치와 그의 신뢰구간, 가령 사망자가 50%되는 시점은 언제인가?

> quantile(fit1,probs=c(0.25,0.5,0.75),conf.int=T) #신뢰구간까지

$quantile

25   50   75

5.5  7.0 12.0

$lower

25  50  75

4.5 6.0 7.0

$upper

25  50  75

7.5  NA  NA

-> 생존확률이 75%인 경우, 즉 사망확률이 25%인 경우의 시점은 5.5일, 생존확률이 50%인 경우 7.0일, 생존확율이 25%인 경우의 시점은 12.0 일

4. 누적한계추정치의 95% 신뢰구간 그래프

> plot(fit1,xlab='time',ylab='Survival function',lwd=2)

> legend(0.2,0.3,c('KM estimate','95% CI'),lty=c(1,2))

-> 점선은 신뢰구간을 의미, 실선은 생존함수추정치.생존확율이 50%인 경우는 약 7.0일임을 알 수 있다.

출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)

비모수적 방법 non-parametric method

6.2.1생명표 방법 life table method

생명표 방법 개념 추가 설명

예] 협심증(angina pectoris)이 있는 2,418명의 남성들에 대한 생존자료이다. 이 생명표에서 경과기관에 따른 생존확률을 구해보자

1. 데이터 입력

> library(KMsurv)

> setwd('C:/Rwork')

> 협심증환자자료 = read.csv('협심증환자자료.csv')

> head(협심증환자자료)

  time censor freq

1  0.5      1  456

2  0.5      0    0

3  1.5      1  226

4  1.5      0   39

5  2.5      1  152

6  2.5      0   22

> attach(협심증환자자료)

2. lifetab 함수 연구

> lifetab

function (tis, ninit, nlost, nevent)

(이하 생략)

-> tis: 시간, 길이는 17, ninit = 전체 데이터, 2418명, nlost = 중도절단, nevent 사건발생 

3. nevent 자료 만들기

> which(censor==1) #사망 데이터

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31

> 집단.사망=협심증환자자료 [which(censor==1),]

> head(집단.사망)

   time censor freq

1   0.5      1  456

3   1.5      1  226

5   2.5      1  152

7   3.5      1  171

9   4.5      1  135

11  5.5      1  125

 4. nlost 자료 만들기

> which(censor==0) #절단 데이터

 [1]  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32

> 집단.절단=협심증환자자료[which(censor==0),]

> head(집단.절단)

   time censor freq

2   0.5      0    0

4   1.5      0   39

6   2.5      0   22

8   3.5      0   23

10  4.5      0   24

12  5.5      0  107

5. ninit전체데이터 자료 만들기

> 사망자수=집단.사망[,3]

> 사망자수

 [1] 456 226 152 171 135 125  83  74  51  42  43  34  18   9   6   0

> 절단자수=집단.절단[,3]

> 절단자수

 [1]   0  39  22  23  24 107 133 102  68  64  45  53  33  27  23  30

> 합=사망자수+절단자수

> 합

 [1] 456 265 174 194 159 232 216 176 119 106  88  87  51  36  29  30

> sum(합) #2418

[1] 2418

6. tis 시간 변수 자료 만들기

> 년 = floor(집단.사망$time) #floor 내림값

> 년

 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

> lt=length(년)

> lt

[1] 16

> length(년)

[1] 17

> 년[lt+1] = NA

> 년[lt+1]

[1] NA

7. 생명표 만들기 

> 생명표=lifetab(년,sum(합),절단자수,사망자수)

> 생명표

      nsubs nlost  nrisk nevent      surv        pdf     hazard     se.surv

0-1    2418     0 2418.0    456 1.0000000 0.18858561 0.20821918 0.000000000

1-2    1962    39 1942.5    226 0.8114144 0.09440394 0.12353102 0.007955134

2-3    1697    22 1686.0    152 0.7170105 0.06464151 0.09440994 0.009179397

3-4    1523    23 1511.5    171 0.6523689 0.07380423 0.11991585 0.009734736

4-5    1329    24 1317.0    135 0.5785647 0.05930618 0.10804322 0.010138361

5-6    1170   107 1116.5    125 0.5192585 0.05813463 0.11859583 0.010304216

6-7     938   133  871.5     83 0.4611239 0.04391656 0.10000000 0.010379949

7-8     722   102  671.0     74 0.4172073 0.04601094 0.11671924 0.010450930

8-9     546    68  512.0     51 0.3711964 0.03697464 0.10483042 0.010578887

9-10    427    64  395.0     42 0.3342218 0.03553750 0.11229947 0.010717477

10-11   321    45  298.5     43 0.2986843 0.04302654 0.15523466 0.010890741

11-12   233    53  206.5     34 0.2556577 0.04209376 0.17941953 0.011124244

12-13   146    33  129.5     18 0.2135639 0.02968456 0.14937759 0.011396799

13-14    95    27   81.5      9 0.1838794 0.02030570 0.11688312 0.011765989

14-15    59    23   47.5      6 0.1635737 0.02066194 0.13483146 0.012259921

15-NA    30    30   15.0      0 0.1429117         NA         NA 0.013300258

           se.pdf   se.hazard

0-1   0.007955134 0.009697769

1-2   0.005975178 0.008201472

2-3   0.005069200 0.007649121

3-4   0.005428013 0.009153696

4-5   0.004945997 0.009285301

5-6   0.005033980 0.010588867

6-7   0.004690538 0.010962697

7-8   0.005175094 0.013545211

8-9   0.005024599 0.014659017

9-10  0.005307615 0.017300846

10-11 0.006269963 0.023601647

11-12 0.006847514 0.030646128

12-13 0.006682743 0.035110295

13-14 0.006514794 0.038894448

14-15 0.008035120 0.054919485

15-NA          NA          NA

-> 5번째 열 surv는 생존 함수, 7번째 열 hazard는 위험 함수. 협심증 환자의 약 19%가 1년 이내, 약 28%가 2년 이내에 사망하는 것으로 추정된다. 따라서 협심증 환자가 2년 이상 생존할 확률은 약 72%가 된다는 것을 알 수 있다. 또한 5년 이상 생존할 확률은 약 52%이다.

#전체 환자 중 생존확율이 70% 되는 지점

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> which.min(abs(생명표$surv-0.7))

[1] 3

> 생명표[3,]

    nsubs nlost nrisk nevent      surv        pdf     hazard     se.surv

2-3  1697    22  1686    152 0.7170105 0.06464151 0.09440994 0.009179397

       se.pdf   se.hazard

2-3 0.0050692 0.007649121

-> abs:절대값, which.min 최소값, 또는 생명표에서 3번째 행을 보면 survival이 70%에 가장 가까이 접근했음을 알 수 있다. 2년 이상 생존할 확률은 72%

#전체 환자 중 생존확율이 50% 되는 지점

> which.min(abs(생명표$surv-0.5))

[1] 6

> 생명표[6,]

    nsubs nlost  nrisk nevent      surv        pdf    hazard    se.surv

5-6  1170   107 1116.5    125 0.5192585 0.05813463 0.1185958 0.01030422

        se.pdf  se.hazard

5-6 0.00503398 0.01058887

-> 5년이 지나면 surv 즉 생존할 확률이 약 52%가 됨을 알 수 있다.

8. 데이터 시각화 – 생존함수 추정치 그래프

 > plot(년[1:lt], 생명표[,5],type='s',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)

> plot(년[1:lt], 생명표[,5],type='o',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)         

9. 데이터 시각화 – 위험함수 추정치 그래프

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> mean(생명표$hazard)

[1] NA

> 생명표$hazard

 [1] 0.20821918 0.12353102 0.09440994 0.11991585 0.10804322 0.11859583

 [7] 0.10000000 0.11671924 0.10483042 0.11229947 0.15523466 0.17941953

[13] 0.14937759 0.11688312 0.13483146         NA

> mean(생명표$hazard,na.rm=T) #NA 값 제외한 평균

[1] 0.1294874

> plot(년[1:lt], 생명표[,7],type='s',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)

> abline(h=0.1294874)

> plot(년[1:lt], 생명표[,7],type='o',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)>

> abline(h=0.1294874)

 출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)

[예] 어떤 약물에 대한 체내 배출연구에서 얻은 자료이다. 연구자는 약의 형태에 따라 체내로부터 배출되는 약물의 양이 달라지는지를 알고자 한다. 그런데 배출되는 약물의 양은 약의 형태뿐만 아니라 배출된 약물을 측정한 시간과 각 개체의 항신진대사 점수에도 영향을 받을 것으로 생각한다. 이러한 경우에는 측정시간과 항신진대사 점수를 2개의 공변량으로 하여 이들을 제어한 약의 형태에 대한 효과를 공분산분석을 통해 알 수 있다.

공변량이 2개 이상인 경우에는 공변량이 하나인 경우의 모형을 그대로 확장해서 각 모수의 추정과 검정을 할 수 있다.

<공분산분석을 위한 두 가지 가정>

1) 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉, 회귀계수가 모든 약의 형태에 대해서 동일해야 하며, 교호작용이 없어야 한다.

2) 공변량x 효과가 0이 아니다. 효과가 0이라면 분산분석을 하면 된다.

2. 공분산분석에서의 검정

1) H0: β1 = 0 항신진대사의 효과가 없음

2) H0: β2 = 0 소요시간의 효과가 없음

이상 2개의 공변량 효과를 제어한 후 배출된 약물량의 모평균이 약의 형태type에 따라 차이가 있는가를 검정할 수 있는데, 이를 검정하기 위한 귀무가설은 아래와 같다.

3) H0 : α1 = α2 = ... = αI

1. library 호출 및 데이터 불러오기

> library(lsmeans)

> setwd('C:/Rwork')

> drug  = read.csv('약물배출량자료.csv')

> head(drug)

  관측번호 약형태 항신진대사점수 소요시간  약물량

1        1      1             37       61 11.3208

2        2      2             37       37 12.9151

3        3      3             45       53 18.8947

4        4      4             41       41 14.6739

5        5      5             57       41  8.6493

6        6      6             49       33  9.5238

> colname<-c('obs','type','x1','x2','y')

> colnames(drug)<-colname

> head(drug)

  obs type x1 x2       y

1   1    1 37 61 11.3208

2   2    2 37 37 12.9151

3   3    3 45 53 18.8947

4   4    4 41 41 14.6739

5   5    5 57 41  8.6493

6   6    6 49 33  9.5238

> attach(drug)

2. 회귀계수의 동일성 검정(교호작용 존재 확인)

공분산 분석을 위한 중요한 가정으로 교호작용이 존재하는지 점검한다. 교호작용이 존재하면 공분산분석이 아닌 분산분석을 실시한다.

> model1 = aov(y~factor(type) + factor(type)*x1 + factor(type)*x2 , data = drug)

> summary(model1)

                Df Sum Sq Mean Sq F value  Pr(>F)   

factor(type)     5  250.3    50.1   2.099   0.113   

x1               1  696.6   696.6  29.206 3.9e-05 ***

x2               1   54.4    54.4   2.282   0.148   

factor(type):x1  5  160.8    32.2   1.349   0.289   

factor(type):x2  5   42.3     8.5   0.355   0.872   

Residuals       18  429.3    23.9                   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘

-> 공변량이 2개 이므로 우변에 * 형태로 입력

귀무가설 H0 : 교호작용이 존재하지 않는다.

대립가설 H1 : 교호작용이 존재한다.

p-value: x1 항신진대사점수 – 0.289

x2 소요시간 – 0.872

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 두 공변량(항신진대사점수, 소요시간)과 약의 형태 사이에 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다. ( = 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일하다.)

3. 이원공분산분석 (Two-way ANOVA)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

             Df Sum Sq Mean Sq F value   Pr(>F)   

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-> 2개의 공변량을 보정하지 않았을 때의 F-값과 p값은 각각 2.216, 0.0808으로 유의수준 5%하에서 유의하지 않다.

귀무가설 H0 : α1 = α2 = ... = αI

대립가설 H1 : 최소한 하나 이상의 약은 효과가 있다.

p-value: 0.0808

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 약의 형태에 따라 배출된 약물량의 차이가 없다.

3-1. 공변량 효과 제어 시 치료법의 효과 검정

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

              Estimate Std. Error t value Pr(>|t|)   

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

x1             -0.7606     0.1375  -5.531 6.51e-06 ***

x2              0.1647     0.1061   1.552 0.131868   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared:  0.6129,    Adjusted R-squared:  0.5161

F-statistic: 6.332 on 7 and 28 DF,  p-value: 0.0001641

귀무가설 H0 : β1 = β2 =  0  공변량(항신진대사, 소요시간)의 효과가 없음

        H0 : α1 = α2 = ... = αI

대립가설 not H0

p-value : 0.0001641

의사결정: 귀무가설을 강하게 기각한다. 5%의 유의수준 하에서 매우 유의하다.

결론: 우리가 세운 모형의 자료에 적합하다는 것을 알 수 있다. 배출된 약물량의 모평균이 약 형태type의 효과와 공변량들의 효과가 없다라고 말할 수 없다.

4. 공변량 효과 제어시 치료법의 효과 검정 - 모형 제곱합 

1)제1종 제곱합(Type I SS): SS(type)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

             Df Sum Sq Mean Sq F value   Pr(>F)   

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

OR

> summary(aov(y~factor(type)+x1+x2,data=drug))

             Df Sum Sq Mean Sq F value   Pr(>F)   

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

-> 결과해석: 약의 형태type가 기여한 제1종 제공합의 p-value는 0.0808로 유의수준 5%하에서 유의하지 않다.

2)제3종 제곱합(Type III SS): SS(type | x1, x2)

> summary(aov(y~x1+x2+factor(type),data=drug))

             Df Sum Sq Mean Sq F value   Pr(>F)   

x1            1  516.5   516.5  22.867 5.03e-05 ***

x2            1   62.8    62.8   2.779   0.1067   

factor(type)  5  422.0    84.4   3.736   0.0102 * 

Residuals    28  632.5    22.6                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘

-> 결과해석: x1,2 given type, 공변량x1,x2가 기여한 상태에서 순수하게 약의 형태type가 기여한 제3종 제공합의 p-value는 0.0102로 유의수준 5%하에서 유의하다.

5. 각 수준 별 추정치를 알아보자. 

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

              Estimate Std. Error t value Pr(>|t|)   

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

x1             -0.7606     0.1375  -5.531 6.51e-06 ***

x2              0.1647     0.1061   1.552 0.131868   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared:  0.6129,    Adjusted R-squared:  0.5161

F-statistic: 6.332 on 7 and 28 DF,  p-value: 0.0001641

-> 출력결과 유의하게 나타나는 것은 3,5,번째 약의 형태의 추정치들이 첫 번째 약의 형태와의 차이를 타나내는데 p-값이 유의수준 5% 하에서 유의하다. 항신진대사 점수(x1)는 p 값이 <0.0001로 유의수준 5%하에서 매우 유의하지만 약의 배출 소요시간(x2)은 p 값이 0.131868로 유위수준 5%하에서 유의하지 않은 것으로 나타났다. 위의 출력 결과를 바탕으로 약의 형태별 공분산모형식을 쓰면 다음과 같다.

ŷ_1j = β^₀ + α^₁ + β^1x_{11j  +} β^2x_21j  = 36.637 + 0 - 0.761 x_11j + 0.165x_21j

ŷ_2j = β^₀ + α^₂ + β^1x_{12j  +} β^2x_21j = 36.637 + 0.897- 0.761 x_12j + 0.165x_22j

ŷ_3j = β^₀ + α^₃ + β^1x_{13j  +} β^2x_21j = 36.637 + 7.908 - 0.761 x_13j + 0.165x_23j