首页 / 技术问答 / R问答列表 / 问答详情

自举生物等效性结果 R

我试图用3个药代动力学参数(AUC0无穷[AUCIFO]、AUClast[AUCLST]、Cmax[CMAX])对R进行生物等效性研究的自举模拟.我想从实际研究数据集中进行替换采样，将样本量从25增加到AUCIFO，执行BE计算(使用BE包中的bw2x2函数)，然后重复AUCIFO0到AUCIFO00次(取决于运行所需的时间)，并 for each 具有90%置信区间的参数生成几何平均比(GMR).

以下是我生成的代表研究日期的reprex数据集(样本量为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"))

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].

我正在try 使用bootstrap包来运行 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])
}

然而，当我try 运行以下命令(小nboot用于故障排除):bootstrap(be_df, 10, theta)时，我收到以下错误消息:Error in UseMethod("filter") : no applicable method for 'filter' applied to an object of class "list"

我认为这与我试图在bootstrap函数中使用数据帧有关.感谢您的帮助！

Solution

感谢@周.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))

推荐答案

要实现replicate个 bootstrap 函数，实际上不需要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)
}

R相关问答推荐

迭代通过1个长度的字符串长字符R

使用Facet_WRAP时更改框图中线的 colored颜色

计算满足R中条件的连续列

无法定义沿边轨迹的 colored颜色 渐变(与值无关)

将多个列值转换为二进制

仅在R中的数据集开始和结束时删除所有 Select 列的具有NA的行

派生程序包｜ ；无法检索'；return()'；的正文

如何计算R glm probit中的线性预测因子？

我如何使用循环来编写冗余的Rmarkdown脚本？

以不同于绘图中元素的方式对GG图图例进行排序

在R中的数据框上使用Apply()函数时，如何保留非数字列？

层次树图的数据树

创建在文本字符串中发现两个不同关键字的实例的数据框

减少雨云面之间的间距并绘制所有统计数据点

随机将数据帧中特定列上的某些行设置为NA

排序R矩阵的行和列

重写时间间隔模糊连接以减少内存消耗

在R中使用ggraph包排列和着色圆

如何在R中添加标识连续日期的新列

R：改进实现简单模型

实用课程推荐
相关教程推荐
休闲君评测 |  TOOLFK工具网 |  chrome插件网 |  古诗文网 |  良许Linux教程网 |  Python之禅 |  红色石头 |  豪评测 |  入党申请书 | 

© 2020 Hi LearnFk 意见反馈 免责声明 关于我们 桂ICP备11002319号-3 友链&广告位+QQ: 1963612630