Bootstrapping生物等效性结果R

Question

我试图用3个药动学参数(AUC0-无穷大AUCIFO，AUClast AUCLST，Cmax CMAX)对R中的生物等效性进行自举模拟。我希望从实际研究数据集中替换样本，将样本大小从25增加到100，执行BE计算(使用来自bw2x2包的BE函数)，然后重复1000至10000次(取决于运行所需的时间)，并为每个参数生成一个几何平均比(GMR)，其CI为90%。

下面是我为表示研究日期(样本大小为25)而生成的reprex数据集：

require(dplyr)
require(BE)

set.seed(123)

aucifo_test <- rnorm(25, 2500, 25)
aucifo_ref <- rnorm(25, 2600, 25)
auclst_test <- rnorm(25, 2400, 25)
auclst_ref <- rnorm(25, 2400, 25)
cmax_test <- rnorm(25, 50, 5)
cmax_ref <- rnorm(25, 55, 5)

SUBJ = c(1:25)
GRP = sample(c(1,2), 25, replace = TRUE) 

be_df <- rbind(data.frame(SUBJ = SUBJ, GRP = GRP, TRT = "R", AUCIFO = aucifo_ref, AUCLST = auclst_ref, CMAX = cmax_ref),
            data.frame(SUBJ = SUBJ, GRP = GRP, TRT = "T",AUCIFO = aucifo_test, AUCLST = auclst_test, CMAX = cmax_test)) %>%
  mutate(PRD = case_when(
    GRP == 1 & TRT == "R" ~ 1,
    GRP == 1 & TRT == "T" ~ 2,
    GRP == 2 & TRT == "R" ~ 2,
    GRP == 2 & TRT == "T" ~ 1
  )) %>%
  mutate(GRP = case_when(
    GRP == 1 ~ "RT",
    GRP == 2 ~ "TR"
  )) %>%
  select(GRP, PRD, SUBJ, TRT, AUCIFO, AUCLST, CMAX)

be_results <- be2x2(be_df, c("AUCIFO", "AUCLST", "CMAX"))

be2x2的输出是一个列表，它可以通过列名进行子设置以获得GMR。以下是BE结果的输出：

print(be_results)

$AUCIFO
$AUCIFO$`Analysis of Variance (log scale)`
                  Sum Sq Df   Mean Sq  F value    Pr(>F)    
SUBJECT        0.0024154 24 0.0001006   1.8132   0.07915 .  
GROUP          0.0001999  1 0.0001999   2.0749   0.16322    
SUBJECT(GROUP) 0.0022155 23 0.0000963   1.7355   0.09686 .  
PERIOD         0.0003246  1 0.0003246   5.8476   0.02392 *  
DRUG           0.0203062  1 0.0203062 365.8577 1.332e-15 ***
ERROR          0.0012766 23 0.0000555                       
TOTAL          0.0245614 49                                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$AUCIFO$`Between and Within Subject Variability`
                                Between Subject Within Subject
Variance Estimate                  2.041146e-05    5.55029e-05
Coefficient of Variation, CV(%)    4.517928e-01    7.45013e-01

$AUCIFO$`Least Square Means (geometric mean)`
                Reference Drug Test Drug
Geometric Means       2601.982  2499.114

$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`
                 Lower Limit Point Estimate Upper Limit
90% CI for Ratio   0.9570003      0.9604654   0.9639432

$AUCIFO$`Sample Size`
                      True Ratio=1 True Ratio=Point Estimate
80% Power Sample Size            2                         2


$AUCLST
$AUCLST$`Analysis of Variance (log scale)`
                  Sum Sq Df    Mean Sq F value  Pr(>F)  
SUBJECT        0.0019802 24 0.00008251  0.9748 0.52554  
GROUP          0.0000025  1 0.00000248  0.0288 0.86668  
SUBJECT(GROUP) 0.0019778 23 0.00008599  1.0159 0.48504  
PERIOD         0.0003203  1 0.00032025  3.7837 0.06408 .
DRUG           0.0001160  1 0.00011600  1.3705 0.25371  
ERROR          0.0019467 23 0.00008464                  
TOTAL          0.0043485 49                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$AUCLST$`Between and Within Subject Variability`
                                Between Subject Within Subject
Variance Estimate                  6.746561e-07   8.464061e-05
Coefficient of Variation, CV(%)    8.213747e-02   9.200228e-01

$AUCLST$`Least Square Means (geometric mean)`
                Reference Drug Test Drug
Geometric Means       2407.201  2399.873

$AUCLST$`90% Confidence Interval of Geometric Mean Ratio (T/R)`
                 Lower Limit Point Estimate Upper Limit
90% CI for Ratio    0.992516      0.9969559    1.001416

$AUCLST$`Sample Size`
                      True Ratio=1 True Ratio=Point Estimate
80% Power Sample Size            2                         2


$CMAX
$CMAX$`Analysis of Variance (log scale)`
                Sum Sq Df  Mean Sq F value   Pr(>F)   
SUBJECT        0.17640 24 0.007350  0.6666 0.834823   
GROUP          0.00032  1 0.000321  0.0419 0.839670   
SUBJECT(GROUP) 0.17608 23 0.007656  0.6943 0.805917   
PERIOD         0.00308  1 0.003078  0.2792 0.602300   
DRUG           0.12204  1 0.122036 11.0680 0.002934 **
ERROR          0.25360 23 0.011026                    
TOTAL          0.55687 49                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$CMAX$`Between and Within Subject Variability`
                                Between Subject Within Subject
Variance Estimate                             0     0.01102599
Coefficient of Variation, CV(%)               0    10.52948160

$CMAX$`Least Square Means (geometric mean)`
                Reference Drug Test Drug
Geometric Means       53.51791  48.47896

$CMAX$`90% Confidence Interval of Geometric Mean Ratio (T/R)`
                 Lower Limit Point Estimate Upper Limit
90% CI for Ratio   0.8608553      0.9058456   0.9531872

$CMAX$`Sample Size`
                      True Ratio=1 True Ratio=Point Estimate
80% Power Sample Size            3                         6

您可以对结果进行子集，以便隔离GMR点估计。例如，如果我想要AUCIFO的点估计，我可以输入be_results$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2]。

我正在尝试使用bootstrap包来运行引导本身。到目前为止，我有以下theta函数：

theta <- function(x) {
  
  temp <- data.frame(SAMP = sample(c(1:25), 100, replace=TRUE), ID = c(1:100)) #get random sample of 100 subject IDs
  
  ref <- temp %>% left_join(., x %>% filter(TRT == "R"), by = c("SAMP" = "SUBJ")) %>% select(-SAMP) #join reference samples subject IDs with study data and replace ID with new value
  
  test <- temp %>% left_join(., x %>% filter(TRT == "T"), by = c("SAMP" = "SUBJ")) %>% select(-SAMP) #join test samples subject IDs with study data and replace ID with new value
  
  be_df <- rbind(ref,test) %>% rename(SUBJ=ID) #join test and reference into single dataset
  
  be_results <- be2x2(be_df, c("AUCIFO", "AUCLST", "CMAX"))
  AUCIFO_GMR <- c(be_results$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  AUCLST_GMR <- c(be_results$AUCLST$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  CMAX_GMR <- c(be_results$CMAX$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
}

但是，当我尝试运行以下代码(用于故障排除的小nboot )：bootstrap(be_df, 10, theta)时，我得到以下错误消息：Error in UseMethod("filter") : no applicable method for 'filter' applied to an object of class "list"

我认为这与我试图在bootstrap函数中使用数据fact这一事实有关。任何帮助都是非常感谢的！

溶液

感谢@jay.sf给出了有用的答案！对代码进行了一些额外的调整以输出数据，如下所示：

theta <- function() {
  
  temp <- data.frame(SAMP = sample(c(1:25), 100, replace=TRUE), ID = c(1:100)) #get random sample of 100 subject IDs
  
  ref <- temp %>% left_join(., x101_boot_be %>% filter(TRT == "R"), by = c("SAMP" = "SUBJ")) %>% select(-SAMP) #join reference samples subject IDs with study data and replace ID with new value
  
  test <- temp %>% left_join(., x101_boot_be %>% filter(TRT == "T"), by = c("SAMP" = "SUBJ")) %>% select(-SAMP) #join test samples subject IDs with study data and replace ID with new value
  
  be_df <- rbind(ref,test) %>% rename(SUBJ=ID) #join test and reference into single dataset
  
  be_results <- be2x2_quiet(be_df, c("AUCIFO", "AUCLST", "CMAX"))
  AUCIFO_GMR <- c(be_results$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  AUCLST_GMR <- c(be_results$AUCLST$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  CMAX_GMR <- c(be_results$CMAX$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  
  output = list(c(AUCIFO_GMR=AUCIFO_GMR, AUCLST_GMR=AUCLST_GMR, CMAX_GMR=CMAX_GMR))
  
  }

boostrap_results <- replicate(1000L, theta())

boostrap_results_df <- as.data.frame(do.call(rbind, boostrap_results))

Answer 1

要使用replicate引导函数，实际上不需要boot包；只需稍微重写一下即可。

theta1 <- function() {
  temp <- data.frame(SAMP= sample(c(1:25), 100, replace=TRUE), ID= c(1:100)) #get random sample of 100 subject IDs
  ref <- temp %>% left_join(., be_df %>% filter(TRT== "R"), by= c("SAMP"= "SUBJ")) %>% select(-SAMP) #join reference samples subject IDs with study data and replace ID with new value
  test <- temp %>% left_join(., be_df %>% filter(TRT== "T"), by= c("SAMP"= "SUBJ")) %>% select(-SAMP) #join test samples subject IDs with study data and replace ID with new value
  be_df <- rbind(ref, test) %>% rename(SUBJ=ID) #join test and reference into single dataset
  be_results <- be2x2_quiet(be_df, c("AUCIFO", "AUCLST", "CMAX"))
  AUCIFO_GMR <- c(be_results$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  AUCLST_GMR <- c(be_results$AUCLST$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  CMAX_GMR <- c(be_results$CMAX$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2])
  c(AUCIFO_GMR, AUCLST_GMR, CMAX_GMR)
}

set.seed(42)
res <- replicate(299L, theta1())

matrixStats::rowQuantiles(res, probs=c(.025, .975))
#           2.5%     97.5%
# [1,] 0.9598403 0.9632898
# [2,] 0.9984093 1.0012312
# [3,] 0.8945873 0.9559859

您也可以使用通常的apply(res, 1, quantile, probs=c(.025, .975))，但matrixStats::rowQuantiles的速度超过50%。

下面是参考的安静版本的b2x2

be2x2_quiet <-  function (Data, Columns = c("AUClast", "Cmax", "Tmax"), rtfName = "", plot=FALSE) {
  if ("data.frame" %in% class(Data)) {
    bedata = Data
  }
  else if ("character" %in% class(Data)) {
    bedata = read.csv(Data)
  }
  else {
    stop("Data should be data.frame or file name!")
  }
  bedata = bedata[order(bedata$GRP, bedata$PRD, bedata$SUBJ), 
  ]
  if (!assert(bedata)) {
    cat("\n Subject count should be balanced!\n")
    return(NULL)
  }
  nCol = length(Columns)
  if (nCol == 0) 
    stop("Input Error. Please, check the arguments!")
  if (rtfName != "") {
    rtf = RTF(rtfName)
    addHeader(rtf, title = "Bioequivalence Test Result")
    addNewLine(rtf)
    addHeader(rtf, "Table of Contents")
    addTOC(rtf)
  }
  Result = vector()
  for (i in 1:nCol) {
    if (plot) {                    ## defaults to FALSE #######################
      plot2x2(bedata, Columns[i])
    }
    if (toupper(Columns[i]) != "TMAX") {
      cResult = test2x2(bedata, Columns[i])
    }
    else {
      cResult = hodges(bedata, Columns[i])
    }
    if (rtfName != "") {
      addPageBreak(rtf)
      addHeader(rtf, title = Columns[i], TOC.level = 1)
      LineResult = capture.output(print(cResult))
      for (j in 1:length(LineResult)) addParagraph(rtf, 
                                                   LineResult[j])
      addPageBreak(rtf)
      addPlot(rtf, plot.fun = plot2x2a, width = 6.5, height = 6.5, 
              res = 300, bedata = bedata, Var = Columns[i])
      addPageBreak(rtf)
      addPlot(rtf, plot.fun = plot2x2b, width = 6.5, height = 6.5, 
              res = 300, bedata = bedata, Var = Columns[i])
    }
    Result = c(Result, list(cResult))
  }
  if (rtfName != "") {
    addPageBreak(rtf)
    addSessionInfo(rtf)
    done(rtf)
  }
  names(Result) = Columns
  return(Result)
}