8 통계처리

8.1 이 장에서는

생물학적 동등성, 용량 비례성을 확인하는 통계 처리 방법을 알아보겠습니다.

library(tidyverse)
library(BE)
library(psych)

8.2 기술통계량 구하기

앞서 3장에서 구한 Theoph_nca를 갖고 기술 통계량 (평균, 표준편차, 최소값, 최대값, skewness, kurtosis)을 구해보겠습니다. psych::describe() 함수를 사용하면 간단히 구할 수 있습니다.

desc_stat_Theoph_nca <- describe(Theoph_nca) %>% 
  select(n, mean, sd, min, max, skew, kurtosis)

knitr::kable(desc_stat_Theoph_nca, digits = 2)
n mean sd min max skew kurtosis
ID* 12 6.50 3.61 1.00 12.00 0.00 -1.50
b0* 12 2.39 0.25 2.03 2.82 0.13 -1.38
CMAX* 12 8.76 1.47 6.44 11.40 0.21 -1.19
CMAXD* 12 0.03 0.00 0.02 0.04 0.21 -1.19
TMAX* 12 1.79 1.11 0.63 3.55 0.70 -1.35
TLAG* 12 0.00 0.00 0.00 0.00 NaN NaN
CLST* 12 1.40 0.72 0.86 3.28 1.57 1.14
CLSTP* 12 1.40 0.72 0.86 3.28 1.58 1.19
TLST* 12 24.20 0.25 23.70 24.65 -0.28 -0.57
LAMZHL* 12 8.18 2.12 6.29 14.30 1.90 2.97
LAMZ* 12 0.09 0.02 0.05 0.11 -0.92 0.40
LAMZLL* 12 7.49 2.40 2.03 9.38 -1.20 -0.03
LAMZUL* 12 24.20 0.25 23.70 24.65 -0.28 -0.57
LAMZNPT* 12 3.83 1.34 3.00 7.00 1.32 0.28
CORRXY* 12 -1.00 0.00 -1.00 -1.00 2.20 3.87
R2* 12 1.00 0.00 0.99 1.00 -2.20 3.87
R2ADJ* 12 1.00 0.00 0.99 1.00 -2.05 3.39
AUCLST* 12 103.81 23.65 73.78 148.92 0.56 -1.12
AUCALL* 12 103.81 23.65 73.78 148.92 0.56 -1.12
AUCIFO* 12 122.19 38.13 84.25 216.61 1.25 0.51
AUCIFOD* 12 0.38 0.12 0.26 0.68 1.25 0.51
AUCIFP* 12 122.18 38.11 84.50 216.61 1.26 0.52
AUCIFPD* 12 0.38 0.12 0.26 0.68 1.26 0.52
AUCPEO* 12 13.54 6.35 8.13 31.25 1.71 2.19
AUCPEP* 12 13.54 6.34 8.16 31.25 1.72 2.23
AUMCLST* 12 883.06 262.98 609.15 1459.07 0.92 -0.42
AUMCIFO* 12 1590.30 1006.57 928.56 4505.53 2.00 2.96
AUMCIFP* 12 1589.85 1006.06 928.49 4505.67 2.01 2.97
AUMCPEO* 12 38.72 11.10 26.50 67.62 1.29 1.10
AUMCPEP* 12 38.72 11.07 26.59 67.62 1.30 1.14
MRTEVLST* 12 8.41 0.59 7.71 9.80 0.99 0.12
MRTEVIFO* 12 12.29 2.96 9.98 20.80 1.90 2.83
MRTEVIFP* 12 12.29 2.95 9.95 20.80 1.91 2.84
VZFO* 12 31.93 6.47 22.22 43.26 0.20 -1.40
VZFP* 12 31.92 6.46 22.22 43.14 0.19 -1.41
CLFO* 12 2.81 0.68 1.48 3.80 -0.45 -0.93
CLFP* 12 2.81 0.68 1.48 3.79 -0.46 -0.93

8.3 Dose Proportionality

DP 처리.

16명의 Cmax와 AUClast가 나온 표입니다. Table 8.1

# setup ----

library(readxl)
library(tidyverse)
library(broom)

dp_data <- # Virtual data from 4 dose groups (N=16)
'Dose,Subject,Cmax,AUClast
50,101,860,2000
50,102,510,2300
50,103,620,2900
50,104,540,2400
100,201,1550,6600
100,202,1440,7400
100,203,2000,7300
100,204,1600,7000
200,301,4100,20400
200,302,2800,9500
200,303,3200,8000
200,304,2550,7070
400,401,4800,22000
400,402,5700,23000
400,403,5800,26700
400,404,5760,28884'

sad_indi_pk <- read_csv(dp_data)
knitr::kable(sad_indi_pk, caption = '16명의 C~max~, AUC~last~')
Table 8.1: 16명의 Cmax, AUClast
Dose Subject Cmax AUClast
50 101 860 2000
50 102 510 2300
50 103 620 2900
50 104 540 2400
100 201 1550 6600
100 202 1440 7400
100 203 2000 7300
100 204 1600 7000
200 301 4100 20400
200 302 2800 9500
200 303 3200 8000
200 304 2550 7070
400 401 4800 22000
400 402 5700 23000
400 403 5800 26700
400 404 5760 28884

그림을 살펴보겠습니다.

sad_indi_pk_log <- sad_indi_pk %>% mutate_all(log)

figA <- ggplot(sad_indi_pk_log, aes(x=Dose, y=Cmax)) +
  geom_smooth(method = 'lm')+
  geom_boxplot(aes(group = Dose), 
               size = 1, 
               outlier.colour = "red", 
               outlier.shape = 1, 
               outlier.size = 3) +
  theme_bw() +
  scale_x_continuous(breaks = c(50, 100, 200, 400)) +
  labs(x = 'Dose (mg)', y = expression('C'[max]*' (ng/mL)'),
       title = expression('C'[max]))
figA

figB <- ggplot(sad_indi_pk_log, aes(x=Dose, y=AUClast)) +
  geom_smooth(method = 'lm')+
  geom_boxplot(aes(group = Dose), 
               size = 1, 
               outlier.colour = "red", 
               outlier.shape = 1, 
               outlier.size = 3) +
  theme_bw() +
  scale_x_continuous(breaks = c(50, 100, 200, 400)) +
  labs(x = 'Dose (mg)', y = expression('AUC'[(0-last)]*' (ng·hr/mL)'),
       title = expression('AUC'[(0-last)]))
figB

lm() 함수를 써서 구할 수 있습니다.

calc_dp <- function(param, fit) {
  bind_cols(fit %>% summary %>% tidy %>% filter(term == 'Dose') %>% select(1, 'estimate', 'std.error'), 
            fit %>% confint(level = 0.95) %>% tidy %>% filter(.rownames == 'Dose'), 
            fit %>% summary %>% glance
            ) %>% 
    filter(term == 'Dose') %>% 
    select(-.rownames, -term) %>% 
    mutate(parameters = param) %>% 
    mutate(est = sprintf('%0.2f (%0.2f)', estimate, std.error)) %>% 
    mutate(ci = sprintf('%0.2f-%0.2f', X2.5.., X97.5..)) %>% 
    select(parameters, est, ci, r.squared, p.value)
}

fit_cmax <- lm(formula = Cmax ~ Dose, data = sad_indi_pk_log)
fit_auclast <- lm(formula = AUClast ~ Dose, data = sad_indi_pk_log)

bind_rows(calc_dp(param = 'Cmax', fit = fit_cmax),
          calc_dp(param = 'AUClast', fit = fit_auclast))
## # A tibble: 2 x 5
##   parameters est         ci        r.squared  p.value
##   <chr>      <chr>       <chr>         <dbl>    <dbl>
## 1 Cmax       1.04 (0.06) 0.90-1.18     0.949 1.80e-10
## 2 AUClast    1.07 (0.09) 0.87-1.27     0.905 1.49e- 8