1
2
3
4
5

library(ggplot2)  # 加载单个包

# 加载多个包的简便方法
packages <- c("dplyr", "tidyr", "readr")
invisible(lapply(packages, library, character.only = TRUE))

1
2
3
4
5
6
7
8

# 更新所有包
update.packages()

# 更新特定包
install.packages("ggplot2")  # 重新安装即可更新

# 检查哪些包需要更新
old.packages()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# 安装renv包
install.packages("renv")

# 初始化项目的包环境
renv::init()

# 保存当前包环境的快照
renv::snapshot()

# 恢复之前保存的包环境
renv::restore()

1
2
3
4
5
6

# 查看当前工作目录
getwd()

# 设置工作目录
setwd("/path/to/your/directory")  # Unix/Mac
setwd("C:/path/to/your/directory")  # Windows

1
2
3
4
5

# 从CRAN安装
install.packages("nlmixr2")

# 或从GitHub安装开发版本
devtools::install_github("nlmixrdevelopment/nlmixr2")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

library(nlmixr2)

# 定义模型
one.cmt <- function() {
  ini({
    tka <- 0.45 # log Ka
    tcl <- 1    # log CL
    tv <- 3.45  # log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.7
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.err)
  })
}

# 拟合模型
fit <- nlmixr(one.cmt, data=pkdata, est="saem")

# 查看结果
summary(fit)

1
2
3
4

# 从LAPKB安装
install.packages("Pmetrics", 
                 repos=c("https://repo.lapkb.com", 
                         "https://cloud.r-project.org"))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

library(Pmetrics)

# 准备数据
data <- PMreadMatrix("pkdata.csv")

# 定义模型
model <- PMbuild("model.txt", data)

# 运行分析
results <- PMrun(model)

# 查看结果
PMplot(results)

1

install.packages("mrgsolve")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

library(mrgsolve)

# 定义模型代码
code <- '
$PARAM CL = 1, VC = 20, KA = 1.2
$CMT GUT CENT
$ODE
dxdt_GUT = -KA * GUT;
dxdt_CENT = KA * GUT - (CL/VC) * CENT;
$TABLE
CP = CENT/VC;
$CAPTURE CP
'

# 编译模型
mod <- mcode("pk1", code)

# 设置剂量
e <- ev(amt=100, cmt=1)  # 100mg口服给药

# 运行模拟
out <- mod %>% ev(e) %>% mrgsim(end=24, delta=0.1)

# 可视化结果
plot(out, CP ~ time)

1

install.packages("RxODE")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

library(RxODE)

# 定义模型
ode <- "
d/dt(depot) = -ka * depot;
d/dt(cent) = ka * depot - (cl/v) * cent;
cp = cent/v;
"

# 编译模型
mod <- RxODE(model = ode)

# 设置参数
theta <- c(ka = 1.2, cl = 1, v = 20)

# 设置初始条件和时间点
inits <- c(depot = 100, cent = 0)
times <- seq(0, 24, by = 0.1)

# 运行模拟
sim <- mod$solve(theta, inits, times)

# 可视化结果
plot(sim$time, sim$cp, type = "l", 
     xlab = "Time (h)", ylab = "Concentration")

1

install.packages("ggplot2")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

library(ggplot2)

# 创建浓度-时间曲线
ggplot(pkdata, aes(x = time, y = conc, color = factor(dose))) +
  geom_point() +
  geom_line() +
  scale_x_continuous(name = "Time (h)") +
  scale_y_continuous(name = "Concentration (ng/mL)") +
  scale_color_discrete(name = "Dose (mg)") +
  theme_bw()

1
2
3

install.packages("tidyverse")  # 安装整个tidyverse
# 或者只安装dplyr
install.packages("dplyr")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

library(dplyr)

# 数据处理示例
pkdata_processed <- pkdata %>%
  filter(time > 0) %>%  # 过滤时间大于0的数据
  mutate(log_conc = log(conc)) %>%  # 创建对数浓度变量
  group_by(subject, dose) %>%  # 按受试者和剂量分组
  summarize(auc = sum(conc * diff(c(0, time))),  # 计算AUC
            cmax = max(conc),  # 计算Cmax
            tmax = time[which.max(conc)])  # 计算Tmax

1
2
3
4
5

# 安装xpose4（基于lattice）
install.packages("xpose4")

# 安装xpose（基于ggplot2）
install.packages("xpose")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

library(xpose)

# 导入nlmixr拟合结果
xpdb <- xpose_data_nlmixr(fit)

# 创建诊断图
dv_vs_pred(xpdb)
dv_vs_ipred(xpdb)
res_vs_pred(xpdb)
res_vs_idv(xpdb)
vpc(xpdb)

1

install.packages("PopED")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

library(PopED)

# 定义模型
ff <- function(model_switch, xt, parameters, poped.db) {
  with(as.list(parameters),{
    y = xt
    MS <- model_switch
    
    # 一室模型
    if(MS==1){
      y = biopharm::one.comp.oral.sd.CL(xt, CL, V, KA)
    }
    
    return(list(y=y, poped.db=poped.db))
  })
}

# 设置模型参数
sfg <- function(x, a, bpop, b, bocc) {
  parameters = c(CL = bpop[1] * exp(b[1]), 
                 V = bpop[2] * exp(b[2]), 
                 KA = bpop[3] * exp(b[3]))
  return(parameters)
}

# 创建PopED数据库
poped.db <- create.poped.database(
  ff_fun = ff,
  fg_fun = sfg,
  bpop = c(CL = 1, V = 20, KA = 1.2),
  notfixed_bpop = c(1, 1, 1),
  d = c(CL = 0.09, V = 0.09, KA = 0.36),
  sigma = 0.04,
  groupsize = 20,
  xt = c(0.5, 1, 2, 4, 8, 12, 24),
  model_switch = 1
)

# 运行优化
output <- poped_optimize(poped.db)

# 查看结果
summary(output)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# 创建数值向量
conc <- c(0.1, 0.5, 2.3, 4.7, 3.2, 1.8, 0.9)

# 创建字符向量
subject_id <- c("S001", "S002", "S003", "S004")

# 创建逻辑向量
above_threshold <- conc > 2.0

# 向量运算
conc * 2  # 所有元素乘以2
log(conc)  # 对所有元素取对数
mean(conc)  # 计算平均值
sd(conc)  # 计算标准差

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# 创建矩阵
dose_matrix <- matrix(c(10, 20, 30, 10, 20, 30), nrow = 3, ncol = 2)

# 为行和列命名
rownames(dose_matrix) <- c("Subject1", "Subject2", "Subject3")
colnames(dose_matrix) <- c("Day1", "Day2")

# 矩阵索引
dose_matrix[1, 2]  # 第1行第2列的元素
dose_matrix[, 1]   # 第1列的所有元素
dose_matrix[2, ]   # 第2行的所有元素

# 矩阵运算
t(dose_matrix)  # 转置
dose_matrix * 2  # 所有元素乘以2

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# 创建数据框
pk_data <- data.frame(
  subject = c("S001", "S001", "S001", "S002", "S002", "S002"),
  time = c(0, 1, 2, 0, 1, 2),
  conc = c(0, 10.5, 5.2, 0, 8.7, 4.3)
)

# 访问数据框元素
pk_data$subject  # 访问subject列
pk_data[, "time"]  # 另一种访问列的方式
pk_data[1, ]  # 第1行
pk_data[pk_data$time > 0, ]  # 条件筛选

# 添加新列
pk_data$log_conc <- log(pk_data$conc + 0.1)  # 添加对数浓度列

# 数据框汇总
summary(pk_data)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# 创建列表
pk_model <- list(
  parameters = c(CL = 1.0, V = 20.0, KA = 1.2),
  data = pk_data,
  description = "One-compartment model with first-order absorption"
)

# 访问列表元素
pk_model$parameters  # 访问parameters元素
pk_model[[1]]  # 另一种访问元素的方式
pk_model$data$conc  # 嵌套访问

1
2
3
4
5

# 导入CSV文件
pk_data <- read.csv("data/pk_data.csv", header = TRUE)

# 导出CSV文件
write.csv(pk_data, "output/processed_pk_data.csv", row.names = FALSE)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# 安装readxl包
install.packages("readxl")

# 导入Excel文件
library(readxl)
pk_data <- read_excel("data/pk_data.xlsx", sheet = "PK_Data")

# 导出Excel文件（使用writexl包）
install.packages("writexl")
library(writexl)
write_xlsx(pk_data, "output/processed_pk_data.xlsx")

1
2
3
4
5
6
7
8
9

# 安装nonmem2rx包
install.packages("nonmem2rx")

# 导入NONMEM数据文件
library(nonmem2rx)
nm_data <- read_nonmem("data/pk.csv")

# 导入NONMEM控制流文件
nm_model <- nonmem2rx("models/run001.mod")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 安装必要的包
install.packages(c("DBI", "RSQLite"))

# 连接SQLite数据库
library(DBI)
con <- dbConnect(RSQLite::SQLite(), "data/pk_database.sqlite")

# 查询数据
pk_data <- dbGetQuery(con, "SELECT * FROM pk_measurements")

# 关闭连接
dbDisconnect(con)

1
2
3
4
5
6
7

# 使用基础R
pk_data_subset <- pk_data[pk_data$time > 0 & pk_data$conc > 0, ]

# 使用dplyr
library(dplyr)
pk_data_subset <- pk_data %>%
  filter(time > 0, conc > 0)

1
2
3
4
5
6

# 使用基础R
pk_data_sorted <- pk_data[order(pk_data$subject, pk_data$time), ]

# 使用dplyr
pk_data_sorted <- pk_data %>%
  arrange(subject, time)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 使用基础R
pk_data$dose_per_kg <- pk_data$dose / pk_data$weight
pk_data$log_conc <- log(pk_data$conc + 0.1)  # 添加0.1避免log(0)

# 使用dplyr
pk_data <- pk_data %>%
  mutate(
    dose_per_kg = dose / weight,
    log_conc = log(conc + 0.1)
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 使用基础R
pk_summary <- aggregate(conc ~ subject + dose, data = pk_data, 
                        FUN = function(x) c(max = max(x), auc = sum(x)))

# 使用dplyr
pk_summary <- pk_data %>%
  group_by(subject, dose) %>%
  summarize(
    cmax = max(conc),
    tmax = time[which.max(conc)],
    auc = sum(conc * diff(c(0, time)))
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# 安装tidyr包
install.packages("tidyr")

# 宽格式转长格式
library(tidyr)
pd_data_long <- pd_data_wide %>%
  pivot_longer(
    cols = c(bp, hr, temp),
    names_to = "parameter",
    values_to = "value"
  )

# 长格式转宽格式
pd_data_wide <- pd_data_long %>%
  pivot_wider(
    names_from = parameter,
    values_from = value
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# 检查缺失值
sum(is.na(pk_data$conc))
colSums(is.na(pk_data))

# 删除含有缺失值的行
pk_data_complete <- na.omit(pk_data)

# 使用dplyr删除特定列中含有缺失值的行
pk_data_complete <- pk_data %>%
  filter(!is.na(conc))

# 替换缺失值
pk_data$conc[is.na(pk_data$conc)] <- 0  # 替换为0
pk_data$conc[is.na(pk_data$conc)] <- min(pk_data$conc, na.rm = TRUE) / 2  # 替换为最小值的一半

# 使用插值替换缺失值
install.packages("zoo")
library(zoo)
pk_data$conc <- na.approx(pk_data$conc, na.rm = FALSE)  # 线性插值

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# if-else语句
if (max(pk_data$conc) > 100) {
  print("High concentration detected")
} else {
  print("Concentration within normal range")
}

# ifelse函数（向量化操作）
pk_data$conc_category <- ifelse(pk_data$conc > 10, "High", "Low")

# 多条件判断
pk_data$dose_level <- case_when(
  pk_data$dose < 10 ~ "Low",
  pk_data$dose >= 10 & pk_data$dose < 50 ~ "Medium",
  pk_data$dose >= 50 ~ "High"
)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

# for循环
subjects <- unique(pk_data$subject)
for (subj in subjects) {
  subj_data <- pk_data[pk_data$subject == subj, ]
  cat("Subject", subj, "has", nrow(subj_data), "observations\n")
}

# while循环
i <- 1
while (i <= length(subjects)) {
  cat("Processing subject", subjects[i], "\n")
  i <- i + 1
}

# apply系列函数（替代循环的更高效方法）
# lapply：返回列表
subject_counts <- lapply(subjects, function(subj) {
  sum(pk_data$subject == subj)
})

# sapply：简化结果
subject_counts <- sapply(subjects, function(subj) {
  sum(pk_data$subject == subj)
})

# tapply：按因子分组应用函数
max_conc_by_subject <- tapply(pk_data$conc, pk_data$subject, max)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

# 创建计算药代动力学参数的函数
calc_pk_params <- function(time, conc, dose = NULL) {
  # 输入检查
  if (length(time) != length(conc)) {
    stop("Time and concentration vectors must have the same length")
  }
  
  # 排序数据（确保按时间顺序）
  ord <- order(time)
  time <- time[ord]
  conc <- conc[ord]
  
  # 计算参数
  cmax <- max(conc)
  tmax <- time[which.max(conc)]
  
  # 使用线性梯形法计算AUC
  auc <- sum(diff(time) * (conc[-1] + conc[-length(conc)]) / 2)
  
  # 如果提供了剂量，计算清除率
  cl <- ifelse(!is.null(dose), dose / auc, NA)
  
  # 返回结果
  return(list(
    cmax = cmax,
    tmax = tmax,
    auc = auc,
    cl = cl
  ))
}

# 使用函数
subj_data <- pk_data[pk_data$subject == "S001", ]
pk_params <- calc_pk_params(subj_data$time, subj_data$conc, subj_data$dose[1])
print(pk_params)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# try-catch结构
fit_model <- function(data) {
  tryCatch(
    {
      # 尝试拟合模型
      model <- nlmixr(one.cmt, data = data, est = "saem")
      return(model)
    },
    error = function(e) {
      # 处理错误
      message("Error fitting model: ", e$message)
      return(NULL)
    },
    warning = function(w) {
      # 处理警告
      message("Warning during model fitting: ", w$message)
      # 继续执行
    },
    finally = {
      # 无论是否发生错误，都会执行
      message("Model fitting attempt completed")
    }
  )
}

# 使用函数
model <- fit_model(pk_data)
if (!is.null(model)) {
  summary(model)
} else {
  message("Model fitting failed, using alternative approach")
  # 替代方法
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# 创建示例PK数据
pk_data <- data.frame(
  ID = rep(1:3, each = 8),
  TIME = rep(c(0, 0.5, 1, 2, 4, 8, 12, 24), 3),
  DV = c(0, 2.1, 4.5, 3.2, 2.0, 1.1, 0.6, 0.2,
         0, 1.8, 3.9, 2.8, 1.7, 0.9, 0.5, 0.1,
         0, 2.4, 5.1, 3.6, 2.3, 1.3, 0.7, 0.3),
  AMT = c(100, 0, 0, 0, 0, 0, 0, 0,
          100, 0, 0, 0, 0, 0, 0, 0,
          100, 0, 0, 0, 0, 0, 0, 0),
  CMT = 1,
  EVID = c(1, 0, 0, 0, 0, 0, 0, 0,
           1, 0, 0, 0, 0, 0, 0, 0,
           1, 0, 0, 0, 0, 0, 0, 0),
  MDV = c(1, 0, 0, 0, 0, 0, 0, 0,
          1, 0, 0, 0, 0, 0, 0, 0,
          1, 0, 0, 0, 0, 0, 0, 0),
  WT = rep(c(70, 65, 80), each = 8),
  AGE = rep(c(35, 42, 28), each = 8),
  SEX = rep(c(1, 0, 1), each = 8)
)

# 查看数据结构
head(pk_data)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

library(dplyr)
library(tidyr)

# 检查并处理缺失值
pk_data_clean <- pk_data %>%
  # 将浓度为NA的观测标记为MDV=1
  mutate(MDV = ifelse(is.na(DV), 1, MDV))

# 检查并处理异常值
pk_data_clean <- pk_data_clean %>%
  # 将浓度小于定量下限(LOQ)的值替换为LOQ/2
  mutate(DV = ifelse(DV < 0.1 & DV > 0, 0.05, DV))

# 添加对数转换的浓度
pk_data_clean <- pk_data_clean %>%
  mutate(LOGDV = ifelse(DV > 0, log(DV), NA))

# 标准化协变量
pk_data_clean <- pk_data_clean %>%
  mutate(
    WT_NORM = WT / 70,  # 体重标准化为70kg
    AGE_CENT = AGE - median(AGE)  # 年龄中心化
  )

# 为不同给药场合添加OCC变量（如果有多次给药）
pk_data_clean <- pk_data_clean %>%
  group_by(ID) %>%
  mutate(OCC = cumsum(EVID == 1)) %>%
  ungroup()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

library(ggplot2)

# 绘制浓度-时间曲线
ggplot(pk_data_clean %>% filter(EVID == 0), 
       aes(x = TIME, y = DV, color = factor(ID))) +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  labs(title = "浓度-时间曲线",
       x = "时间 (h)",
       y = "浓度 (ng/mL)",
       color = "受试者") +
  theme_bw()

# 计算并可视化主要PK参数
pk_params <- pk_data_clean %>%
  filter(EVID == 0) %>%
  group_by(ID) %>%
  summarize(
    Cmax = max(DV),
    Tmax = TIME[which.max(DV)],
    AUClast = sum(diff(c(0, TIME)) * (DV[-1] + DV[-length(DV)]) / 2),
    T_half = NA  # 半衰期需要通过模型估计
  )

# 查看PK参数
print(pk_params)

# 可视化PK参数与协变量的关系
ggplot(pk_params %>% 
         left_join(pk_data_clean %>% 
                     select(ID, WT, AGE, SEX) %>% 
                     distinct(), by = "ID"), 
       aes(x = WT, y = Cmax, color = factor(SEX))) +
  geom_point() +
  labs(title = "Cmax vs 体重",
       x = "体重 (kg)",
       y = "Cmax (ng/mL)",
       color = "性别") +
  theme_bw()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# 安装PKNCA包
install.packages("PKNCA")

# 加载包
library(PKNCA)
library(dplyr)

# 准备数据
nca_data <- pk_data_clean %>%
  filter(EVID == 0) %>%
  select(ID, TIME, DV, AMT)

# 设置NCA参数
conc_obj <- PKNCAconc(nca_data, 
                      conc ~ TIME | ID,
                      exclude = DV < 0)

dose_data <- pk_data_clean %>%
  filter(EVID == 1) %>%
  select(ID, TIME, AMT)

dose_obj <- PKNCAdose(dose_data,
                      dose ~ TIME | ID)

# 定义要计算的NCA参数
nca_parameters <- data.frame(
  parameter = c("auclast", "cmax", "tmax", "half.life", "cl.obs"),
  method = c("linear", "max", "max", "calc", "calc")
)

# 执行NCA分析
nca_obj <- PKNCAdata(conc_obj, dose_obj)
nca_results <- pk.nca(nca_obj, nca_parameters)

# 查看结果
summary(nca_results)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

# 提取NCA结果
nca_summary <- as.data.frame(nca_results$result)

# 合并协变量数据
nca_with_cov <- nca_summary %>%
  filter(parameter %in% c("auclast", "cmax", "tmax", "half.life")) %>%
  select(ID = subject, parameter, value) %>%
  pivot_wider(names_from = parameter, values_from = value) %>%
  left_join(pk_data_clean %>% 
              select(ID, WT, AGE, SEX) %>% 
              distinct(), by = "ID")

# 可视化NCA参数
# AUC vs 体重
ggplot(nca_with_cov, aes(x = WT, y = auclast, color = factor(SEX))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "AUC vs 体重",
       x = "体重 (kg)",
       y = "AUC (h*ng/mL)",
       color = "性别") +
  theme_bw()

# Cmax vs 体重
ggplot(nca_with_cov, aes(x = WT, y = cmax, color = factor(SEX))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Cmax vs 体重",
       x = "体重 (kg)",
       y = "Cmax (ng/mL)",
       color = "性别") +
  theme_bw()

# 半衰期分布
ggplot(nca_with_cov, aes(x = half.life)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  labs(title = "半衰期分布",
       x = "半衰期 (h)",
       y = "频数") +
  theme_bw()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

# 假设我们有两种制剂（测试制剂T和参比制剂R）的NCA结果
be_data <- data.frame(
  subject = rep(1:12, each = 2),
  formulation = rep(c("T", "R"), 12),
  sequence = rep(rep(c("TR", "RT"), each = 2), 6),
  period = rep(c(1, 2, 2, 1), 6),
  AUC = c(rnorm(12, mean = 100, sd = 20), rnorm(12, mean = 105, sd = 20)),
  Cmax = c(rnorm(12, mean = 15, sd = 3), rnorm(12, mean = 16, sd = 3))
)

# 计算几何平均比和90%置信区间
library(bear)

# 对数转换
be_data$logAUC <- log(be_data$AUC)
be_data$logCmax <- log(be_data$Cmax)

# 使用混合效应模型进行BE分析
library(lme4)

# AUC的BE分析
auc_model <- lmer(logAUC ~ formulation + sequence + period + (1|subject), 
                  data = be_data)
auc_summary <- summary(auc_model)

# 提取固定效应
auc_fixed <- fixef(auc_model)
auc_T_R_ratio <- exp(auc_fixed["formulationT"])  # T相对于R的比值

# 计算90%置信区间
auc_ci <- confint(auc_model, parm = "formulationT", level = 0.9)
auc_lower_ci <- exp(auc_ci[1])
auc_upper_ci <- exp(auc_ci[2])

# Cmax的BE分析
cmax_model <- lmer(logCmax ~ formulation + sequence + period + (1|subject), 
                   data = be_data)
cmax_summary <- summary(cmax_model)

# 提取固定效应
cmax_fixed <- fixef(cmax_model)
cmax_T_R_ratio <- exp(cmax_fixed["formulationT"])  # T相对于R的比值

# 计算90%置信区间
cmax_ci <- confint(cmax_model, parm = "formulationT", level = 0.9)
cmax_lower_ci <- exp(cmax_ci[1])
cmax_upper_ci <- exp(cmax_ci[2])

# 汇总BE结果
be_results <- data.frame(
  Parameter = c("AUC", "Cmax"),
  Ratio = c(auc_T_R_ratio, cmax_T_R_ratio) * 100,
  Lower_90CI = c(auc_lower_ci, cmax_lower_ci) * 100,
  Upper_90CI = c(auc_upper_ci, cmax_upper_ci) * 100,
  BE_Conclusion = c(
    ifelse(auc_lower_ci >= 0.8 & auc_upper_ci <= 1.25, "BE", "Not BE"),
    ifelse(cmax_lower_ci >= 0.8 & cmax_upper_ci <= 1.25, "BE", "Not BE")
  )
)

# 打印BE结果
print(be_results)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

# 安装并加载nlmixr2
install.packages("nlmixr2")
library(nlmixr2)

# 准备数据
# 确保数据包含ID, TIME, DV, AMT, CMT, EVID, MDV等必要列

# 定义一室模型（口服给药）
one_cmt_model <- function() {
  ini({
    # 固定效应参数（群体典型值）
    tka <- log(1.0)  # log(Ka)
    tcl <- log(5.0)  # log(CL)
    tv <- log(80.0)  # log(V)
    
    # 随机效应参数（个体间变异）
    eta.ka ~ 0.4
    eta.cl ~ 0.3
    eta.v ~ 0.2
    
    # 残差误差模型
    add.err <- 0.1  # 加性误差
    prop.err <- 0.2  # 比例误差
  })
  model({
    # 参数转换（从对数尺度）
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    
    # 定义微分方程
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    
    # 或使用内置的线性房室模型函数
    # linCmt() ~ add(add.err) + prop(prop.err)
    
    # 连接模型（预测值与观测值的关系）
    cp = central / v
    cp ~ add(add.err) + prop(prop.err)
  })
}

# 拟合模型
fit_one_cmt <- nlmixr(one_cmt_model, data = pk_data_clean, est = "saem")

# 查看拟合结果
summary(fit_one_cmt)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

# 定义二室模型（口服给药）
two_cmt_model <- function() {
  ini({
    # 固定效应参数
    tka <- log(1.0)    # log(Ka)
    tcl <- log(5.0)    # log(CL)
    tv1 <- log(50.0)   # log(V1)
    tq <- log(10.0)    # log(Q)
    tv2 <- log(100.0)  # log(V2)
    
    # 随机效应参数
    eta.ka ~ 0.4
    eta.cl ~ 0.3
    eta.v1 ~ 0.2
    eta.q ~ 0.3
    eta.v2 ~ 0.3
    
    # 残差误差模型
    add.err <- 0.1
    prop.err <- 0.2
  })
  model({
    # 参数转换
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v1 <- exp(tv1 + eta.v1)
    q <- exp(tq + eta.q)
    v2 <- exp(tv2 + eta.v2)
    
    # 定义微分方程
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v1 * central - q/v1 * central + q/v2 * peripheral
    d/dt(peripheral) = q/v1 * central - q/v2 * peripheral
    
    # 或使用内置的线性房室模型函数
    # linCmt() ~ add(add.err) + prop(prop.err)
    
    # 连接模型
    cp = central / v1
    cp ~ add(add.err) + prop(prop.err)
  })
}

# 拟合模型
fit_two_cmt <- nlmixr(two_cmt_model, data = pk_data_clean, est = "saem")

# 查看拟合结果
summary(fit_two_cmt)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70