Home > Enterprise >  Bootstrapping bioequivalence results R
Bootstrapping bioequivalence results R

Time:07-16

I am trying to run a bootstrap simulation of a bioequivalence study in R with 3 pharmacokinetic parameters (AUC0-infinity [AUCIFO], AUClast [AUCLST], Cmax [CMAX]). I want to sample with replacement from the actual study dataset to increase the sample size from 25 to 100, perform the BE calculation (using the bw2x2 function from the BE package), and then repeat that 1000 to 10000 times (depending on how long it takes to run) and generate a geometric mean ratio (GMR) for each parameter with 90%CI.

Here is a reprex dataset that I generated to represent the study date (sample size of 25):

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

The output of be2x2 is a list, which can be subsetted by the column name to obtain the GMR. Here is the output of the BE results:

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

You can subset the results such that you isolate the GMR point estimate. For example if I want the point estimate for AUCIFO, I can enter be_results$AUCIFO$`90% Confidence Interval of Geometric Mean Ratio (T/R)`[2].

I am trying to use the bootstrap package to run the bootstrapping itself. So far I have the following theta function:

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])
}

However, when I try to run the following (small nboot for troubleshooting): bootstrap(be_df, 10, theta), I get the following error message: Error in UseMethod("filter") : no applicable method for 'filter' applied to an object of class "list"

I am thinking this has something to do with the fact that I am trying to use a dataframe in the bootstrap function. Any help is appreciated!

CodePudding user response:

To replicate your bootstrap function you actually don't need the boot package; just re-write it a little.

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

You could also use the usual apply(res, 1, quantile, probs=c(.025, .975)) but matrixStats::rowQuantiles is > 50% faster.


And here is the referenced quiet version of 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)
}
  • Related