12.6 생물학적 동등성 시험의 분산 분석
sasLM package는 실제 data이면서 다음과 같이 출판된 BEdata가 포함되어 있다.
Bae KS, Kang SH. Bioequivalence data analysis for the case of separate hospitalization. Transl Clin Pharmacol. 2017;25(2):93-100.일반적인 생물학적 동등성 시험(줄여서 생동시험, BE study)에서는 대상자(subject)를, RT와 TR, 두 군으로 나눈다. RT군에 배정된 대상자는 먼저 대조약(reference drug)을 투약 받고, 약물이 제거되도록 일정 기간의 휴약기(wash-out period, 대개 일주일)가 지난 후에 시험약(test drug)을 투약 받는다. 반면 TR군에 배정된 대상자는 먼저 시험약을 투약 받고, 휴약기가 지난 후에 대조약을 투약 받는다. 즉, 대상자는 특정 약물이 아니라, 투약 받는 순서(sequence)에 무작위로 배정되며, 한 대상자는 두 가지 약물을 모두 복용하게 된다. 매 투약마다 투약 직전 시점부터 (대체로) 24-72시간 동안 12-14회의 채혈이 이루어지며, 채혈한 시점의 약물 농도를 측정한다. 측정된 약물농도로부터 AUC(area under the time-concentration curve, 농도-시간 곡선 아래 면적)와 Cmax(maximal concentration, 최고 농도)를 구하게 된다. R Package에 주어진 BEdata는 이 중에 논란이 있었던 Cmax 자료만을 제시하고 있고, 한 대상자는 두 번의 Cmax 자료를 가지고 있다. 한 사람은 하나의 순서군(sequence group, RT 또는 TR)에만 속하게 되므로, 대상자(subject)는 순서군(sequence group)에 내포(nest)되었다고 한다. 이를 SAS에서는 SUBJ(SEQ)와 같이 표현하고, R에서는 SUBJ %in% SEQ 라고 표현할 수 있다. R에서는 nested effect를 표현하는 단순한 방법으로 SEQ/SUBJ 가 있다. 이는 SEQ + SEQ:SUBJ 와 같다.
이런 문법을 이용해서 GLM 함수로 처리하면 다음과 같다.
GLM(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata)$ANOVA
Response : log(CMAX)
Df Sum Sq Mean Sq F value Pr(>F)
MODEL 48 23.1924 0.48317 5.6278 4.395e-08 ***
RESIDUALS 42 3.6059 0.08585
CORRECTED TOTAL 90 26.7983
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$Fitness
Root MSE log(CMAX) Mean Coef Var R-square Adj R-sq
0.2930098 6.071036 4.826355 0.8654428 0.7116631
$`Type I`
Df Sum Sq Mean Sq F value Pr(>F)
SEQ 1 0.6454 0.64544 7.5178 0.008938 **
SEQ:SUBJ 45 22.4395 0.49866 5.8081 3.359e-08 ***
PRD 1 0.0969 0.09686 1.1281 0.294242
TRT 1 0.0106 0.01057 0.1231 0.727410
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type II`
Df Sum Sq Mean Sq F value Pr(>F)
SEQ 1 0.6440 0.64395 7.5005 0.009011 **
SEQ:SUBJ 45 22.5232 0.50052 5.8298 3.173e-08 ***
PRD 1 0.0996 0.09958 1.1599 0.287632
TRT 1 0.0106 0.01057 0.1231 0.727410
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type III`
Df Sum Sq Mean Sq F value Pr(>F)
SEQ 1 0.3368 0.33679 3.9228 0.05421 .
SEQ:SUBJ 45 22.5232 0.50052 5.8298 3.173e-08 ***
PRD 1 0.0996 0.09958 1.1599 0.28763
TRT 1 0.0106 0.01057 0.1231 0.72741
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1그런데, 위 결과에서 그냥 Type 3 SS를 읽어서는 안 된다. 그 이유는 SEQ effect에는 SUBJ가 nest되어 있고, SUBJ effect는 random effect로 간주되기 때문에 SEQ effect 변동성의 비교 기준은 MSE가 아닌 SEQ:SUBJ term이 되기 때문이다.
즉, SEQ 행의 0.05421이라는 p-value를 사용해선 안 된다.
따라서, 이를 반영한 ANOVA는 다음과 같다.
T3test(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, Error="SEQ:SUBJ")$SEQ
Error: 0.9669*SEQ:SUBJ + 0.03311*MSE
Df Sum Sq Mean Sq F value Pr(>F)
SEQ 1.000 0.3368 0.33679 0.6919 0.4099
Error 45.529 22.1626 0.48679
$`The Rest Terms`
Df Sum Sq Mean Sq F value Pr(>F)
PRD 1 0.0996 0.09958 1.1599 0.2876
TRT 1 0.0106 0.01057 0.1231 0.7274
SEQ:SUBJ 45 22.5232 0.50052 5.8298 3.173e-08 ***
RESIDUALS 42 3.6059 0.08585
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SEQ effect에 대한 p-value는 이제 0.4099로 이전 GLM 결과의 0.05421보다 크며, 이것이 더 정확한 p-value이다.
결과에서 보듯이 오차항(error term)이 순수하게 SEQ:SUBJ만으로 이루어져 있지 않고, MSE가 약 3.3% 포함되게 되었다. 이는 unbalanced data이기 때문이다. 이 비율만 보고 싶으면, EMS 함수를 사용하면 된다.
EMS(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata) SEQ PRD TRT SEQ:SUBJ
SEQ 43.90909 0.00000 0.00000 1.869318
PRD 0.00000 43.90909 0.00000 0.000000
TRT 0.00000 0.00000 43.90909 0.000000
SEQ:SUBJ 0.00000 0.00000 0.00000 1.933333위의 표에서 행의 이름은 MS(mean square, SS/Df)들의 이름이다. 열의 이름은 fixed factor인 경우에 quadratic function(Q)이며, random factor인 경우에는 분산이다.
첫 번째 행과 마지막 행의 뜻은 다음과 같다.
\(MS_{SEQ} = 43.90909 Q_{SEQ} + 1.869318 \sigma_{SEQ:SUBJ}^2\)
\(MS_{SEQ:SUBJ} = 1.933333 \sigma_{SEQ:SUBJ}^2\)
여기에서 \(Q_{SEQ} = \sum (y_i - \bar{y})^2\)로 \(y_i\)는 i-th level의 고유값(상수)이고, \(\bar{y}\)는 level 고유값들의 평균이다.
결과표에서 EMS들의 비율이 정수가 아닌 것은 unbalanced data이기 때문이다.
따라서, 위의 경우에는 그냥 MS들만의 비(ratio)로서 \(Q_{SEQ}\)의 유의성을 판정할 수 없게 된다.
분모의 분산을 보정하는 한 가지 방법으로 SAS는 1.869318/1.933333 = 0.9669, 즉, 96.69%의 \(MS_{SEQ:SUBJ}\)만 쓰고, 부족분을 MSE에서 채운다.
그런데, 이 연구가 실제로는 대상자 모집에 어려움이 있어서 같은 날에 입원하여 시험하지 못하고, 서로 다른 3일에 걸쳐 입원하였다. 이 정보는 ADM column에 표시되어 있다. 입원일은 임상시험 계획 당시나 계획서 승인 시점에서 알 수가 없는 값이며, 대상자 모집 상태에 따라서, 달라지므로 이는 random effect로 보는 것이 타당하다. 마찬가지로, 임상시험이 아닌 다른 일반적인 시험/실험에서도 실험일은 모두 random effect로 간주한다.
하나의 입원일 안에서 다시 RT군과 TR군으로 나누어지므로 입원 SEQ도 입원일에 nest되어 있다고 볼 수도 있다. (그렇지 않다고 볼 수도 있지만, 통계분석가들이 대개 이렇게 처리하기 때문에 여기에서도 일단 그것을 따른다.)
그렇다면 GLM과 T3test는 다음과 같이 바뀌게 된다.
GLM(log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT, BEdata)$ANOVA
Response : log(CMAX)
Df Sum Sq Mean Sq F value Pr(>F)
MODEL 48 23.1924 0.48317 5.6278 4.395e-08 ***
RESIDUALS 42 3.6059 0.08585
CORRECTED TOTAL 90 26.7983
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$Fitness
Root MSE log(CMAX) Mean Coef Var R-square Adj R-sq
0.2930098 6.071036 4.826355 0.8654428 0.7116631
$`Type I`
Df Sum Sq Mean Sq F value Pr(>F)
ADM 2 0.1768 0.08839 1.0296 0.365994
ADM:SEQ 3 1.4592 0.48641 5.6655 0.002369 **
ADM:SEQ:SUBJ 41 21.4489 0.52314 6.0934 2.233e-08 ***
PRD 1 0.0969 0.09686 1.1281 0.294242
TRT 1 0.0106 0.01057 0.1231 0.727410
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type II`
Df Sum Sq Mean Sq F value Pr(>F)
ADM 2 0.1757 0.08783 1.0230 0.368290
ADM:SEQ 3 1.4636 0.48785 5.6823 0.002329 **
ADM:SEQ:SUBJ 41 21.5279 0.52507 6.1158 2.111e-08 ***
PRD 1 0.0996 0.09958 1.1599 0.287632
TRT 1 0.0106 0.01057 0.1231 0.727410
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type III`
Df Sum Sq Mean Sq F value Pr(>F)
ADM 2 0.1440 0.07201 0.8387 0.43936
ADM:SEQ 3 0.6899 0.22997 2.6786 0.05916 .
ADM:SEQ:SUBJ 41 21.5279 0.52507 6.1158 2.111e-08 ***
PRD 1 0.0996 0.09958 1.1599 0.28763
TRT 1 0.0106 0.01057 0.1231 0.72741
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1T3test(log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT, BEdata, Error="ADM:SEQ")[1]$ADM
Error: 0.854*ADM:SEQ + 0.146*MSE
Df Sum Sq Mean Sq F value Pr(>F)
ADM 2.0000 0.14402 0.07201 0.3447 0.7307
Error 3.3943 0.70916 0.20892 T3test(log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT, BEdata, Error="ADM:SEQ:SUBJ")[2:3]$`ADM:SEQ`
Error: 0.9183*ADM:SEQ:SUBJ + 0.08165*MSE
Df Sum Sq Mean Sq F value Pr(>F)
ADM:SEQ 3.000 0.6899 0.22997 0.4701 0.7047
Error 42.192 20.6407 0.48921
$`The Rest Terms`
Df Sum Sq Mean Sq F value Pr(>F)
PRD 1 0.0996 0.09958 1.1599 0.2876
TRT 1 0.0106 0.01057 0.1231 0.7274
ADM:SEQ:SUBJ 41 21.5279 0.52507 6.1158 2.111e-08 ***
RESIDUALS 42 3.6059 0.08585
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1위 결과에서 교호작용항(interaction term)을 추가할 수 있으나, 이는 random effect로 간주하여야 하고, 이것의 크기가 통계적으로 유의한 경우에도 factor 내의 level별 평균을 구하는 것은 의미가 없다.
만약 위 결과로 test drug과 reference drug의 geometric mean ratio(GMR, Test/Reference)에 대한 90% 신뢰구간을 구하면 다음과 같다.
PDIFF(log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT, BEdata, "TRT", conf.level=0.9) # log scale Estimate Lower CL Upper CL Std. Error t value Df Pr(>|t|)
R - T -0.021944 -0.12712 0.083236 0.062535 -0.3509 42 0.7274위에서 Test - Reference가 아니라, Reference - Test로 나왔기 때문에 순서를 바꾸기 위해 rev=T라는 option을 추가해준다.
exp(PDIFF(log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT, BEdata, "TRT", conf.level=0.9, rev=T)[1:3])[1] 1.0221866 0.9201339 1.1355579위 출력된 숫자는 순서대로 point estimate (점 추정치), lower limit (신뢰구간 하한), upper limit (신뢰구간 상한)이다.