Chapter 6. Analysis of Data from Longitudinal Clinical Trials

[cardiomoon]

Chapter 6. Analysis of Data from Longitudinal Clinical Trials

이번 장에서는 longitudinal clinical trial의 데이타를 R을 이용해 분석할 것이다. Longitudinal clinical trial의 데이타의 특징은 covariate가 있고 오랜 시간에 걸쳐 데이타가 측정된다는 것이다.따라서 이러한 자료의 분석은 치료 및 covariate의 효과와 함께 시간에 따른 변화를 분석하여야 한다. 

6.3.1. Diastolic Blood Pressure Data

먼저 3장, 4장에서와 같이 데이타를 불러들여 dat에 저장한다.

> dat=read.csv("table31.csv")

> #diff=DBP5-DBP1 이라는 field를 추가한다

> dat$diff=dat$DBP5-dat$DBP1

이 데이타의 요약정보는 다음과 같이 볼수 있다.

6.3.1.1.  예비분석

데이타프레임인 dat는 “wide” form으로 “long” form으로 만들기 위해 다음의 R 코드를 사용한다.

> # call reshape 40개 자료의 변환  =>  시간별로 분리된 200개의 자료

> Dat = reshape(dat, direction="long",

     varying=c("DBP1","DBP2","DBP3","DBP4","DBP5"), 

     idvar = c("Subject","TRT","Age","Sex","diff"),sep="")

> colnames(Dat) = c("Subject","TRT","Age","Sex","diff","Time","DBP")

> Dat$Time = as.factor(Dat$Time)

> # show the first 6 observations

> head(Dat)

               Subject TRT Age Sex diff Time DBP

1.A.43.F.-9.1        1   A  43   F   -9    1 114

2.A.51.M.-15.1       2   A  51   M  -15    1 116

3.A.48.F.-21.1       3   A  48   F  -21    1 119

4.A.42.F.-14.1       4   A  42   F  -14    1 115

5.A.49.M.-11.1       5   A  49   M  -11    1 116

6.A.47.M.-15.1       6   A  47   M  -15    1 117

예비분석으로 40명 환자의  DBP의 시간에 따른 변화를 그래프로 그려본다(그림 6.1). xyplot 그래프함수는 lattice 라이브러리에 있으므로 라이브러리를 사용한다는 것을 R에게 알려준다.

> # call "lattice" library and using "xyplot" 

> library(lattice)

> # call xyplot to make the plot

> print(xyplot(DBP~Time|as.factor(Subject),type="l",groups=TRT,

   strip=strip.custom(bg="white"), lty=c(1,8), lwd=2,

   layout=c(10,4), Dat))

그림6.1 DBP as a Function of Time for Each Patient

시간에 따른 혈압변동의 추세를 보다 잘 보여주기 위해 모든 환자의 자료를 치료 군 별로 모아서 그림 6.2와 같은 그래프를 그린다. 다음과 같은 R 코드를 사용한다.

> # For better presentation of trend

> print(xyplot(DBP~Time|TRT,type="l",Dat,

   groups=as.factor(Subject),

   strip=strip.custom(bg="white")))

그림 6.2  DBP as a function of Time Grouped for all patients

그림 6.2를 보면 명백한 경향이 나타난다 : 즉  평균 DBP는 치료군 A에서 치료군 B보다 더 빨리 떨어진다. 또한 혈압 감소의 정도 또한 치료군 A에서 더 크다. 그러나 떨어지는 속도와 정도는 40명 환자 개개인에 따라 다르다.

떨어지는 추세를 보다 잘 표현하기 위해 boxplot을 그려본다.(그림 6.3)

> # call function "bwplot" to make boxplots

> print(bwplot(DBP~as.factor(Time)|TRT,Dat, xlab="Time", strip=strip.custom(bg="white")))

그림 6.3 Boxplot by Treatment

그림 6.3에서도 우리는 치료군 A의 혈압 감소 속도와 정도가 더 크다는 것을 알 수 있다.  감소 경향의 정확한 분석을 위해 그림 6.1과 그림 6.2에서 보이는 각 환자 그래프의 기울기와 절편을 구한다. 첫번째 환자의 기울기와 절편을 다음과 같이 구한다.

> # 첫번째  subject의 regression plot을 그리고 기울기,절편을 구한다.

>  mod = lm(DBP~as.numeric(Time), Dat[Dat$Subject==1,]) 

>  mod

Call:

lm(formula = DBP ~ as.numeric(Time), data = Dat[Dat$Subject == 

    1, ])

Coefficients:

     (Intercept)  as.numeric(Time)  

           118.4              -2.4  

For 루프를 사용하여 40명 환자 그래프를 선형회귀분석을 통하여 기울기와 절편을 구한다. 다음의 R 코드를 사용한다.

> # 40명 반복

> num.Subj = 40 

> # initiate the intercept and slope 

> intercept = slope = numeric(num.Subj) 

> # loop-over 

> for(i in 1:num.Subj){

    # fit regression model 

        mod = lm(DBP~as.numeric(Time), Dat[Dat$Subject==i,]) 

    # extract the intercept and slope 

       intercept[i] = coef(mod)[1] 

       slope[i] = coef(mod)[2] 

  }

> # make a dataframe "dat.coef" 

> dat.coef = data.frame(Subject=dat$Subject,TRT=dat$TRT,

  Intercept = intercept, Slope=slope) 

> # print it out 

> dat.coef

   Subject TRT Intercept Slope

1        1   A     118.4  -2.4

2        2   A     121.0  -4.0

3        3   A     125.7  -5.3

4        4   A     119.6  -3.2

5        5   A     117.8  -3.0

6        6   A     121.0  -3.8

7        7   A     119.4  -4.0

8        8   A     125.1  -4.9

9        9   A     117.7  -2.5

10      10   A     120.7  -4.3

11      11   A     120.3  -3.5

12      12   A     121.2  -3.4

13      13   A     123.5  -4.1

14      14   A     124.7  -5.1

15      15   A     118.3  -3.3

16      16   A     118.2  -3.2

17      17   A     121.0  -3.6

18      18   A     124.2  -4.2

19      19   A     119.1  -3.7

20      20   A     122.4  -3.8

21      21   B     115.0  -0.6

22      22   B     117.7  -1.7

23      23   B     116.6  -1.4

24      24   B     113.9   0.1

25      25   B     117.4  -1.8

26      26   B     116.4  -1.2

27      27   B     120.1  -0.9

28      28   B     119.9  -1.3

29      29   B     116.2  -1.6

30      30   B     119.6  -1.6

31      31   B     116.3  -0.5

32      32   B     118.9  -2.1

33      33   B     122.3  -1.7

34      34   B     117.7  -1.1

35      35   B     119.9  -1.7

36      36   B     119.9  -1.7

37      37   B     117.9  -1.9

38      38   B     116.9  -0.9

39      39   B     116.3  -0.5

40      40   B     116.8  -0.6

dat.coef 데이타를 보고 절편이 120mmHg 근처에서 변동하며 기울기는 -2.5mmHg/month 근처에서 변동하는 것을 알 수 있다. 이들 환자의 기울기와 절편을 bivariate plot으로 나타내보면 그림 6.4와 같다. 이 그래프를 보면 치료군에 따른 감소의 경향을 명백하게 보여준다.

기울기는 치료군 A에서 보다 가파르며 평균 기저혈압은 절편에서 보이는 것처럼 120mmHg정도이다.  기울기와 절편에 내재되어 있는 변이는 무작위효과(random-effects)를 보여주는 것이다.

> # Make histogram for both intercept and slope 

> int.hist = hist(intercept,plot=F)

> slope.hist = hist(slope,plot=F) 

> # make layout for plotting

> top = max(c(int.hist$counts, slope.hist$counts)) 

> nf = layout(matrix(c(2,0,1,3),2,2,byrow=T), c(3,1), c(1,3),T) 

> par(mar=c(5,4,1,1))

> # plot the intercept and slope 

> plot(Slope~Intercept,las=1,dat.coef,xlab="Intercept", ylab="Slope",

   pch=as.character(TRT)) 

> par(mar=c(0,4,1,1)) 

> # add the intercept histogram 

> barplot(int.hist$counts, axes=FALSE, ylim=c(0, top), space=0) 

> par(mar=c(5,0,1,1))

> # add the slope histogram 

> barplot(slope.hist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE)

그림 6.4.Bivariate Plot for Intercept and Slope

이제 기울기와 절편의 관계를 선형회귀를 통해 구해본다. 상호작용이 있는 것과 없는 것 두가지 모델을 사용해본다.

> # Slope and Intercept relationship by a linear regression

> # fit model 1 with interaction

> mod1.coef=lm(Slope~Intercept*TRT,dat.coef)

> summary(mod1.coef)

Call:

lm(formula = Slope ~ Intercept * TRT, data = dat.coef)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.66359 -0.29475 -0.03143  0.34701  0.75317 

Coefficients:

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

(Intercept)    28.33737    4.68741   6.045 6.04e-07 ***

Intercept      -0.26539    0.03874  -6.850 5.17e-08 ***

TRTB           -8.29639    7.37956  -1.124    0.268    

Intercept:TRTB  0.08475    0.06198   1.367    0.180    

---

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

Residual standard error: 0.4293 on 36 degrees of freedom

Multiple R-squared: 0.919, Adjusted R-squared: 0.9122 

F-statistic: 136.1 on 3 and 36 DF,  p-value: < 2.2e-16 

> # fit model 2 without interaction

> mod2.coef=lm(Slope~Intercept+TRT,dat.coef)

> summary(mod2.coef)

Call:

lm(formula = Slope ~ Intercept + TRT, data = dat.coef)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.73316 -0.38494  0.02806  0.33483  0.87272 

Coefficients:

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

(Intercept) 24.33216    3.70220   6.572 1.06e-07 ***

Intercept   -0.23228    0.03059  -7.592 4.68e-09 ***

TRTB         1.79136    0.16831  10.643 8.16e-13 ***

---

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

Residual standard error: 0.4343 on 37 degrees of freedom

Multiple R-squared: 0.9147, Adjusted R-squared: 0.9101 

F-statistic: 198.5 on 2 and 37 DF,  p-value: < 2.2e-16 

모델 1에서 보면 상호작용(Intercept:TRTB) 은 통계적으로 유의하지 않다. 모델 2에서 보면 두 치료군간에(TRTB) 유의한 차가 있다. 두 치료군 간의 차이를 t-test를 이용하여 더 분석해보면 

> t.test(Slope~TRT,dat.coef)

Welch Two Sample t-test

data:  Slope by TRT 

t = -11.673, df = 35.556, p-value = 1.019e-13

alternative hypothesis: true difference in means is not equal to 0 

95 percent confidence interval:

 -2.96976 -2.09024 

sample estimates:

mean in group A mean in group B 

         -3.765          -1.235 

> t.test(Intercept~TRT,dat.coef)

Welch Two Sample t-test

data:  Intercept by TRT 

t = 4.3669, df = 36.266, p-value = 0.0001008

alternative hypothesis: true difference in means is not equal to 0 

95 percent confidence interval:

 1.703495 4.656505 

sample estimates:

mean in group A mean in group B 

        120.965         117.785 

Longitudinal clincial trial 데이타 분석에서 이 장에서 본 것 같은 예비분석은 반응특징분석(response feature analysis)이라고 한다. 이 분석은 각각의 환자로부터 단순하면서도 예비적인 자료의 탐색과 요약을 통하여 기본적인 특징을 추출해준다. 이러한 특징 문석은 단순한 결론을 내리기 위한 기본적인 정보의 요약과 함께 추가적인 분석의 방향 또한 제시해준다. 그러나 다른 특징이 포함된다면 정보를 잃어버릴 수 있다. 보다 효율적인 분석은 자료의 모든 정보를 종합적으로 사용하는 것인데 다음 절에서 다룬다.

6.3.1.2 Longitudinal Modeling

이 분석에서는 가장 잘 맞는 모델을 찾기 위해 최근에 update된 R package인 lme4를 사용하여 일련의 모델을 적용해본다. 먼저 “Model 1”은 TRT-by-Time 상호작용을 기울기와 절편에 대한 무작위효과로 분석해보고 “Model 2”에서는 절편에 대한 효과만을 적용하여 비교한다.

 > # 6.3.1.2  Loongitudinal Modeling

> # load the library lme4 

> install.packages("lme4")  # 처음에 한번만 실행한다

> library(lme4)

> # Fit Model 1 

> mod1DBP = lmer(DBP~TRT*as.numeric(Time)+(as.numeric(Time)|Subject), Dat)  

> # Fit Model 2 

> mod2DBP = lmer(DBP~TRT*as.numeric(Time)+(1|Subject), Dat)

> # model comparison 

> anova(mod1DBP, mod2DBP)

Data: Dat

Models:

mod2DBP: DBP ~ TRT * as.numeric(Time) + (1 | Subject)

mod1DBP: DBP ~ TRT * as.numeric(Time) + (as.numeric(Time) | Subject)

        Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)

mod2DBP  6 884.57 904.36 -436.28                         

mod1DBP  8 886.49 912.87 -435.24 2.0808      2     0.3533

모델 1은 “ TRT”와 “Time”의 상호작용이 기울기와 절편에 무작위 효과로 작용하는 것이고 모델 2는 기울기에 대한 무작위효과는 없다고 보고 절편만 무작위효과를 넣은 것이다. Anova 테스트 결과 p 값은 0.35이며 이것은 두 모델이 통계적으로 다르지 않다는 것을 의미한다. 따라서 보다 단순한 모델2가 권장된다. 이 모델을 요약하면 다음과 같다.

> summary(mod2DBP)

Linear mixed model fit by REML 

Formula: DBP ~ TRT * as.numeric(Time) + (1 | Subject) 

   Data: Dat 

   AIC   BIC logLik deviance REMLdev

 889.6 909.4 -438.8    872.6   877.6

Random effects:

 Groups   Name        Variance Std.Dev.

 Subject  (Intercept) 1.3774   1.1736  

 Residual             3.8146   1.9531  

Number of obs: 200, groups: Subject, 40

Fixed effects:

                      Estimate Std. Error t value

(Intercept)           120.9650     0.5279  229.15

TRTB                   -3.1800     0.7465   -4.26

as.numeric(Time)       -3.7650     0.1381  -27.26

TRTB:as.numeric(Time)   2.5300     0.1953   12.95

Correlation of Fixed Effects:

            (Intr) TRTB   as.(T)

TRTB        -0.707              

as.nmrc(Tm) -0.785  0.555       

TRTB:s.n(T)  0.555 -0.785 -0.707

 TRT와 Time의 상호작용을 다른 두 개의 모델(모델 3 과 모델 4)을 더 사용하여 분석해본다.

> # fit Model 3 

> mod3DBP = lmer(DBP~TRT+as.numeric(Time)+(as.numeric(Time)|Subject), Dat) 

> # fit Model 4 

> mod4DBP = lmer(DBP~TRT+as.numeric(Time)+(1|Subject), Dat) 

> # model comparison 

> anova(mod3DBP, mod4DBP)

Data: Dat

Models:

mod4DBP: DBP ~ TRT + as.numeric(Time) + (1 | Subject)

mod3DBP: DBP ~ TRT + as.numeric(Time) + (as.numeric(Time) | Subject)

        Df    AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    

mod4DBP  5 999.55 1016.05 -494.78                             

mod3DBP  7 945.38  968.47 -465.69 58.172      2  2.334e-13 ***

---

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

모델 3에서는 모델2에서 TRT와 Time의 상호작용을 제거한 대신 기울기와 절편에 대한 무작위효과를 유지하기 위해 상호작용을 무작위효과로 바꾼 것이다. 모델 4는 모델 3을 보다 단순화 한 것으로 절편에 대한 무작위효과를 본 것이다. 두 모델의 비교에서 p 값이 낮으므로 모델 3이 모델 4와 통계적으로 유의하게 다르다는 것을 알 수 있다. 따라서 covariate의 효과를 보기위해 모델 3을 사용한다. 모델 3에 대한 요약은 다음과 같다.

  > summary(mod3DBP)

Linear mixed model fit by REML 

Formula: DBP ~ TRT + as.numeric(Time) + (as.numeric(Time) | Subject) 

   Data: Dat 

   AIC   BIC logLik deviance REMLdev

 947.4 970.5 -466.7    931.4   933.4

Random effects:

 Groups   Name             Variance Std.Dev. Corr   

 Subject  (Intercept)      12.6154  3.5518          

          as.numeric(Time)  1.7455  1.3212   -0.947 

 Residual                   3.5350  1.8802          

Number of obs: 200, groups: Subject, 40

Fixed effects:

                 Estimate Std. Error t value

(Intercept)      117.6482     0.6815  172.63

TRTB               3.4535     0.4553    7.58

as.numeric(Time)  -2.5000     0.2291  -10.91

Correlation of Fixed Effects:

            (Intr) TRTB  

TRTB        -0.334       

as.nmrc(Tm) -0.881  0.000

이와 같은 분석을 통해, 최종 모델 3은 다음과 같이 요약된다.

 DBP_ijk = ( β₀ + γ_ok ) + ( β₁ + γ_1k) X Time_j + TRT_i+ ε_ijk 

i는 치료군 A 또는 B, j는 1부터 5 까지의 시간, k는 환자(1부터 40) 을 의미한다. 고정 효과는  β₀  = 117.6이고 β₁ = -2.5이며 이는 DBP가 한달에 2.5 mmHg의 속도로 감소한다는 것을 뜻한다. 치료군 B와 치료군 A의 감소 속도의 차이는 3.45mmHg/month 이며 모든 인수가 통계적으로 유의하다. 무작위 효과는 다음과 같다.

모델의 진단은 모델에 맞춘후 잔차를 이용해 할 수 있다. 무작위효과의 정규성 검정을 위해 QQ plot을 다음과 같이 사용한다.

> # call "qqmath" from "lattice" for residual plot 

> print(qqmath(~resid(mod3DBP)|TRT,Dat,

   strip=strip.custom(bg="white"), xlab="Theoretical Normal", ylab="Residuals"))

그림 6.5 QQ Plot for Model 3

그림 6.5에서 보는 것처럼 QQ plot이 각 군별로 직선을 형성하고 있어 정규성 가정에 문제가 없는 것으로 보여진다.

또한 잔차를 치료군별로 시간에 대한 함수로써 나타내면 그림 6.6과 같으며 이 역시 Model 3이 합당하다는 것을 보여준다. 

 # Boxplot of residuals for all subjects by treatment

> print(bwplot(resid(mod3DBP)~as.factor(Time)|TRT,Dat, strip=strip.custom(bg="white"),   

   xlab="Time",ylim=c(-5,5), ylab="Residuals"))

그림 6.6 Boxplot of Residuals for Alll Subjects by Treatment

Covariate인 나이와 성별의 효과를 보기 위해 다른 모든 효과와 함께 모델 5와 6을 사용하여 검정해본다. 모델 5에서는 모델 3에 나이에 따른 효과를 추가하였다.

> # fit Model 3 include ``Age" effect 

> mod5DBP = lmer(DBP~TRT+as.numeric(Time)+Age+(as.numeric(Time)|Subject), Dat)

> # call anova to test ``Age" effect 

> anova(mod3DBP, mod5DBP)

Data: Dat

Models:

mod3DBP: DBP ~ TRT + as.numeric(Time) + (as.numeric(Time) | Subject)

mod5DBP: DBP ~ TRT + as.numeric(Time) + Age + (as.numeric(Time) | Subject)

        Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)    

mod3DBP  7 945.38 968.47 -465.69                             

mod5DBP  8 911.03 937.41 -447.51 36.354      1  1.645e-09 ***

---

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

anova검정에서 p 값은 1.6X10^-9으로 통계적으로 유의한 결과이다. 성별에 의한 효과는 다음과 같이 검정한다.

> # fit Model 6 including ``Age" and ``Sex" 

> mod6DBP = lmer(DBP~TRT+as.numeric(Time)+Age+Sex+(as.numeric(Time)|Subject), Dat)

> # test the ``Sex" effect 

> anova(mod5DBP, mod6DBP)

Data: Dat

Models:

mod5DBP: DBP ~ TRT + as.numeric(Time) + Age + (as.numeric(Time) | Subject)

mod6DBP: DBP ~ TRT + as.numeric(Time) + Age + Sex + (as.numeric(Time) | 

mod6DBP:     Subject)

        Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)

mod5DBP  8 911.03 937.41 -447.51                         

mod6DBP  9 912.38 942.06 -447.19 0.6507      1     0.4199

이 경우 p 값은 0.41로써 성별은 유의한 효과가 없는 것으로 나타났다. 이 결과는 제 4장의 결과와 같다.

[박윤정]

저는 선생님께서 쓰신 책을 사려고 합니다. 아침 강의 듣고 바로 가입했습니당. 대단하십니다. 그 바쁜 와중에

[cardiomoon]

헉... 벌써.... ㅅ..ㅅ ...강의때 졸지 않는 분 거의 한손에 꼽을듯했는데요....ㅋㅋㅋ