引言

药代动力学（Pharmacokinetics, PK）和药效动力学（Pharmacodynamics, PD）建模是现代药物研发和临床应用中不可或缺的工具。非线性混合效应建模（Nonlinear Mixed-Effects Modeling, NLME）作为其中的核心方法，能够同时描述群体和个体水平的药物行为，为药物剂量优化、临床试验设计和特殊人群用药提供科学依据。

nlmixr2是一个用于非线性混合效应建模的开源R软件包，它结合了多种建模方法，支持参数估计、模型诊断和仿真等工作。作为NONMEM、Monolix和Phoenix NLME等商业软件的开源替代品，nlmixr2具有以下优势：

开源免费：作为R包，nlmixr2完全开源且免费使用，降低了研究成本
集成R生态系统：可以无缝集成R的数据处理、统计分析和可视化功能
多种估计方法：提供多种参数估计算法，包括FOCE、SAEM、FOCEI等
灵活的模型定义：支持多种模型定义方式，包括微分方程和闭式解
丰富的诊断工具：集成了全面的模型诊断和评估功能
强大的仿真能力：通过与rxode2包的集成，提供高效的模型仿真功能

这是一份Manus整理的教程，虽然绝大部分内容我都看过，但难免有一本正经的胡说八道的地方，教程的正确性需要有智慧的眼光来看待。本教程旨在为药代动力学/药效动力学研究者和临床药师提供全面的nlmixr2使用指南，从基本概念到高级应用，帮助读者掌握这一强大工具的使用方法。无论您是初学者还是有经验的建模者，本教程都将为您提供有价值的参考。

接下来，我们将从nlmixr2的安装和设置开始，逐步深入探讨非线性混合效应建模的理论基础、参数估计方法、模型诊断技术、模型仿真方法，并通过一个完整的实际案例展示nlmixr2的应用流程。

让我们开始这段非线性混合效应建模的探索之旅！

nlmixr2安装和设置指南

本指南将详细介绍如何在不同操作系统上安装和设置nlmixr2软件包，包括必要的依赖项和可能遇到的问题解决方案。

安装前的准备工作

nlmixr2是一个用于非线性混合效应建模的R软件包，在安装前需要确保以下条件：

已安装R（推荐使用最新版本，至少R 4.0.0以上）
已安装RStudio（可选，但推荐使用）
具有安装R包的权限
安装了适当的编译工具（因为nlmixr2和rxode2需要编译器支持）

不同操作系统的安装步骤

Windows系统

安装R和RStudio
- 从R官网下载并安装最新版本的R
- 从RStudio官网下载并安装RStudio
安装Rtools
- nlmixr2需要Rtools进行编译
- 从Rtools官网下载与您的R版本对应的Rtools
- 安装时确保选择"Add rtools to system PATH"选项
安装nlmixr2
- 打开RStudio，在控制台中运行以下命令：
1
install.packages("nlmixr2", dependencies = TRUE)

对于R 4.0.x和R 4.1.x版本

需要降级symengine包：

1
2
3
4
5
6
# 首先安装remotes包（如果尚未安装）
install.packages("remotes")
# 然后安装特定版本的symengine
remotes::install_version("symengine", version = "0.1.6")
# 最后安装nlmixr2
install.packages("nlmixr2", dependencies = TRUE)

macOS系统

安装R和RStudio
- 从R官网下载并安装最新版本的R
- 从RStudio官网下载并安装RStudio
安装编译工具
- 从App Store安装Xcode
- 安装gfortran：
  - 从mac.r-project.org/tools/下载并安装gfortran
  - 添加gfortran目录到系统路径，在终端中运行：
  1
  export PATH=$PATH:/usr/local/gfortran/bin
  - 为了永久添加到路径，可以将上述命令添加到~/.bash_profile或~/.zshrc文件中
安装nlmixr2
- 打开RStudio，在控制台中运行以下命令：
1
install.packages("nlmixr2", dependencies = TRUE)

对于R 4.0.x和R 4.1.x版本

需要降级symengine包：

1
2
3
4
5
6
# 首先安装remotes包（如果尚未安装）
install.packages("remotes")
# 然后安装特定版本的symengine
remotes::install_version("symengine", version = "0.1.6")
# 最后安装nlmixr2
install.packages("nlmixr2", dependencies = TRUE)

Linux系统（Ubuntu/Debian为例）

安装R和RStudio

安装R：

1
2
3
4
5
6
7
8
9
# 添加CRAN仓库
sudo apt update
sudo apt install --no-install-recommends software-properties-common dirmngr
wget -qO- https://cloud.r-project.org/bin/linux/ubuntu/marutter_pubkey.asc | sudo tee -a /etc/apt/trusted.gpg.d/cran_ubuntu_key.asc
sudo add-apt-repository "deb https://cloud.r-project.org/bin/linux/ubuntu $(lsb_release -cs)-cran40/"

# 安装R
sudo apt update
sudo apt install r-base r-base-dev

安装RStudio：

1
2
3
# 下载最新版本的RStudio（请检查最新版本链接）
wget https://download1.rstudio.org/desktop/bionic/amd64/rstudio-2023.06.0-421-amd64.deb
sudo apt install ./rstudio-2023.06.0-421-amd64.deb

安装编译工具和依赖项

1
2
sudo apt update
sudo apt install build-essential libcurl4-openssl-dev libxml2-dev libssl-dev gfortran

安装nlmixr2
- 打开RStudio，在控制台中运行以下命令：
1
install.packages("nlmixr2", dependencies = TRUE)

对于R 4.0.x和R 4.1.x版本

需要降级symengine包：

1
2
3
4
5
6
# 首先安装remotes包（如果尚未安装）
install.packages("remotes")
# 然后安装特定版本的symengine
remotes::install_version("symengine", version = "0.1.6")
# 最后安装nlmixr2
install.packages("nlmixr2", dependencies = TRUE)

验证安装

无论使用哪种操作系统，安装完成后都应该验证nlmixr2是否正确安装：

1
2
3
4
5
# 加载nlmixr2包
library(nlmixr2)

# 检查安装
nlmixr2::nlmixr2CheckInstall()

如果安装正确，应该会看到成功消息。如果出现错误，请参考下面的常见问题解决方案。

安装开发版本（可选）

如果您需要最新的开发版本，可以通过以下方式安装：

使用R universe（推荐）

这是安装开发版本最快的方式，因为它为最新和上一个版本的R提供了mac和windows的二进制文件（无需等待编译）：

1
2
3
4
5
install.packages(c("dparser", "nlmixr2data", "lotri", "rxode2ll",
                   "rxode2", "nlmixr2est", "nlmixr2extra", "nlmixr2plot",
                   "nlmixr2"),
                 repos = c('https://nlmixr2.r-universe.dev',
                           'https://cloud.r-project.org'))

使用remotes包

这确保获取最新的开发版本：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# 安装remotes包（如果尚未安装）
install.packages("remotes")

# 安装nlmixr2及其依赖项
remotes::install_github("nlmixr2/dparser-R")
remotes::install_github("nlmixr2/nlmixr2data")
remotes::install_github("nlmixr2/lotri")
remotes::install_github("nlmixr2/rxode2ll")
remotes::install_github("nlmixr2/rxode2")
remotes::install_github("nlmixr2/nlmixr2est")
remotes::install_github("nlmixr2/nlmixr2extra")
remotes::install_github("nlmixr2/nlmixr2plot")
remotes::install_github("nlmixr2/nlmixr2")

安装支持包（可选）

nlmixr2生态系统中有许多有用的支持包：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
install.packages(c("xpose.nlmixr2",    # 基于xpose的额外拟合优度图
                   "nlmixr2targets",   # 简化与targets包的工作
                   "babelmixr2",       # 在nlmixr2模型和NONMEM、Monolix之间转换
                   "nonmem2rx",        # 将NONMEM转换为rxode2/nlmixr2模型
                   "nlmixr2lib",       # 模型库和模型修改函数
                   "nlmixr2rpt"),      # nlmixr2的自动Microsoft Word和PowerPoint报告
                 repos = c('https://nlmixr2.r-universe.dev',
                           'https://cloud.r-project.org'))

# 一些额外的包
remotes::install_github("ggPMXdevelopment/ggPMX")       # 拟合优度图
remotes::install_github("RichardHooijmaijers/shinyMixR") # Shiny运行管理器

常见问题解决方案

1. 编译器问题

问题：安装过程中出现与编译器相关的错误。

解决方案：

Windows：确保已正确安装Rtools并添加到系统路径
macOS：确保已安装Xcode和gfortran，并正确设置路径
Linux：确保已安装所有必要的开发工具和库

2. 依赖包安装失败

问题：某些依赖包安装失败。

解决方案：

尝试单独安装失败的依赖包
检查是否有特定版本要求，可能需要先安装特定版本的依赖包
对于R 4.0.x和R 4.1.x，确保先降级symengine包

3. 内存不足错误

问题：安装过程中出现内存不足错误。

解决方案：

增加R的内存限制：

1
memory.limit(size = 8000)  # 将内存限制设置为8GB（仅适用于Windows）

关闭其他内存密集型应用程序
在具有更多RAM的计算机上尝试安装

4. 路径问题

问题：安装后无法加载nlmixr2或其依赖项。

解决方案：

检查R库路径：
1
.libPaths()
确保安装路径在R的库路径中
尝试在管理员/超级用户模式下安装

5. 版本兼容性问题

问题：nlmixr2与其他包或R版本不兼容。

解决方案：

更新R到最新版本
确保安装了兼容版本的依赖包
对于旧版本的R，尝试安装旧版本的nlmixr2

下一步

成功安装nlmixr2后，您可以开始学习基本概念并尝试一些简单的示例。在下一节中，我们将介绍非线性混合效应建模的基本概念，并展示如何使用nlmixr2进行简单的药代动力学模型拟合。

非线性混合效应建模基本概念

非线性混合效应建模（Nonlinear Mixed-Effects Modeling，NLME）是一种强大的统计方法，广泛应用于药代动力学（PK）、药效动力学（PD）和暴露-反应关系的分析。本节将详细介绍NLME的基本概念、数学基础以及在药物研发中的应用，为使用nlmixr2软件包进行建模奠定理论基础。

什么是非线性混合效应建模？

非线性混合效应建模是一种统计方法，用于分析具有非线性特性的纵向数据，特别是当数据来自多个受试者或实验单位时。在药物研发领域，这种方法通常被称为群体药代动力学/药效动力学（Population PK/PD）建模。

“非线性"指的是模型参数与观测结果之间的关系不是线性的，这在描述药物在体内的行为时非常常见。例如，药物的吸收、分布和消除通常遵循非线性过程。

“混合效应"则指模型同时包含固定效应（fixed effects）和随机效应（random effects）：

固定效应：代表整个群体的平均特性，如典型药物清除率
随机效应：描述个体间变异（between-subject variability，BSV）和个体内变异（within-subject variability，WSV）

非线性混合效应模型的数学表达

一个基本的非线性混合效应模型可以表示为：

$$y_{ij} = f(x_{ij}, \theta_i) + \varepsilon_{ij}$$

其中：

$y_{ij}$ 是第 $i$ 个受试者在第 $j$ 个时间点的观测值（如血药浓度）
$f$ 是描述系统行为的非线性函数
$x_{ij}$ 是自变量（如时间、剂量等）
$\theta_i$ 是第 $i$ 个受试者的参数向量
$\varepsilon_{ij}$ 是残差误差

个体参数 $\theta_i$ 可以进一步分解为：

$$\theta_i = g(\theta_{pop}, \eta_i, z_i)$$

其中：

$\theta_{pop}$ 是群体参数（固定效应）
$\eta_i$ 是描述个体间变异的随机效应，通常假设服从多元正态分布 $\eta_i \sim N(0, \Omega)$
$z_i$ 是协变量（如体重、年龄、性别等）
$g$ 是参数模型函数，通常采用指数形式以确保参数为正值

残差误差 $\varepsilon_{ij}$ 通常假设服从正态分布 $\varepsilon_{ij} \sim N(0, \sigma^2)$，可以采用不同的误差模型：

加性误差：$\varepsilon_{ij} \sim N(0, \sigma^2)$
比例误差：$\varepsilon_{ij} \sim N(0, (\sigma \cdot f(x_{ij}, \theta_i))^2)$
混合误差：结合加性和比例误差

常见的药代动力学模型

一室模型

最简单的PK模型是一室模型，将人体视为单一均质的区室，药物在其中均匀分布。

对于静脉注射给药，一室模型的微分方程为：

$$\frac{dA}{dt} = -k \cdot A$$

其中 $A$ 是药物量，$k$ 是消除速率常数。

解析解为：

$$C(t) = \frac{Dose}{V} \cdot e^{-k \cdot t}$$

其中 $C(t)$ 是时间 $t$ 的血药浓度，$V$ 是分布容积。

二室模型

二室模型将人体分为中央室（血液和灌注良好的组织）和外周室（灌注较差的组织）。

微分方程组为：

$$\frac{dA_1}{dt} = -k_{10} \cdot A_1 - k_{12} \cdot A_1 + k_{21} \cdot A_2$$

$$\frac{dA_2}{dt} = k_{12} \cdot A_1 - k_{21} \cdot A_2$$

其中：

$A_1$ 和 $A_2$ 分别是中央室和外周室的药物量
$k_{10}$ 是从中央室的消除速率常数
$k_{12}$ 和 $k_{21}$ 是室间转运速率常数

口服给药模型

对于口服给药，需要考虑吸收过程。一个简单的一阶吸收一室模型的微分方程为：

$$\frac{dA_{gut}}{dt} = -k_a \cdot A_{gut}$$

$$\frac{dA_{central}}{dt} = k_a \cdot A_{gut} - k \cdot A_{central}$$

其中 $k_a$ 是吸收速率常数。

药效动力学模型

药效动力学模型描述药物浓度与效应之间的关系。常见的PD模型包括：

Emax模型（Hill方程）

$$E = E_0 + \frac{E_{max} \cdot C^n}{EC_{50}^n + C^n}$$

其中：

$E$ 是药物效应
$E_0$ 是基线效应
$E_{max}$ 是最大效应
$C$ 是药物浓度
$EC_{50}$ 是产生50%最大效应的浓度
$n$ 是Hill系数，描述浓度-效应曲线的陡峭程度

间接效应模型

间接效应模型描述药物通过影响内源性物质的产生或消除来产生效应：

$$\frac{dR}{dt} = k_{in} \cdot (1 \pm \frac{E_{max} \cdot C}{EC_{50} + C}) - k_{out} \cdot R$$

其中：

$R$ 是效应物质的量
$k_{in}$ 是产生速率
$k_{out}$ 是消除速率

协变量模型

协变量模型用于解释个体间变异的部分原因，常见的协变量包括体重、年龄、性别、肾功能等。

一个典型的协变量模型例子是体重对清除率的影响：

$$CL_i = CL_{pop} \cdot (\frac{WT_i}{WT_{ref}})^{0.75} \cdot e^{\eta_{CL,i}}$$

其中：

$CL_i$ 是个体 $i$ 的清除率
$CL_{pop}$ 是群体典型清除率
$WT_i$ 是个体 $i$ 的体重
$WT_{ref}$ 是参考体重（通常为70kg）
$\eta_{CL,i}$ 是个体 $i$ 的随机效应

nlmixr2中的模型表示

nlmixr2提供了灵活的语法来定义非线性混合效应模型。一个基本的nlmixr2模型包括两个主要部分：

初始值部分（ini）：定义模型参数的初始估计值和随机效应
模型部分（model）：定义模型结构，包括微分方程、观测模型和误差模型

下面是一个简单的一室模型的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
27
28
29
30
31
32
33
one_cmt <- function() {
  ini({
    # 固定效应参数（群体参数）的初始值
    tka <- log(1.0)   # 吸收速率常数的对数
    tcl <- log(0.1)   # 清除率的对数
    tv <- log(10)     # 分布容积的对数
    
    # 随机效应（个体间变异）
    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
    
    # 计算血药浓度
    cp = central / v
    
    # 误差模型
    cp ~ add(add.err) + prop(prop.err)
  })
}

在这个例子中：

tka、tcl和tv是群体参数的对数值
eta.ka、eta.cl和eta.v是描述个体间变异的随机效应
微分方程描述了药物从给药部位（depot）到中央室（central）的转运和消除
观测模型将中央室的药物量转换为血药浓度
误差模型结合了加性和比例误差

非线性混合效应建模的优势

与传统的非房室分析（NCA）或单独分析每个受试者的数据相比，NLME具有以下优势：

数据整合：可以同时分析来自多个受试者、多个研究的数据
处理稀疏数据：即使每个受试者的采样点很少，也能得到可靠的估计
定量个体间变异：可以量化参数的变异程度
协变量分析：可以识别影响药物行为的因素
预测能力：可以预测不同剂量、不同人群的药物行为
模拟能力：可以模拟临床试验，优化试验设计

非线性混合效应建模的应用

NLME在药物研发的各个阶段都有广泛应用：

早期药物开发：
- 首次人体试验（FIH）剂量选择
- 确定药物的PK特性
- 评估剂量线性性
临床药理学：
- 特殊人群研究（如肾功能不全患者）
- 药物相互作用研究
- 暴露-反应关系分析
临床试验优化：
- 剂量方案优化
- 采样时间优化
- 临床试验模拟
监管申报：
- 支持剂量选择和给药方案
- 特殊人群剂量调整建议
- 标签说明支持

实际应用示例：华法林PK/PD模型

华法林是一种常用的抗凝血药物，其PK/PD关系复杂。下面是一个使用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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
warfarin_pkpd <- function() {
  ini({
    # PK参数
    tka <- log(1.0)     # 吸收速率常数的对数
    tcl <- log(0.1)     # 清除率的对数
    tv <- log(10)       # 分布容积的对数
    
    # PD参数
    tEC50 <- log(2)     # EC50的对数
    tEmax <- log(100)   # Emax的对数
    
    # 随机效应
    eta.ka ~ 0.4
    eta.cl ~ 0.3
    eta.v ~ 0.2
    eta.EC50 ~ 0.3
    eta.Emax ~ 0.2
    
    # 误差模型
    prop.err <- 0.1
    add.err.pd <- 0.2
  })
  model({
    # 个体PK参数
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    
    # 个体PD参数
    EC50 <- exp(tEC50 + eta.EC50)
    Emax <- exp(tEmax + eta.Emax)
    
    # PK模型
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    
    # 计算浓度
    cp = central / v
    
    # PD模型（Emax模型）
    effect = Emax * cp / (EC50 + cp)
    
    # 误差模型
    cp ~ prop(prop.err)
    effect ~ add(add.err.pd)
  })
}

这个模型包括：

一个一室PK模型，描述华法林的吸收和消除
一个Emax PD模型，描述华法林浓度与抗凝效应的关系
个体间变异和残差误差模型

总结

非线性混合效应建模是一种强大的统计方法，可以同时分析来自多个受试者的数据，量化个体间变异，并识别影响药物行为的因素。nlmixr2提供了一个灵活的R框架来实现这些模型，支持从简单的一室模型到复杂的PK/PD模型的各种应用。

在下一节中，我们将详细介绍nlmixr2中的参数估计方法，包括最大似然估计、条件估计和贝叶斯方法，以及如何选择和配置适当的估计算法。

nlmixr2中的参数估计方法

参数估计是非线性混合效应建模中的核心环节，它决定了模型的准确性和可靠性。nlmixr2提供了多种参数估计方法，以适应不同类型的模型和数据特点。本节将详细介绍nlmixr2中可用的各种参数估计算法、它们的理论基础、优缺点以及如何在实际应用中选择和配置适当的估计方法。

参数估计的基本原理

在非线性混合效应模型中，参数估计的目标是找到最能解释观测数据的模型参数值。这通常涉及最大化似然函数（或等价地，最小化目标函数）。

对于一个基本的非线性混合效应模型：

$$y_{ij} = f(x_{ij}, \theta_i) + \varepsilon_{ij}$$ $$\theta_i = g(\theta_{pop}, \eta_i, z_i)$$ $$\eta_i \sim N(0, \Omega)$$ $$\varepsilon_{ij} \sim N(0, \sigma^2)$$

参数估计的目标是找到群体参数 $\theta_{pop}$、协变量效应参数、随机效应协方差矩阵 $\Omega$ 和残差误差参数 $\sigma^2$ 的最优值。

nlmixr2中的估计方法概述

nlmixr2提供了多种估计方法，每种方法都有其特定的算法、优势和适用场景。主要的估计方法包括：

FOCE（First-Order Conditional Estimation）：条件一阶估计方法
SAEM（Stochastic Approximation Expectation Maximization）：随机近似期望最大化算法
FOCEI（First-Order Conditional Estimation with Interaction）：带交互项的条件一阶估计
nlme：基于R的nlme包的估计方法
posthoc：后验估计方法
FOCEi：带交互项的条件一阶估计（nlmixr2的实现）
FOCE：条件一阶估计（nlmixr2的实现）
FO：一阶估计方法

下面我们将详细介绍这些方法的原理、特点和使用方法。

FOCE和FOCEI方法

原理

FOCE（First-Order Conditional Estimation）和FOCEI（First-Order Conditional Estimation with Interaction）是最常用的参数估计方法之一，也是NONMEM中的标准方法。

这些方法基于对似然函数的近似：

对随机效应 $\eta$ 进行泰勒展开，保留一阶项
对每个个体，条件于当前的 $\eta$ 估计值计算似然函数
迭代优化群体参数和个体参数

FOCEI比FOCE多考虑了随机效应与残差误差之间的交互作用，通常在非线性程度较高的模型中表现更好。

nlmixr2中的使用

在nlmixr2中，可以通过以下方式使用FOCE或FOCEI方法：

 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
# 定义模型
mod <- function() {
  ini({
    tka <- log(1.0)
    tcl <- log(0.1)
    tv <- log(10)
    eta.ka ~ 0.4
    eta.cl ~ 0.3
    eta.v ~ 0.2
    add.err <- 0.1
  })
  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
    cp = central / v
    cp ~ add(add.err)
  })
}

# 使用FOCEI方法拟合模型
fit_focei <- nlmixr(mod, data, est="focei")

# 使用FOCE方法拟合模型
fit_foce <- nlmixr(mod, data, est="foce")

控制参数

FOCEI和FOCE方法可以通过foceiControl()和foceControl()函数进行配置：

1
2
3
4
5
6
7
8
9
fit_focei <- nlmixr(mod, data, est="focei", 
                    control=foceiControl(
                      maxOuterIterations=100,  # 最大外层迭代次数
                      maxInnerIterations=1000, # 最大内层迭代次数
                      atol=1e-8,              # 绝对容差
                      rtol=1e-8,              # 相对容差
                      covMethod="r,s",        # 协方差矩阵计算方法
                      interaction=TRUE        # 是否包含交互项
                    ))

主要控制参数包括：

maxOuterIterations：最大外层迭代次数
maxInnerIterations：最大内层迭代次数
atol和rtol：ODE求解器的绝对和相对容差
covMethod：协方差矩阵计算方法（“r,s"表示使用R矩阵和S矩阵）
interaction：是否包含交互项（FOCEI为TRUE，FOCE为FALSE）

优缺点

优点：

计算效率较高，特别是对于简单模型
与NONMEM兼容，结果可比较
对于中等复杂度的模型表现良好

缺点：

对于高度非线性模型可能不稳定
对初始值敏感
可能陷入局部最优解

SAEM方法

原理

SAEM（Stochastic Approximation Expectation Maximization）是一种基于EM算法的随机近似方法，特别适用于复杂的非线性混合效应模型。

SAEM算法的基本步骤：

初始化：设定初始参数值
模拟步骤：对每个个体，从条件分布中抽样随机效应
随机近似：更新参数估计，使用当前和历史估计的加权平均
重复：迭代直到收敛

SAEM不需要对似然函数进行线性化近似，因此对于高度非线性模型更稳健。

nlmixr2中的使用

在nlmixr2中，可以通过以下方式使用SAEM方法：

1
2
# 使用SAEM方法拟合模型
fit_saem <- nlmixr(mod, data, est="saem")

控制参数

SAEM方法可以通过saemControl()函数进行配置：

1
2
3
4
5
6
7
8
9
fit_saem <- nlmixr(mod, data, est="saem", 
                   control=saemControl(
                     nBurn=200,         # 燃烧期迭代次数
                     nEm=300,           # EM迭代次数
                     nmc=3,             # 每次迭代的蒙特卡洛样本数
                     seed=12345,        # 随机数种子
                     print=50,          # 每50次迭代打印一次
                     covMethod="linFim" # 协方差矩阵计算方法
                   ))

主要控制参数包括：

nBurn：燃烧期迭代次数，在此期间算法探索参数空间
nEm：EM迭代次数，在此期间算法收敛到最终估计
nmc：每次迭代的蒙特卡洛样本数
seed：随机数种子，用于结果复现
print：打印频率
covMethod：协方差矩阵计算方法

优缺点

优点：

对高度非线性模型更稳健
不需要对似然函数进行线性化近似
不易陷入局部最优解
计算效率高，特别是对于复杂模型

缺点：

是一种随机算法，结果可能有轻微变化
需要调整多个控制参数以获得最佳性能
对于简单模型可能计算量过大

nlme方法

原理

nlme方法基于R的nlme包，使用伪似然方法进行参数估计。它是一种成熟的非线性混合效应模型估计方法，但在nlmixr2中主要用于简单模型或初步分析。

nlmixr2中的使用

1
2
# 使用nlme方法拟合模型
fit_nlme <- nlmixr(mod, data, est="nlme")

控制参数

nlme方法可以通过nlmeControl()函数进行配置：

1
2
3
4
5
6
fit_nlme <- nlmixr(mod, data, est="nlme", 
                   control=nlmeControl(
                     msVerbose=TRUE,     # 详细输出
                     returnObject=TRUE,  # 返回nlme对象
                     pnlsTol=0.001       # pnls容差
                   ))

优缺点

优点：

基于成熟的R包nlme
对于简单模型计算快速
与R的nlme包兼容

缺点：

功能有限，不支持某些高级特性
对于复杂模型可能不稳定
在nlmixr2中主要用于初步分析

其他估计方法

posthoc方法

posthoc方法用于在已知群体参数的情况下估计个体参数（随机效应）。它通常用于模型诊断或基于先前估计的群体参数进行个体预测。

1
2
3
4
5
6
7
8
9
# 使用posthoc方法估计个体参数
fit_posthoc <- nlmixr(mod, data, est="posthoc", 
                      control=posthocControl(
                        thetaFixed=c(1.0, 0.1, 10.0), # 固定群体参数
                        omega=matrix(c(0.16, 0, 0, 
                                      0, 0.09, 0, 
                                      0, 0, 0.04), 3, 3), # 固定随机效应协方差矩阵
                        sigma=0.01                       # 固定残差误差
                      ))

FO方法

FO（First-Order）是最简单的估计方法，它在随机效应的均值（通常为零）处线性化模型。这种方法计算速度快，但准确性较低，主要用于初步分析或非常简单的模型。

1
2
# 使用FO方法拟合模型
fit_fo <- nlmixr(mod, data, est="fo")

如何选择合适的估计方法

选择合适的估计方法取决于多种因素，包括：

模型复杂度：
- 简单模型：FOCE、FOCEI或FO可能足够
- 复杂模型：SAEM通常更稳健
数据特性：
- 稀疏数据：SAEM可能表现更好
- 丰富数据：FOCEI通常足够
计算资源：
- 有限资源：FOCE或FOCEI可能更快
- 充足资源：SAEM可能提供更准确的结果
与其他软件比较的需求：
- 需要与NONMEM比较：使用FOCEI
- 需要与Monolix比较：使用SAEM
模型开发阶段：
- 初步探索：FO或nlme可能足够快
- 最终模型：FOCEI或SAEM提供更准确的结果

实际应用示例

下面是一个使用不同估计方法拟合同一模型的比较示例：

 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
library(nlmixr2)
library(ggplot2)

# 加载示例数据
data(theo_sd)
theo_data <- theo_sd

# 定义一个一室模型
one_cmt_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  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
    cp = central / v
    cp ~ add(add.err)
  })
}

# 使用不同方法拟合模型
fit_focei <- nlmixr(one_cmt_model, theo_data, est="focei")
fit_foce <- nlmixr(one_cmt_model, theo_data, est="foce")
fit_saem <- nlmixr(one_cmt_model, theo_data, est="saem")
fit_fo <- nlmixr(one_cmt_model, theo_data, est="fo")

# 比较参数估计结果
compare_fits <- function(fits, names) {
  params <- c("tka", "tcl", "tv", "add.err")
  results <- data.frame(Parameter = params)
  
  for (i in 1:length(fits)) {
    fit <- fits[[i]]
    name <- names[i]
    estimates <- fixef(fit)[params]
    results[[name]] <- estimates
  }
  
  return(results)
}

fits <- list(fit_focei, fit_foce, fit_saem, fit_fo)
names <- c("FOCEI", "FOCE", "SAEM", "FO")
comparison <- compare_fits(fits, names)
print(comparison)

# 比较目标函数值
obj_values <- data.frame(
  Method = names,
  OBJ = sapply(fits, function(fit) fit$objf)
)
print(obj_values)

# 可视化比较
plot_comparison <- function(fit, name) {
  p <- plot(fit) + ggtitle(paste("Method:", name))
  return(p)
}

plots <- lapply(1:length(fits), function(i) plot_comparison(fits[[i]], names[i]))
gridExtra::grid.arrange(grobs=plots, ncol=2)

参数估计的诊断和评估

在选择估计方法后，评估估计质量非常重要。nlmixr2提供了多种诊断工具：

参数估计的标准误差：

1
summary(fit)  # 显示参数估计及其标准误差

协方差矩阵和相关矩阵：

1
2
vcov(fit)     # 显示参数估计的协方差矩阵
cov2cor(vcov(fit))  # 显示参数估计的相关矩阵

条件数：

1
kappa(fit)    # 显示协方差矩阵的条件数，评估参数估计的稳定性

收敛诊断：

1
fit$message   # 显示优化器的收敛消息

目标函数值分解：

1
2
fit$objf      # 总目标函数值
fit$objfComponents  # 目标函数的各个组成部分

高级参数估计技巧

1. 使用多种初始值

对于复杂模型，使用多种初始值可以避免陷入局部最优解：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# 定义不同的初始值
ini1 <- c(tka=log(1.0), tcl=log(0.1), tv=log(10))
ini2 <- c(tka=log(2.0), tcl=log(0.2), tv=log(15))

# 使用不同初始值拟合
fit1 <- nlmixr(mod, data, est="focei", 
               control=foceiControl(theta=ini1))
fit2 <- nlmixr(mod, data, est="focei", 
               control=foceiControl(theta=ini2))

# 比较目标函数值
if (fit1$objf < fit2$objf) {
  best_fit <- fit1
} else {
  best_fit <- fit2
}

2. 分步估计策略

对于复杂模型，分步估计可能更稳定：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# 第一步：固定随机效应，仅估计群体参数
fit1 <- nlmixr(mod, data, est="focei", 
               control=foceiControl(
                 skipCov=TRUE,  # 跳过协方差计算
                 omegaIsChol=FALSE,  # omega不是Cholesky分解
                 omega=matrix(0, 3, 3)  # 固定omega为零矩阵
               ))

# 第二步：使用第一步的估计作为初始值，估计所有参数
fit2 <- nlmixr(mod, data, est="focei", 
               control=foceiControl(
                 theta=fixef(fit1)  # 使用第一步的固定效应估计
               ))

3. 使用不同的优化算法

nlmixr2支持多种优化算法，可以根据模型特点选择：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# 使用不同的优化器
fit_nlminb <- nlmixr(mod, data, est="focei", 
                     control=foceiControl(
                       optimMethod="nlminb"  # 使用nlminb优化器
                     ))

fit_lbfgsb <- nlmixr(mod, data, est="focei", 
                     control=foceiControl(
                       optimMethod="L-BFGS-B"  # 使用L-BFGS-B优化器
                     ))

总结

nlmixr2提供了多种参数估计方法，每种方法都有其特定的优势和适用场景。选择合适的估计方法取决于模型复杂度、数据特性、计算资源和研究目的。

FOCEI：适用于中等复杂度的模型，与NONMEM兼容
SAEM：适用于复杂的非线性模型，更稳健但计算量大
其他方法（FO、FOCE、nlme、posthoc）：在特定场景下有用

在实际应用中，建议尝试多种估计方法并比较结果，结合模型诊断工具评估估计质量，并根据研究目的选择最合适的方法。

在下一节中，我们将详细介绍nlmixr2中的模型诊断技术，包括残差分析、VPC（视觉预测检验）和引导法等，以评估模型的适合度和预测能力。

nlmixr2中的模型诊断技术

模型诊断是药代动力学/药效动力学建模过程中的关键步骤，用于评估模型的适合度、预测能力和假设的合理性。良好的模型诊断可以帮助识别模型的缺陷，指导模型改进，并增强对模型预测的信心。本节将详细介绍nlmixr2中可用的各种模型诊断工具、图形和统计方法，以及如何解释诊断结果。

模型诊断的目的

模型诊断的主要目的包括：

评估模型的适合度：模型是否能够准确描述观测数据
检验模型假设：如随机效应的正态性、残差的独立性等
识别异常值和有影响的观测值：可能影响参数估计的数据点
评估模型的预测能力：模型对未见数据的预测准确性
比较不同模型：为模型选择提供客观依据

nlmixr2中的基本诊断工具

1. 模型拟合摘要

拟合模型后，首先应查看模型的基本摘要信息：

 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
library(nlmixr2)
data(theo_sd)

# 定义模型
one_cmt_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  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
    cp = central / v
    cp ~ add(add.err)
  })
}

# 拟合模型
fit <- nlmixr(one_cmt_model, theo_sd, est="focei")

# 查看摘要
summary(fit)

摘要信息包括：

固定效应参数估计及其标准误差
随机效应参数（方差-协方差矩阵）
残差误差参数
目标函数值和信息准则（AIC、BIC）
条件数和其他诊断统计量

2. 参数估计的不确定性

评估参数估计的不确定性是模型诊断的重要部分：

1
2
3
4
5
6
7
8
# 查看参数估计的协方差矩阵
vcov(fit)

# 查看参数估计的相关矩阵
cov2cor(vcov(fit))

# 查看参数估计的置信区间
confint(fit)

关键诊断指标包括：

标准误差：较大的标准误差表明参数估计不精确
相关系数：参数间高度相关（>0.9或<-0.9）可能表明过度参数化
条件数：高条件数（>1000）表明参数估计不稳定

3. 基本拟合图

nlmixr2提供了多种基本拟合图来评估模型性能：

1
2
3
4
5
6
7
8
# 基本诊断图
plot(fit)

# 个体拟合图
plot(fit, type="individual")

# 观测值vs预测值图
plot(fit, type="pred")

基本拟合图包括：

观测值vs群体预测值（PRED）：评估群体模型的整体拟合
观测值vs个体预测值（IPRED）：评估个体模型的拟合
条件加权残差（CWRES）vs时间：评估模型随时间的拟合
条件加权残差vs预测值：评估模型在不同浓度水平的拟合
个体拟合图：评估每个受试者的模型拟合

高级诊断技术

1. 条件加权残差（CWRES）分析

CWRES是评估模型适合度的重要工具，理想情况下应该服从标准正态分布：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# 添加CWRES
fit <- addCwres(fit)

# CWRES vs 时间
plot(fit, type="cwres")

# CWRES的QQ图
plot(fit, type="qq")

# CWRES的直方图
hist(fit$CWRES, breaks=20, main="CWRES Histogram", 
     xlab="CWRES", freq=FALSE)
curve(dnorm(x, 0, 1), add=TRUE, col="red", lwd=2)

# CWRES的统计检验
shapiro.test(fit$CWRES)  # 正态性检验

CWRES诊断的关键点：

CWRES应该随机分布在0附近，无明显趋势
CWRES的范围通常在-2到2之间，超过±3的值可能是异常值
CWRES的QQ图应接近对角线
CWRES的直方图应接近标准正态分布

2. 归一化预测分布误差（NPDE）

NPDE是一种基于模拟的诊断方法，特别适用于稀疏数据和非线性模型：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# 添加NPDE
fit <- addNpde(fit)

# NPDE vs 时间
plot(fit, type="npde")

# NPDE的QQ图和直方图
plot(fit, type="npde-qq")
plot(fit, type="npde-hist")

# NPDE的统计检验
shapiro.test(fit$NPDE)  # 正态性检验
t.test(fit$NPDE)        # 均值是否为0
var.test(fit$NPDE, rnorm(length(fit$NPDE)))  # 方差是否为1

NPDE诊断的关键点：

NPDE应该服从标准正态分布（均值0，方差1）
NPDE vs 时间或预测值不应显示任何趋势
NPDE的QQ图应接近对角线
NPDE的直方图应接近标准正态分布

3. 视觉预测检验（VPC）

VPC是一种强大的模型评估工具，通过比较观测数据和模型模拟数据的分布来评估模型的预测能力：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# 基本VPC
vpc_fit <- vpcPlot(fit, n=500, show=list(obs_dv=TRUE))
print(vpc_fit)

# 分层VPC（按协变量分组）
vpc_strat <- vpcPlot(fit, n=500, stratify=c("SEX"), 
                    show=list(obs_dv=TRUE))
print(vpc_strat)

# 预测校正VPC
vpc_pc <- vpcPlot(fit, n=500, pred_corr=TRUE, 
                 show=list(obs_dv=TRUE))
print(vpc_pc)

VPC的关键点：

观测数据的分位数（通常是5%、50%、95%）应该落在相应模拟分位数的置信区间内
观测数据点应该均匀分布在模拟区域内
任何系统性偏差都表明模型可能存在问题
预测校正VPC对于非线性模型或有协变量影响的模型特别有用

4. 引导法（Bootstrap）

引导法通过重复抽样来评估参数估计的稳健性和不确定性：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# 执行引导分析（这可能需要较长时间）
boot_fit <- bootplot(fit, nboot=100, seed=12345)

# 查看引导结果
print(boot_fit)

# 绘制引导参数分布
plot(boot_fit)

# 计算引导置信区间
boot_ci <- boot_fit$confint
print(boot_ci)

引导法诊断的关键点：

参数的引导分布应该接近正态分布
原始估计值应该接近引导分布的中位数
引导失败率应该较低（<10%）
引导置信区间提供了参数不确定性的非参数估计

5. 协变量筛选和模型比较

在模型开发过程中，需要比较不同模型并评估协变量的影响：

 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
# 定义基础模型和协变量模型
base_model <- function() {
  # 基础模型定义...
}

cov_model <- function() {
  # 包含协变量的模型定义...
}

# 拟合两个模型
fit_base <- nlmixr(base_model, data, est="focei")
fit_cov <- nlmixr(cov_model, data, est="focei")

# 比较模型
compare_models <- function(fits, names) {
  results <- data.frame(
    Model = names,
    OBJ = sapply(fits, function(fit) fit$objf),
    AIC = sapply(fits, function(fit) fit$AIC),
    BIC = sapply(fits, function(fit) fit$BIC)
  )
  return(results)
}

fits <- list(fit_base, fit_cov)
names <- c("Base Model", "Covariate Model")
comparison <- compare_models(fits, names)
print(comparison)

# 似然比检验（对于嵌套模型）
dof <- length(fixef(fit_cov)) - length(fixef(fit_base))
p_value <- 1 - pchisq(fit_base$objf - fit_cov$objf, df=dof)
cat("Likelihood Ratio Test p-value:", p_value, "\n")

模型比较的关键点：

目标函数值（OBJ）的显著降低表明模型改进
对于嵌套模型，可以使用似然比检验（差异在统计上是否显著）
对于非嵌套模型，可以使用AIC或BIC（较小值表示更好的模型）
模型复杂性增加应该带来显著的拟合改善

使用xpose.nlmixr2进行高级诊断

xpose.nlmixr2包提供了与xpose兼容的诊断工具，可以生成更多高级诊断图：

 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(xpose)
library(xpose.nlmixr2)

# 将nlmixr2拟合对象转换为xpose对象
xp_fit <- xpose_data_nlmixr(fit)

# 基本诊断图
dv_vs_pred(xp_fit)
dv_vs_ipred(xp_fit)
res_vs_pred(xp_fit)
res_vs_idv(xp_fit)

# 残差分布图
res_distrib(xp_fit)
res_qq(xp_fit)

# 参数分布图
prm_distrib(xp_fit)
prm_qq(xp_fit)

# 协变量关系图
cov_splom(xp_fit)

# 个体随机效应图
ranef_qq(xp_fit)
ranef_vs_cov(xp_fit, covariates=c("WT", "AGE"))

xpose.nlmixr2的优势：

提供更多样化和可定制的诊断图
与xpose包兼容，可以使用xpose的所有功能
支持交互式图形和分面图

诊断结果的解释和模型改进

基于诊断结果，可以识别模型的问题并进行相应改进：

1. 结构模型问题

如果诊断显示结构模型问题（如CWRES vs 时间有明显趋势），可能的解决方案包括：

增加或减少房室数量
修改吸收模型（如从一阶吸收改为零阶吸收或迟滞吸收）
添加非线性清除或分布过程
考虑时变参数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# 修改为二室模型示例
two_cmt_model <- function() {
  ini({
    # 参数初始值...
  })
  model({
    # 参数转换...
    
    # 二室模型微分方程
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - (cl/v) * central - (q/v) * central + (q/v2) * peripheral
    d/dt(peripheral) = (q/v) * central - (q/v2) * peripheral
    
    cp = central / v
    cp ~ add(add.err)
  })
}

2. 统计模型问题

如果诊断显示统计模型问题（如异方差性残差），可能的解决方案包括：

修改残差误差模型（如从加性改为比例或组合）
修改随机效应结构（如添加或移除协方差）
考虑不同的分布假设

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# 修改误差模型示例
error_model <- function() {
  ini({
    # 参数初始值...
    prop.err <- 0.1
    add.err <- 0.5
  })
  model({
    # 模型定义...
    
    # 组合误差模型
    cp ~ add(add.err) + prop(prop.err)
  })
}

3. 协变量模型问题

如果诊断显示协变量关系问题（如随机效应与协变量相关），可能的解决方案包括：

添加新的协变量关系
修改现有协变量关系的函数形式
考虑协变量之间的交互作用

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 添加体重对清除率影响的示例
cov_model <- function() {
  ini({
    # 参数初始值...
    beta_cl_wt <- 0.75  # 体重对清除率的影响系数
  })
  model({
    # 基本参数...
    
    # 添加协变量影响
    cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt
    
    # 其余模型...
  })
}

4. 数据问题

如果诊断显示数据问题（如异常值），可能的解决方案包括：

检查并纠正数据录入错误
考虑排除或降权异常值
添加缺失数据处理策略

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# 使用IF条件处理异常值的示例
outlier_model <- function() {
  ini({
    # 参数初始值...
  })
  model({
    # 模型定义...
    
    # 使用不同误差模型处理可能的异常值
    if (DVID == 1) {
      cp ~ add(add.err1)
    } else {
      cp ~ add(add.err2)
    }
  })
}

模型诊断的工作流程

一个完整的模型诊断工作流程通常包括以下步骤：

基本诊断：
- 检查参数估计的合理性和精确度
- 查看基本拟合图（观测值vs预测值，残差图）
- 评估目标函数值和信息准则
残差分析：
- 检查CWRES和NPDE的分布和趋势
- 识别潜在的模型误设和异常值
预测检验：
- 执行VPC分析评估模型预测能力
- 考虑分层VPC以评估协变量影响
稳健性分析：
- 使用引导法评估参数估计的稳健性
- 检查参数相关性和条件数
模型比较：
- 比较不同模型结构和复杂度
- 使用统计检验和信息准则选择最佳模型
模型改进：
- 基于诊断结果修改模型
- 重复诊断过程直到获得满意的模型

实际应用示例：茶碱PK模型诊断

下面是一个使用nlmixr2对茶碱PK数据进行模型诊断的完整示例：

 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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
library(nlmixr2)
library(ggplot2)
library(xpose.nlmixr2)
library(gridExtra)

# 加载茶碱数据
data(theo_sd)
theo_data <- theo_sd

# 定义一室模型
one_cmt_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  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
    cp = central / v
    cp ~ add(add.err)
  })
}

# 拟合模型
fit <- nlmixr(one_cmt_model, theo_data, est="focei")

# 1. 基本诊断
cat("Model summary:\n")
print(summary(fit))

cat("\nParameter estimates and confidence intervals:\n")
print(confint(fit))

# 2. 残差分析
fit <- addCwres(fit)
fit <- addNpde(fit)

# 创建诊断图
p1 <- plot(fit, type="pred")
p2 <- plot(fit, type="ipred")
p3 <- plot(fit, type="cwres")
p4 <- plot(fit, type="npde")

grid.arrange(p1, p2, p3, p4, ncol=2)

# CWRES和NPDE的统计检验
cat("\nCWRES normality test:\n")
print(shapiro.test(fit$CWRES))

cat("\nNPDE normality test:\n")
print(shapiro.test(fit$NPDE))

cat("\nNPDE mean test (should be close to 0):\n")
print(t.test(fit$NPDE))

cat("\nNPDE variance test (should be close to 1):\n")
print(var.test(fit$NPDE, rnorm(length(fit$NPDE))))

# 3. 视觉预测检验
vpc_result <- vpcPlot(fit, n=500, show=list(obs_dv=TRUE))
print(vpc_result)

# 4. 个体拟合图
ind_plots <- plot(fit, type="individual", ncol=3, nrow=3)
print(ind_plots)

# 5. 使用xpose进行高级诊断
xp_fit <- xpose_data_nlmixr(fit)

p5 <- dv_vs_pred(xp_fit)
p6 <- dv_vs_ipred(xp_fit)
p7 <- res_vs_pred(xp_fit)
p8 <- res_vs_idv(xp_fit)

grid.arrange(p5, p6, p7, p8, ncol=2)

p9 <- res_distrib(xp_fit)
p10 <- res_qq(xp_fit)
p11 <- prm_distrib(xp_fit)
p12 <- ranef_qq(xp_fit)

grid.arrange(p9, p10, p11, p12, ncol=2)

# 6. 诊断结论
cat("\nDiagnostic conclusions:\n")
cat("1. Parameter estimates are reasonable with acceptable precision.\n")
cat("2. Residual plots show no major trends, suggesting adequate structural model.\n")
cat("3. CWRES and NPDE distributions are approximately normal.\n")
cat("4. VPC shows good agreement between observations and predictions.\n")
cat("5. Individual fits capture the observed data well.\n")
cat("6. Overall, the one-compartment model appears adequate for this dataset.\n")

常见问题和解决方案

1. 残差显示明显趋势

问题：CWRES vs 时间或预测值显示明显趋势。

解决方案：

检查结构模型是否合适（如房室数量）
考虑非线性过程（如饱和清除）
检查时变参数的可能性

2. 个体拟合不佳

问题：某些个体的拟合明显不佳。

解决方案：

检查这些个体的特殊特征（如极端协变量值）
考虑添加相关协变量
检查是否存在数据问题或异常值

3. VPC显示系统性偏差

问题：VPC中观测分位数落在模拟置信区间外。

解决方案：

检查剂量或采样时间的系统性错误
考虑更复杂的误差模型
尝试预测校正VPC
检查是否需要分层VPC（按协变量分组）

4. 参数估计不精确

问题：参数估计的标准误差很大或置信区间很宽。

解决方案：

检查模型是否过度参数化
考虑固定某些难以估计的参数
使用更多信息性先验（如基于文献的先验）
收集更多数据或优化采样设计

5. 随机效应与协变量相关

问题：随机效应与协变量显示相关性。

解决方案：

将相关协变量纳入模型
尝试不同的协变量关系函数形式
检查是否有未考虑的交互作用

总结

模型诊断是药代动力学/药效动力学建模过程中不可或缺的一部分。nlmixr2提供了丰富的诊断工具，从基本的残差分析到高级的视觉预测检验和引导法，帮助研究者全面评估模型性能并指导模型改进。

良好的模型诊断实践包括：

使用多种诊断工具评估模型的不同方面
结合图形和统计方法进行全面评估
基于诊断结果系统地改进模型
记录诊断过程和决策理由

在下一节中，我们将详细介绍nlmixr2中的模型仿真方法，包括如何使用已建立的模型进行临床试验设计、剂量优化和特殊人群预测。

nlmixr2中的模型仿真方法

模型仿真是药代动力学/药效动力学（PK/PD）建模中的关键应用，它允许研究者利用已建立的模型预测未来或假设场景下的药物行为。在药物开发过程中，仿真可以用于优化临床试验设计、探索不同剂量方案、预测特殊人群的药物暴露，以及支持监管决策。本节将详细介绍如何使用nlmixr2进行模型仿真，包括基本仿真技术、高级应用和实际案例。

模型仿真的基本原理

模型仿真基于以下基本步骤：

建立模型：使用观测数据拟合PK/PD模型并估计参数
设计仿真场景：定义要探索的剂量、给药方案、人群特征等
生成仿真数据：基于模型和设计的场景生成预测数据
分析仿真结果：评估结果并得出结论

仿真可以在两个层面进行：

确定性仿真：仅使用群体典型参数值进行预测，不考虑变异性
随机性仿真：考虑参数变异性（个体间变异、个体内变异和残差误差），生成多个虚拟受试者或试验

nlmixr2中的基本仿真功能

nlmixr2通过与rxode2包的集成提供了强大的仿真功能。以下是基本仿真的步骤：

1. 从拟合模型中提取参数

首先，我们需要从已拟合的模型中提取参数：

 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
library(nlmixr2)
library(rxode2)
library(ggplot2)

# 加载示例数据
data(theo_sd)

# 定义并拟合模型
one_cmt_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  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
    cp = central / v
    cp ~ add(add.err)
  })
}

fit <- nlmixr(one_cmt_model, theo_sd, est="focei")

# 提取参数
params <- fixef(fit)  # 固定效应参数
omega <- fit$omega    # 随机效应协方差矩阵
sigma <- fit$sigma    # 残差误差参数

2. 创建仿真模型

接下来，我们需要创建用于仿真的rxode2模型：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# 创建rxode2模型
sim_model <- rxode2({
  # 参数
  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
  
  # 输出变量
  cp = central / v
})

3. 设置仿真场景

然后，我们需要定义仿真场景，包括剂量、给药时间和观测时间：

1
2
3
4
5
6
7
# 创建剂量事件表
dose_regimen <- et(timeUnits="hours", amtUnits="mg") %>%
  et(amt=100, cmt=1, ii=24, addl=9) %>%  # 100mg，每24小时一次，共10次
  et(0:240)  # 观测时间点，从0到240小时

# 查看剂量事件表
print(dose_regimen)

4. 执行确定性仿真

首先，我们可以执行确定性仿真，只使用群体典型参数值：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# 设置参数值（使用拟合模型的固定效应参数）
theta <- c(
  tka = params["tka"],
  tcl = params["tcl"],
  tv = params["tv"]
)

# 执行确定性仿真
sim_det <- sim_model$solve(theta, dose_regimen)

# 绘制结果
plot(sim_det, cp ~ time, main="确定性仿真", xlab="时间 (小时)", ylab="血药浓度")

5. 执行随机性仿真

接下来，我们可以执行随机性仿真，考虑参数变异性：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# 设置仿真参数
nsub <- 100  # 虚拟受试者数量
seed <- 12345  # 随机数种子

# 创建参数样本
set.seed(seed)
sim_params <- rxode2::rxInits(fit, nsub)

# 执行随机性仿真
sim_stoch <- sim_model$solve(sim_params, dose_regimen)

# 绘制结果
ggplot(sim_stoch, aes(x=time, y=cp, group=id)) +
  geom_line(alpha=0.2) +
  geom_line(data=sim_det, aes(x=time, y=cp), color="red", linewidth=1) +
  labs(title="随机性仿真", x="时间 (小时)", y="血药浓度") +
  theme_bw()

高级仿真技术

1. 添加协变量效应

在实际应用中，我们通常需要考虑协变量对药物行为的影响：

 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
# 定义包含协变量的模型
cov_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    beta_cl_wt <- 0.75  # 体重对清除率的影响系数
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt  # 体重对清除率的影响
    v <- exp(tv + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    cp = central / v
    cp ~ add(add.err)
  })
}

# 假设我们有体重数据
theo_sd$WT <- rnorm(nrow(theo_sd), 70, 15)  # 模拟体重数据

# 拟合模型
fit_cov <- nlmixr(cov_model, theo_sd, est="focei")

# 创建用于仿真的rxode2模型
sim_cov_model <- rxode2({
  # 参数
  ka <- exp(tka + eta.ka)
  cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt
  v <- exp(tv + eta.v)
  
  # 微分方程
  d/dt(depot) = -ka * depot
  d/dt(central) = ka * depot - cl/v * central
  
  # 输出变量
  cp = central / v
})

# 创建不同体重的虚拟人群
nsub <- 300
set.seed(seed)
sim_params_cov <- rxode2::rxInits(fit_cov, nsub)

# 添加不同体重组
wt_groups <- c(50, 70, 90)
sim_params_cov$WT <- rep(wt_groups, each=nsub/length(wt_groups))
sim_params_cov$WT_GROUP <- factor(sim_params_cov$WT)

# 执行仿真
sim_cov <- sim_cov_model$solve(sim_params_cov, dose_regimen)

# 绘制不同体重组的结果
ggplot(sim_cov, aes(x=time, y=cp, group=interaction(id, WT_GROUP), color=WT_GROUP)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=WT_GROUP), fun=mean, geom="line", linewidth=1) +
  labs(title="不同体重组的药物浓度", x="时间 (小时)", y="血药浓度", color="体重 (kg)") +
  theme_bw()

2. 剂量-暴露-反应关系仿真

PK/PD模型的一个重要应用是评估剂量-暴露-反应关系：

  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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# 定义PK/PD模型
pkpd_model <- function() {
  ini({
    # PK参数
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    
    # PD参数
    tEC50 <- log(2)
    tEmax <- log(100)
    
    # 随机效应
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    eta.EC50 ~ 0.3
    eta.Emax ~ 0.2
    
    # 误差模型
    prop.err <- 0.1
    add.err.pd <- 0.2
  })
  model({
    # PK参数
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    
    # PD参数
    EC50 <- exp(tEC50 + eta.EC50)
    Emax <- exp(tEmax + eta.Emax)
    
    # PK模型
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    
    # 计算浓度
    cp = central / v
    
    # PD模型（Emax模型）
    effect = Emax * cp / (EC50 + cp)
    
    # 误差模型
    cp ~ prop(prop.err)
    effect ~ add(add.err.pd)
  })
}

# 假设我们已经拟合了模型
# fit_pkpd <- nlmixr(pkpd_model, pkpd_data, est="focei")

# 为了演示，我们手动设置参数
theta_pkpd <- c(
  tka = log(1.5),
  tcl = log(0.08),
  tv = log(0.5),
  tEC50 = log(2),
  tEmax = log(100)
)

omega_pkpd <- matrix(0, 5, 5)
diag(omega_pkpd) <- c(0.36, 0.09, 0.01, 0.09, 0.04)

# 创建用于仿真的rxode2模型
sim_pkpd_model <- rxode2({
  # PK参数
  ka <- exp(tka + eta.ka)
  cl <- exp(tcl + eta.cl)
  v <- exp(tv + eta.v)
  
  # PD参数
  EC50 <- exp(tEC50 + eta.EC50)
  Emax <- exp(tEmax + eta.Emax)
  
  # PK模型
  d/dt(depot) = -ka * depot
  d/dt(central) = ka * depot - cl/v * central
  
  # 计算浓度和效应
  cp = central / v
  effect = Emax * cp / (EC50 + cp)
})

# 创建不同剂量方案
dose_levels <- c(50, 100, 200)
dose_regimens <- lapply(dose_levels, function(dose) {
  et(timeUnits="hours", amtUnits="mg") %>%
    et(amt=dose, cmt=1, ii=24, addl=9) %>%
    et(0:240)
})

# 生成参数样本
nsub <- 100
set.seed(seed)
sim_params_pkpd <- list(
  tka = theta_pkpd["tka"] + rnorm(nsub, 0, sqrt(omega_pkpd[1,1])),
  tcl = theta_pkpd["tcl"] + rnorm(nsub, 0, sqrt(omega_pkpd[2,2])),
  tv = theta_pkpd["tv"] + rnorm(nsub, 0, sqrt(omega_pkpd[3,3])),
  tEC50 = theta_pkpd["tEC50"] + rnorm(nsub, 0, sqrt(omega_pkpd[4,4])),
  tEmax = theta_pkpd["tEmax"] + rnorm(nsub, 0, sqrt(omega_pkpd[5,5])),
  eta.ka = rep(0, nsub),
  eta.cl = rep(0, nsub),
  eta.v = rep(0, nsub),
  eta.EC50 = rep(0, nsub),
  eta.Emax = rep(0, nsub)
)

# 执行不同剂量的仿真
sim_results <- list()
for (i in 1:length(dose_levels)) {
  sim_results[[i]] <- sim_pkpd_model$solve(sim_params_pkpd, dose_regimens[[i]])
  sim_results[[i]]$dose <- dose_levels[i]
}

# 合并结果
sim_pkpd_all <- do.call(rbind, sim_results)
sim_pkpd_all$dose <- factor(sim_pkpd_all$dose)

# 绘制PK和PD结果
p1 <- ggplot(sim_pkpd_all, aes(x=time, y=cp, group=interaction(id, dose), color=dose)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=dose), fun=mean, geom="line", linewidth=1) +
  labs(title="不同剂量的药物浓度", x="时间 (小时)", y="血药浓度", color="剂量 (mg)") +
  theme_bw()

p2 <- ggplot(sim_pkpd_all, aes(x=time, y=effect, group=interaction(id, dose), color=dose)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=dose), fun=mean, geom="line", linewidth=1) +
  labs(title="不同剂量的药物效应", x="时间 (小时)", y="效应", color="剂量 (mg)") +
  theme_bw()

gridExtra::grid.arrange(p1, p2, ncol=2)

3. 特殊人群仿真

模型仿真可以用于预测特殊人群（如肾功能不全患者）的药物行为：

  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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# 定义包含肾功能影响的模型
renal_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    beta_cl_crcl <- 0.75  # 肌酐清除率对药物清除率的影响系数
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.err <- 0.5
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (CRCL/100)^beta_cl_crcl  # 肌酐清除率对清除率的影响
    v <- exp(tv + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    cp = central / v
    cp ~ add(add.err)
  })
}

# 假设我们已经拟合了模型
# fit_renal <- nlmixr(renal_model, renal_data, est="focei")

# 为了演示，我们手动设置参数
theta_renal <- c(
  tka = log(1.5),
  tcl = log(0.08),
  tv = log(0.5),
  beta_cl_crcl = 0.75
)

omega_renal <- matrix(0, 3, 3)
diag(omega_renal) <- c(0.36, 0.09, 0.01)

# 创建用于仿真的rxode2模型
sim_renal_model <- rxode2({
  # 参数
  ka <- exp(tka + eta.ka)
  cl <- exp(tcl + eta.cl) * (CRCL/100)^beta_cl_crcl
  v <- exp(tv + eta.v)
  
  # 微分方程
  d/dt(depot) = -ka * depot
  d/dt(central) = ka * depot - cl/v * central
  
  # 输出变量
  cp = central / v
  auc = central / v / cl  # 计算AUC
})

# 定义不同肾功能组
renal_groups <- list(
  normal = c(mean=100, sd=15),      # 正常肾功能
  mild = c(mean=70, sd=10),         # 轻度肾功能不全
  moderate = c(mean=40, sd=7.5),    # 中度肾功能不全
  severe = c(mean=20, sd=5)         # 重度肾功能不全
)

# 生成参数样本
nsub_per_group <- 50
set.seed(seed)
sim_params_renal <- list()

for (group in names(renal_groups)) {
  params <- list(
    tka = rep(theta_renal["tka"], nsub_per_group),
    tcl = rep(theta_renal["tcl"], nsub_per_group),
    tv = rep(theta_renal["tv"], nsub_per_group),
    beta_cl_crcl = rep(theta_renal["beta_cl_crcl"], nsub_per_group),
    eta.ka = rnorm(nsub_per_group, 0, sqrt(omega_renal[1,1])),
    eta.cl = rnorm(nsub_per_group, 0, sqrt(omega_renal[2,2])),
    eta.v = rnorm(nsub_per_group, 0, sqrt(omega_renal[3,3])),
    CRCL = rnorm(nsub_per_group, renal_groups[[group]]["mean"], renal_groups[[group]]["sd"])
  )
  
  # 确保CRCL不为负
  params$CRCL <- pmax(params$CRCL, 5)
  
  params$group <- rep(group, nsub_per_group)
  sim_params_renal[[group]] <- params
}

# 合并所有组的参数
sim_params_all <- list()
for (param in names(sim_params_renal[[1]])) {
  sim_params_all[[param]] <- unlist(lapply(sim_params_renal, `[[`, param))
}

# 创建剂量事件表
dose_regimen_renal <- et(timeUnits="hours", amtUnits="mg") %>%
  et(amt=100, cmt=1, ii=24, addl=6) %>%  # 100mg，每24小时一次，共7次
  et(0:168)  # 观测时间点，从0到168小时

# 执行仿真
sim_renal <- sim_renal_model$solve(sim_params_all, dose_regimen_renal)
sim_renal$group <- factor(sim_params_all$group, 
                         levels=names(renal_groups),
                         labels=c("正常", "轻度", "中度", "重度"))

# 绘制不同肾功能组的结果
p3 <- ggplot(sim_renal, aes(x=time, y=cp, group=interaction(id, group), color=group)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=group), fun=mean, geom="line", linewidth=1) +
  labs(title="不同肾功能组的药物浓度", x="时间 (小时)", y="血药浓度", color="肾功能") +
  theme_bw()

# 计算并比较AUC
auc_summary <- aggregate(auc ~ id + group, data=sim_renal, FUN=max)
p4 <- ggplot(auc_summary, aes(x=group, y=auc, fill=group)) +
  geom_boxplot() +
  labs(title="不同肾功能组的AUC比较", x="肾功能", y="AUC", fill="肾功能") +
  theme_bw()

gridExtra::grid.arrange(p3, p4, ncol=2)

临床试验设计优化

模型仿真可以用于优化临床试验设计，包括样本量估计、采样时间优化和剂量选择：

1. 采样时间优化

 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
# 定义用于仿真的PK模型
pk_model <- rxode2({
  # 参数
  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
  
  # 输出变量
  cp = central / v
})

# 设置参数值
theta_pk <- c(
  tka = log(1.5),
  tcl = log(0.08),
  tv = log(0.5)
)

omega_pk <- matrix(0, 3, 3)
diag(omega_pk) <- c(0.36, 0.09, 0.01)

# 生成参数样本
nsub <- 100
set.seed(seed)
sim_params_pk <- list(
  tka = theta_pk["tka"] + rnorm(nsub, 0, sqrt(omega_pk[1,1])),
  tcl = theta_pk["tcl"] + rnorm(nsub, 0, sqrt(omega_pk[2,2])),
  tv = theta_pk["tv"] + rnorm(nsub, 0, sqrt(omega_pk[3,3])),
  eta.ka = rep(0, nsub),
  eta.cl = rep(0, nsub),
  eta.v = rep(0, nsub)
)

# 创建剂量事件表
dose_regimen_pk <- et(timeUnits="hours", amtUnits="mg") %>%
  et(amt=100, cmt=1) %>%
  et(0:24)  # 每小时一个时间点

# 执行仿真
sim_pk <- pk_model$solve(sim_params_pk, dose_regimen_pk)

# 计算每个时间点的浓度变异系数（CV）
cv_by_time <- aggregate(cp ~ time, data=sim_pk, FUN=function(x) sd(x)/mean(x)*100)

# 绘制CV随时间的变化
p5 <- ggplot(cv_by_time, aes(x=time, y=cp)) +
  geom_line() +
  geom_point() +
  labs(title="浓度变异系数随时间的变化", x="时间 (小时)", y="变异系数 (%)") +
  theme_bw()

# 找出信息量最大的采样时间点（CV最大的点）
optimal_times <- cv_by_time$time[order(cv_by_time$cp, decreasing=TRUE)[1:6]]
cat("最优采样时间点：", optimal_times, "\n")

# 绘制浓度-时间曲线，标记最优采样时间点
p6 <- ggplot(sim_pk, aes(x=time, y=cp, group=id)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=1), fun=mean, geom="line", color="red", linewidth=1) +
  geom_vline(xintercept=optimal_times, linetype="dashed", color="blue") +
  labs(title="浓度-时间曲线与最优采样时间点", x="时间 (小时)", y="血药浓度") +
  theme_bw()

gridExtra::grid.arrange(p5, p6, ncol=2)

2. 剂量选择优化

 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
# 定义目标暴露范围
target_range <- c(5, 10)  # 目标Cmax范围

# 创建不同剂量方案
dose_levels <- seq(50, 200, by=25)
dose_regimens <- lapply(dose_levels, function(dose) {
  et(timeUnits="hours", amtUnits="mg") %>%
    et(amt=dose, cmt=1) %>%
    et(0:24)
})

# 执行不同剂量的仿真
sim_results_dose <- list()
for (i in 1:length(dose_levels)) {
  sim_results_dose[[i]] <- pk_model$solve(sim_params_pk, dose_regimens[[i]])
  sim_results_dose[[i]]$dose <- dose_levels[i]
}

# 合并结果
sim_dose_all <- do.call(rbind, sim_results_dose)
sim_dose_all$dose <- factor(sim_dose_all$dose)

# 计算每个剂量的Cmax
cmax_by_dose <- aggregate(cp ~ id + dose, data=sim_dose_all, FUN=max)

# 计算每个剂量达到目标范围的概率
prob_in_range <- aggregate(cp ~ dose, data=cmax_by_dose, 
                          FUN=function(x) mean(x >= target_range[1] & x <= target_range[2])*100)

# 绘制不同剂量的Cmax分布
p7 <- ggplot(cmax_by_dose, aes(x=dose, y=cp)) +
  geom_boxplot() +
  geom_hline(yintercept=target_range, linetype="dashed", color="red") +
  labs(title="不同剂量的Cmax分布", x="剂量 (mg)", y="Cmax") +
  theme_bw()

# 绘制达到目标范围的概率
p8 <- ggplot(prob_in_range, aes(x=dose, y=cp, group=1)) +
  geom_line() +
  geom_point() +
  labs(title="达到目标范围的概率", x="剂量 (mg)", y="概率 (%)") +
  theme_bw()

gridExtra::grid.arrange(p7, p8, ncol=2)

# 找出最优剂量（达到目标范围概率最高的剂量）
optimal_dose <- as.numeric(as.character(prob_in_range$dose[which.max(prob_in_range$cp)]))
cat("最优剂量：", optimal_dose, "mg\n")

群体模拟和临床试验模拟

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
 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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
# 定义临床试验设计
n_subjects <- 60  # 受试者数量
n_groups <- 3     # 剂量组数量
subjects_per_group <- n_subjects / n_groups

# 创建受试者信息
set.seed(seed)
subjects <- data.frame(
  ID = 1:n_subjects,
  GROUP = rep(1:n_groups, each=subjects_per_group),
  WT = rnorm(n_subjects, 70, 15),
  SEX = rbinom(n_subjects, 1, 0.5),
  AGE = rnorm(n_subjects, 40, 10)
)

# 根据组别分配剂量
subjects$DOSE <- 50 * subjects$GROUP

# 创建访视计划
visits <- data.frame(
  VISIT = 1:5,
  TIME = c(0, 24, 48, 72, 96)  # 每天一次访视
)

# 创建采样时间
sampling_times <- c(0, 0.5, 1, 2, 4, 8, 12, 24)

# 创建完整的试验设计
trial_design <- expand.grid(
  ID = subjects$ID,
  VISIT = visits$VISIT
)
trial_design <- merge(trial_design, subjects, by="ID")
trial_design <- merge(trial_design, visits, by="VISIT")

# 创建剂量事件表
dose_events <- data.frame()
for (id in unique(subjects$ID)) {
  subject_data <- subjects[subjects$ID == id, ]
  dose <- subject_data$DOSE
  
  # 每次访视给药
  for (visit in 1:4) {  # 最后一次访视不给药
    time <- visits$TIME[visit]
    dose_events <- rbind(dose_events, data.frame(
      ID = id,
      TIME = time,
      AMT = dose,
      CMT = 1,
      EVID = 1,
      GROUP = subject_data$GROUP,
      WT = subject_data$WT,
      SEX = subject_data$SEX,
      AGE = subject_data$AGE
    ))
  }
  
  # 添加采样时间
  for (visit in 1:5) {
    base_time <- visits$TIME[visit]
    for (st in sampling_times) {
      if (visit < 5 || st <= 24) {  # 最后一次访视只采样到24小时
        dose_events <- rbind(dose_events, data.frame(
          ID = id,
          TIME = base_time + st,
          AMT = 0,
          CMT = 1,
          EVID = 0,
          GROUP = subject_data$GROUP,
          WT = subject_data$WT,
          SEX = subject_data$SEX,
          AGE = subject_data$AGE
        ))
      }
    }
  }
}

# 按ID和时间排序
dose_events <- dose_events[order(dose_events$ID, dose_events$TIME), ]

# 创建rxode2事件表
trial_et <- et(dose_events)

# 定义包含协变量的模型
trial_model <- rxode2({
  # 参数
  CL = exp(TVCL + 0.75*log(WT/70) + 0.2*(SEX==1) + eta.cl)
  V = exp(TVV + 0.75*log(WT/70) + eta.v)
  KA = exp(TVKA + eta.ka)
  
  # 微分方程
  d/dt(depot) = -KA * depot
  d/dt(central) = KA * depot - CL/V * central
  
  # 输出变量
  cp = central / V
})

# 设置参数值
theta_trial <- c(
  TVCL = log(0.08),
  TVV = log(0.5),
  TVKA = log(1.5)
)

omega_trial <- matrix(0, 3, 3)
diag(omega_trial) <- c(0.09, 0.01, 0.36)

# 生成参数样本
sim_params_trial <- list(
  TVCL = rep(theta_trial["TVCL"], n_subjects),
  TVV = rep(theta_trial["TVV"], n_subjects),
  TVKA = rep(theta_trial["TVKA"], n_subjects),
  eta.cl = rnorm(n_subjects, 0, sqrt(omega_trial[1,1])),
  eta.v = rnorm(n_subjects, 0, sqrt(omega_trial[2,2])),
  eta.ka = rnorm(n_subjects, 0, sqrt(omega_trial[3,3]))
)

# 执行临床试验仿真
sim_trial <- trial_model$solve(sim_params_trial, trial_et)

# 添加残差误差
set.seed(seed)
sim_trial$DV <- sim_trial$cp * (1 + rnorm(nrow(sim_trial), 0, 0.1))  # 10%比例误差

# 将剂量组转换为因子
sim_trial$GROUP <- factor(sim_trial$GROUP, 
                         levels=1:3, 
                         labels=c("低剂量", "中剂量", "高剂量"))

# 绘制不同剂量组的浓度-时间曲线
p9 <- ggplot(sim_trial[sim_trial$EVID == 0, ], 
            aes(x=TIME, y=DV, group=interaction(ID, GROUP), color=GROUP)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=GROUP), fun=mean, geom="line", linewidth=1) +
  labs(title="临床试验模拟：不同剂量组的浓度-时间曲线", 
       x="时间 (小时)", y="血药浓度", color="剂量组") +
  theme_bw()

# 计算每个受试者的AUC
auc_by_subject <- aggregate(cp ~ ID + GROUP, data=sim_trial, FUN=function(x) sum(x) * 24/length(x))

# 绘制不同剂量组的AUC比较
p10 <- ggplot(auc_by_subject, aes(x=GROUP, y=cp, fill=GROUP)) +
  geom_boxplot() +
  labs(title="临床试验模拟：不同剂量组的AUC比较", 
       x="剂量组", y="AUC", fill="剂量组") +
  theme_bw()

gridExtra::grid.arrange(p9, p10, ncol=2)

模型仿真的最佳实践

1. 仿真设计

明确仿真目标：确定仿真的具体目的（如剂量选择、采样优化等）
选择合适的模型：使用经过充分验证的模型进行仿真
考虑参数不确定性：在仿真中纳入参数估计的不确定性
设计合理的场景：设计能够回答研究问题的仿真场景

2. 仿真执行

使用足够的样本量：通常需要数百或数千个虚拟受试者
设置随机数种子：确保结果可重现
考虑计算效率：对于复杂模型，可能需要优化计算策略
验证仿真代码：在执行大规模仿真前验证代码的正确性

3. 结果分析

使用适当的统计方法：根据仿真目标选择合适的统计分析方法
可视化结果：使用图形直观地展示仿真结果
量化不确定性：报告结果的置信区间或预测区间
进行敏感性分析：评估关键假设对结果的影响

4. 报告和解释

清晰描述方法：详细说明模型、参数和仿真设计
透明报告结果：报告所有相关结果，不仅是支持结论的结果
讨论局限性：明确指出仿真的假设和局限性
提供实际建议：基于仿真结果提出具体的建议

实际应用案例：个体化给药方案设计

下面是一个使用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
 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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
# 定义包含协变量的PK模型
ind_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.08)
    tv <- log(0.5)
    beta_cl_wt <- 0.75
    beta_cl_age <- -0.2
    beta_cl_crcl <- 0.5
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.err <- 0.1
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt * (AGE/40)^beta_cl_age * (CRCL/100)^beta_cl_crcl
    v <- exp(tv + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    cp = central / v
    cp ~ prop(prop.err)
  })
}

# 假设我们已经拟合了模型
# fit_ind <- nlmixr(ind_model, ind_data, est="focei")

# 为了演示，我们手动设置参数
theta_ind <- c(
  tka = log(1.5),
  tcl = log(0.08),
  tv = log(0.5),
  beta_cl_wt = 0.75,
  beta_cl_age = -0.2,
  beta_cl_crcl = 0.5
)

omega_ind <- matrix(0, 3, 3)
diag(omega_ind) <- c(0.36, 0.09, 0.01)

# 创建用于仿真的rxode2模型
sim_ind_model <- rxode2({
  # 参数
  ka <- exp(tka + eta.ka)
  cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt * (AGE/40)^beta_cl_age * (CRCL/100)^beta_cl_crcl
  v <- exp(tv + eta.v)
  
  # 微分方程
  d/dt(depot) = -ka * depot
  d/dt(central) = ka * depot - cl/v * central
  
  # 输出变量
  cp = central / v
})

# 定义目标治疗窗
target_min <- 5  # 最小有效浓度
target_max <- 15  # 最大安全浓度

# 定义患者特征
patient_characteristics <- data.frame(
  ID = 1:4,
  WT = c(60, 80, 70, 90),
  AGE = c(30, 50, 70, 40),
  CRCL = c(120, 90, 50, 100),
  Description = c("年轻轻体重", "中年标准体重", "老年肾功能不全", "中年超重")
)

# 创建参数样本
sim_params_ind <- list()
for (i in 1:nrow(patient_characteristics)) {
  params <- list(
    tka = theta_ind["tka"],
    tcl = theta_ind["tcl"],
    tv = theta_ind["tv"],
    beta_cl_wt = theta_ind["beta_cl_wt"],
    beta_cl_age = theta_ind["beta_cl_age"],
    beta_cl_crcl = theta_ind["beta_cl_crcl"],
    eta.ka = 0,
    eta.cl = 0,
    eta.v = 0,
    WT = patient_characteristics$WT[i],
    AGE = patient_characteristics$AGE[i],
    CRCL = patient_characteristics$CRCL[i]
  )
  sim_params_ind[[i]] <- params
}

# 定义标准剂量方案
standard_dose <- 100  # mg
standard_regimen <- et(timeUnits="hours", amtUnits="mg") %>%
  et(amt=standard_dose, cmt=1, ii=24, addl=6) %>%
  et(0:168)

# 执行标准剂量方案的仿真
sim_results_standard <- list()
for (i in 1:length(sim_params_ind)) {
  sim_results_standard[[i]] <- sim_ind_model$solve(sim_params_ind[[i]], standard_regimen)
  sim_results_standard[[i]]$ID <- i
  sim_results_standard[[i]]$Description <- patient_characteristics$Description[i]
}

# 合并结果
sim_standard_all <- do.call(rbind, sim_results_standard)
sim_standard_all$Description <- factor(sim_standard_all$Description)

# 绘制标准剂量方案的结果
p11 <- ggplot(sim_standard_all, aes(x=time, y=cp, color=Description)) +
  geom_line() +
  geom_hline(yintercept=c(target_min, target_max), linetype="dashed", color="red") +
  labs(title="标准剂量方案 (100mg 每日一次)", 
       x="时间 (小时)", y="血药浓度", color="患者特征") +
  theme_bw()

# 根据患者特征调整剂量
adjusted_doses <- c(80, 100, 150, 80)  # 根据患者特征调整的剂量

# 创建个体化剂量方案
ind_regimens <- list()
for (i in 1:length(adjusted_doses)) {
  ind_regimens[[i]] <- et(timeUnits="hours", amtUnits="mg") %>%
    et(amt=adjusted_doses[i], cmt=1, ii=24, addl=6) %>%
    et(0:168)
}

# 执行个体化剂量方案的仿真
sim_results_ind <- list()
for (i in 1:length(sim_params_ind)) {
  sim_results_ind[[i]] <- sim_ind_model$solve(sim_params_ind[[i]], ind_regimens[[i]])
  sim_results_ind[[i]]$ID <- i
  sim_results_ind[[i]]$Description <- patient_characteristics$Description[i]
  sim_results_ind[[i]]$Dose <- adjusted_doses[i]
}

# 合并结果
sim_ind_all <- do.call(rbind, sim_results_ind)
sim_ind_all$Description <- factor(sim_ind_all$Description)

# 绘制个体化剂量方案的结果
p12 <- ggplot(sim_ind_all, aes(x=time, y=cp, color=Description)) +
  geom_line() +
  geom_hline(yintercept=c(target_min, target_max), linetype="dashed", color="red") +
  labs(title="个体化剂量方案", 
       x="时间 (小时)", y="血药浓度", color="患者特征") +
  theme_bw() +
  facet_wrap(~paste0(Description, " (", Dose, "mg)"), ncol=2)

gridExtra::grid.arrange(p11, p12, ncol=1)

# 计算在治疗窗内的时间百分比
calc_time_in_range <- function(data, min, max) {
  in_range <- data$cp >= min & data$cp <= max
  return(mean(in_range) * 100)
}

# 比较标准剂量和个体化剂量的效果
comparison <- data.frame(
  Patient = patient_characteristics$Description,
  Standard_Dose = sapply(split(sim_standard_all, sim_standard_all$ID), 
                        function(x) calc_time_in_range(x, target_min, target_max)),
  Individual_Dose = sapply(split(sim_ind_all, sim_ind_all$ID), 
                          function(x) calc_time_in_range(x, target_min, target_max)),
  Adjusted_Dose = adjusted_doses
)

# 打印比较结果
print(comparison)

总结

模型仿真是药代动力学/药效动力学建模的强大应用，可以帮助研究者优化临床试验设计、探索不同剂量方案、预测特殊人群的药物行为，以及支持监管决策。nlmixr2通过与rxode2的集成提供了全面的仿真功能，从基本的确定性仿真到复杂的临床试验模拟。

在实际应用中，模型仿真应该基于经过充分验证的模型，考虑参数不确定性，设计合理的仿真场景，并使用适当的统计方法分析结果。通过遵循最佳实践，模型仿真可以为药物开发提供宝贵的见解，减少试错成本，加速药物上市进程。

在下一节中，我们将通过一个完整的实际案例，展示如何使用nlmixr2进行从数据准备、模型构建、参数估计、模型诊断到模型应用的全过程，整合前面所学的所有知识。

nlmixr2实际案例分析：从数据准备到模型应用

本节将通过一个完整的实际案例，展示如何使用nlmixr2进行从数据准备、模型构建、参数估计、模型诊断到模型应用的全过程。我们将使用茶碱（theophylline）数据集作为示例，这是一个经典的药代动力学数据集，包含12名受试者在单次口服剂量后的血药浓度测量数据。

1. 数据准备与探索

首先，我们需要加载必要的包并准备数据：

 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
# 加载必要的包
library(nlmixr2)
library(ggplot2)
library(dplyr)
library(xpose.nlmixr2)
library(gridExtra)

# 加载茶碱数据集
data(theo_sd)
head(theo_sd)

# 数据结构概览
str(theo_sd)
summary(theo_sd)

# 检查数据完整性
any(is.na(theo_sd))  # 检查是否有缺失值

# 基本数据可视化
ggplot(theo_sd, aes(x=TIME, y=DV, group=ID, color=factor(ID))) +
  geom_point() +
  geom_line() +
  labs(title="茶碱血药浓度-时间曲线", 
       x="时间 (小时)", 
       y="血药浓度 (mg/L)",
       color="受试者ID") +
  theme_bw()

# 检查剂量信息
dose_info <- theo_sd %>% 
  filter(EVID == 1) %>% 
  select(ID, TIME, AMT)
print(dose_info)

# 计算基本PK参数（非房室分析）
nca_results <- data.frame(ID = unique(theo_sd$ID))

# 计算Cmax和Tmax
cmax_tmax <- theo_sd %>%
  filter(EVID == 0) %>%
  group_by(ID) %>%
  summarize(Cmax = max(DV),
            Tmax = TIME[which.max(DV)])

# 合并结果
nca_results <- merge(nca_results, cmax_tmax, by="ID")
print(nca_results)

通过初步数据探索，我们可以了解数据的基本结构和特点：

数据包含12名受试者
每名受试者接受单次口服剂量（320mg）
采样时间点为0、0.25、0.5、1、2、3.5、5、7、9、12、24小时
没有缺失数据

2. 模型构建

基于数据特点和药物特性，我们将构建一个一室模型来描述茶碱的药代动力学：

 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
# 定义一室模型
one_cmt_model <- function() {
  ini({
    # 固定效应参数（群体参数）的初始值
    tka <- log(1.5)    # 吸收速率常数的对数
    tcl <- log(0.04)   # 清除率的对数
    tv <- log(0.5)     # 分布容积的对数
    
    # 随机效应（个体间变异）
    eta.ka ~ 0.6       # 吸收速率常数的变异
    eta.cl ~ 0.3       # 清除率的变异
    eta.v ~ 0.1        # 分布容积的变异
    
    # 残差误差模型
    prop.err <- 0.1    # 比例误差
  })
  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
    
    # 计算血药浓度
    cp = central / v
    
    # 误差模型
    cp ~ prop(prop.err)
  })
}

这个模型包括：

一个吸收室（depot）和一个中央室（central）
一阶吸收过程（ka）
一阶消除过程（cl/v）
个体间变异（eta参数）
比例残差误差模型

3. 参数估计

接下来，我们使用不同的估计方法拟合模型，并比较结果：

 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
# 使用FOCEI方法拟合模型
fit_focei <- nlmixr(one_cmt_model, theo_sd, est="focei")

# 使用SAEM方法拟合模型
fit_saem <- nlmixr(one_cmt_model, theo_sd, est="saem")

# 查看FOCEI拟合结果
summary(fit_focei)

# 查看SAEM拟合结果
summary(fit_saem)

# 比较两种方法的参数估计结果
compare_params <- function(fit1, fit2, method1, method2) {
  params <- c("tka", "tcl", "tv", "prop.err")
  results <- data.frame(Parameter = params)
  
  results[[method1]] <- fixef(fit1)[params]
  results[[paste0(method1, "_RSE")]] <- 
    100 * sqrt(diag(fit1$cov.fixed)[params]) / abs(fixef(fit1)[params])
  
  results[[method2]] <- fixef(fit2)[params]
  results[[paste0(method2, "_RSE")]] <- 
    100 * sqrt(diag(fit2$cov.fixed)[params]) / abs(fixef(fit2)[params])
  
  return(results)
}

param_comparison <- compare_params(fit_focei, fit_saem, "FOCEI", "SAEM")
print(param_comparison)

# 比较目标函数值和信息准则
obj_comparison <- data.frame(
  Method = c("FOCEI", "SAEM"),
  OBJ = c(fit_focei$objf, fit_saem$objf),
  AIC = c(fit_focei$AIC, fit_saem$AIC),
  BIC = c(fit_focei$BIC, fit_saem$BIC)
)
print(obj_comparison)

# 选择最佳拟合模型（假设FOCEI结果更好）
best_fit <- fit_focei

通过比较不同估计方法的结果，我们可以选择最佳拟合模型。在这个例子中，我们假设FOCEI方法提供了更好的拟合结果。

4. 模型诊断

接下来，我们对选定的模型进行全面诊断：

 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
# 添加条件加权残差和归一化预测分布误差
best_fit <- addCwres(best_fit)
best_fit <- addNpde(best_fit)

# 基本诊断图
basic_gof <- plot(best_fit)
print(basic_gof)

# 残差分析
p1 <- plot(best_fit, type="cwres")
p2 <- plot(best_fit, type="npde")
grid.arrange(p1, p2, ncol=2)

# 个体拟合图
ind_plots <- plot(best_fit, type="individual", ncol=3, nrow=4)
print(ind_plots)

# 视觉预测检验（VPC）
vpc_result <- vpcPlot(best_fit, n=500, show=list(obs_dv=TRUE))
print(vpc_result)

# 使用xpose进行高级诊断
xp_fit <- xpose_data_nlmixr(best_fit)

p3 <- dv_vs_pred(xp_fit)
p4 <- dv_vs_ipred(xp_fit)
p5 <- res_vs_pred(xp_fit)
p6 <- res_vs_idv(xp_fit)

grid.arrange(p3, p4, p5, p6, ncol=2)

p7 <- res_distrib(xp_fit)
p8 <- res_qq(xp_fit)
p9 <- prm_distrib(xp_fit)
p10 <- ranef_qq(xp_fit)

grid.arrange(p7, p8, p9, p10, ncol=2)

# 统计检验
cat("CWRES正态性检验:\n")
print(shapiro.test(best_fit$CWRES))

cat("NPDE正态性检验:\n")
print(shapiro.test(best_fit$NPDE))

cat("NPDE均值检验（应接近0）:\n")
print(t.test(best_fit$NPDE))

cat("NPDE方差检验（应接近1）:\n")
print(var.test(best_fit$NPDE, rnorm(length(best_fit$NPDE))))

# 引导法分析（可选，计算时间较长）
# boot_fit <- bootplot(best_fit, nboot=100, seed=12345)
# print(boot_fit)
# plot(boot_fit)

通过全面的模型诊断，我们可以评估模型的适合度和预测能力：

观测值vs预测值图显示模型拟合良好
残差分析显示残差分布合理，无明显趋势
个体拟合图显示模型能够很好地描述个体数据
VPC显示模型预测与观测数据一致
统计检验确认残差分布合理

5. 模型改进

基于诊断结果，我们可以考虑改进模型。例如，我们可以尝试添加协变量效应：

 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
# 假设我们有体重数据（为演示目的，生成模拟数据）
set.seed(12345)
theo_sd$WT <- rnorm(nrow(theo_sd), 70, 15)
theo_sd$WT <- ave(theo_sd$WT, theo_sd$ID, FUN=function(x) rep(x[1], length(x)))

# 定义包含体重协变量的模型
cov_model <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(0.04)
    tv <- log(0.5)
    beta_cl_wt <- 0.75  # 体重对清除率的影响系数
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.err <- 0.1
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt  # 体重对清除率的影响
    v <- exp(tv + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    cp = central / v
    cp ~ prop(prop.err)
  })
}

# 拟合协变量模型
fit_cov <- nlmixr(cov_model, theo_sd, est="focei")

# 查看拟合结果
summary(fit_cov)

# 比较基础模型和协变量模型
model_comparison <- data.frame(
  Model = c("Base Model", "Covariate Model"),
  OBJ = c(best_fit$objf, fit_cov$objf),
  AIC = c(best_fit$AIC, fit_cov$AIC),
  BIC = c(best_fit$BIC, fit_cov$BIC),
  Parameters = c(length(fixef(best_fit)), length(fixef(fit_cov)))
)
print(model_comparison)

# 似然比检验
dof <- length(fixef(fit_cov)) - length(fixef(best_fit))
p_value <- 1 - pchisq(best_fit$objf - fit_cov$objf, df=dof)
cat("似然比检验 p值:", p_value, "\n")

# 如果协变量模型显著更好，则更新最佳模型
if (p_value < 0.05) {
  best_fit <- fit_cov
  cat("协变量模型显著优于基础模型\n")
} else {
  cat("协变量模型未显著改善拟合，保留基础模型\n")
}

通过模型比较，我们可以确定是否需要包含协变量效应。在这个例子中，我们假设体重对清除率有显著影响，因此选择协变量模型作为最终模型。

6. 模型应用：剂量优化

最后，我们使用最终模型进行应用，例如剂量优化：

  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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
# 创建rxode2模型用于仿真
library(rxode2)

# 从最终模型中提取参数
params <- fixef(best_fit)
omega <- best_fit$omega
sigma <- best_fit$sigma

# 创建仿真模型
if ("beta_cl_wt" %in% names(params)) {
  # 协变量模型
  sim_model <- rxode2({
    # 参数
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (WT/70)^beta_cl_wt
    v <- exp(tv + eta.v)
    
    # 微分方程
    d/dt(depot) = -ka * depot
    d/dt(central) = ka * depot - cl/v * central
    
    # 输出变量
    cp = central / v
  })
} else {
  # 基础模型
  sim_model <- rxode2({
    # 参数
    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
    
    # 输出变量
    cp = central / v
  })
}

# 定义目标治疗窗
target_min <- 5  # 最小有效浓度
target_max <- 15  # 最大安全浓度

# 创建不同剂量方案
dose_levels <- seq(200, 500, by=50)
dose_regimens <- lapply(dose_levels, function(dose) {
  et(timeUnits="hours", amtUnits="mg") %>%
    et(amt=dose, cmt=1, ii=12, addl=9) %>%  # 每12小时一次，共10次
    et(0:240)  # 观测时间点，从0到240小时
})

# 生成虚拟人群
nsub <- 100
set.seed(12345)

if ("beta_cl_wt" %in% names(params)) {
  # 协变量模型参数
  sim_params <- list(
    tka = rep(params["tka"], nsub),
    tcl = rep(params["tcl"], nsub),
    tv = rep(params["tv"], nsub),
    beta_cl_wt = rep(params["beta_cl_wt"], nsub),
    eta.ka = rnorm(nsub, 0, sqrt(omega[1,1])),
    eta.cl = rnorm(nsub, 0, sqrt(omega[2,2])),
    eta.v = rnorm(nsub, 0, sqrt(omega[3,3])),
    WT = rnorm(nsub, 70, 15)  # 生成体重数据
  )
} else {
  # 基础模型参数
  sim_params <- list(
    tka = rep(params["tka"], nsub),
    tcl = rep(params["tcl"], nsub),
    tv = rep(params["tv"], nsub),
    eta.ka = rnorm(nsub, 0, sqrt(omega[1,1])),
    eta.cl = rnorm(nsub, 0, sqrt(omega[2,2])),
    eta.v = rnorm(nsub, 0, sqrt(omega[3,3]))
  )
}

# 执行不同剂量的仿真
sim_results <- list()
for (i in 1:length(dose_levels)) {
  sim_results[[i]] <- sim_model$solve(sim_params, dose_regimens[[i]])
  sim_results[[i]]$dose <- dose_levels[i]
}

# 合并结果
sim_all <- do.call(rbind, sim_results)
sim_all$dose <- factor(sim_all$dose)

# 绘制不同剂量的浓度-时间曲线
p11 <- ggplot(sim_all, aes(x=time, y=cp, group=interaction(id, dose), color=dose)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=dose), fun=mean, geom="line", linewidth=1) +
  geom_hline(yintercept=c(target_min, target_max), linetype="dashed", color="red") +
  labs(title="不同剂量的茶碱浓度-时间曲线", 
       x="时间 (小时)", y="血药浓度 (mg/L)", color="剂量 (mg)") +
  theme_bw() +
  ylim(0, 30)

# 计算每个剂量的Cmin和Cmax（稳态）
steady_state <- sim_all %>%
  filter(time >= 120) %>%  # 只考虑稳态（假设5天后达到稳态）
  group_by(id, dose) %>%
  summarize(Cmin = min(cp),
            Cmax = max(cp),
            Cavg = mean(cp))

# 计算每个剂量达到目标范围的概率
prob_in_range <- steady_state %>%
  group_by(dose) %>%
  summarize(
    Prob_Cmin_Above_Min = mean(Cmin >= target_min) * 100,
    Prob_Cmax_Below_Max = mean(Cmax <= target_max) * 100,
    Prob_In_Range = mean(Cmin >= target_min & Cmax <= target_max) * 100
  )

# 绘制达到目标范围的概率
p12 <- ggplot(prob_in_range, aes(x=dose, y=Prob_In_Range, group=1)) +
  geom_line() +
  geom_point() +
  labs(title="不同剂量达到治疗窗的概率", 
       x="剂量 (mg)", y="概率 (%)") +
  theme_bw()

grid.arrange(p11, p12, ncol=2)

# 找出最优剂量（达到目标范围概率最高的剂量）
optimal_dose <- as.numeric(as.character(
  prob_in_range$dose[which.max(prob_in_range$Prob_In_Range)]))
cat("最优剂量：", optimal_dose, "mg 每12小时一次\n")

# 打印剂量优化结果表格
print(prob_in_range)

通过模型仿真和剂量优化，我们可以确定最佳给药方案，使血药浓度维持在目标治疗窗内的概率最大化。

7. 特殊人群剂量调整

最后，我们可以使用模型为特殊人群（如不同体重的患者）提供个体化剂量建议：

  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
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
# 定义不同体重组
weight_groups <- c(50, 70, 90)
weight_labels <- c("低体重 (50kg)", "标准体重 (70kg)", "高体重 (90kg)")

# 为每个体重组创建参数
sim_params_wt <- list()

for (i in 1:length(weight_groups)) {
  if ("beta_cl_wt" %in% names(params)) {
    # 协变量模型参数
    sim_params_wt[[i]] <- list(
      tka = rep(params["tka"], nsub),
      tcl = rep(params["tcl"], nsub),
      tv = rep(params["tv"], nsub),
      beta_cl_wt = rep(params["beta_cl_wt"], nsub),
      eta.ka = rnorm(nsub, 0, sqrt(omega[1,1])),
      eta.cl = rnorm(nsub, 0, sqrt(omega[2,2])),
      eta.v = rnorm(nsub, 0, sqrt(omega[3,3])),
      WT = rep(weight_groups[i], nsub)
    )
  } else {
    # 基础模型参数（体重不影响）
    sim_params_wt[[i]] <- list(
      tka = rep(params["tka"], nsub),
      tcl = rep(params["tcl"], nsub),
      tv = rep(params["tv"], nsub),
      eta.ka = rnorm(nsub, 0, sqrt(omega[1,1])),
      eta.cl = rnorm(nsub, 0, sqrt(omega[2,2])),
      eta.v = rnorm(nsub, 0, sqrt(omega[3,3]))
    )
  }
}

# 使用最优剂量创建给药方案
optimal_regimen <- et(timeUnits="hours", amtUnits="mg") %>%
  et(amt=optimal_dose, cmt=1, ii=12, addl=9) %>%
  et(0:240)

# 执行不同体重组的仿真
sim_results_wt <- list()
for (i in 1:length(weight_groups)) {
  sim_results_wt[[i]] <- sim_model$solve(sim_params_wt[[i]], optimal_regimen)
  sim_results_wt[[i]]$weight_group <- weight_labels[i]
}

# 合并结果
sim_all_wt <- do.call(rbind, sim_results_wt)
sim_all_wt$weight_group <- factor(sim_all_wt$weight_group, levels=weight_labels)

# 绘制不同体重组的浓度-时间曲线
p13 <- ggplot(sim_all_wt, aes(x=time, y=cp, group=interaction(id, weight_group), color=weight_group)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=weight_group), fun=mean, geom="line", linewidth=1) +
  geom_hline(yintercept=c(target_min, target_max), linetype="dashed", color="red") +
  labs(title=paste0("不同体重组使用相同剂量 (", optimal_dose, "mg) 的浓度-时间曲线"), 
       x="时间 (小时)", y="血药浓度 (mg/L)", color="体重组") +
  theme_bw() +
  ylim(0, 30)

# 计算每个体重组的Cmin和Cmax（稳态）
steady_state_wt <- sim_all_wt %>%
  filter(time >= 120) %>%
  group_by(id, weight_group) %>%
  summarize(Cmin = min(cp),
            Cmax = max(cp),
            Cavg = mean(cp))

# 计算每个体重组达到目标范围的概率
prob_in_range_wt <- steady_state_wt %>%
  group_by(weight_group) %>%
  summarize(
    Prob_Cmin_Above_Min = mean(Cmin >= target_min) * 100,
    Prob_Cmax_Below_Max = mean(Cmax <= target_max) * 100,
    Prob_In_Range = mean(Cmin >= target_min & Cmax <= target_max) * 100
  )

print(prob_in_range_wt)

# 为不同体重组计算调整剂量
if ("beta_cl_wt" %in% names(params)) {
  # 根据体重调整剂量
  adjusted_doses <- round(optimal_dose * (weight_groups/70)^params["beta_cl_wt"] / 50) * 50
} else {
  # 不调整剂量
  adjusted_doses <- rep(optimal_dose, length(weight_groups))
}

# 创建调整剂量的给药方案
adjusted_regimens <- list()
for (i in 1:length(adjusted_doses)) {
  adjusted_regimens[[i]] <- et(timeUnits="hours", amtUnits="mg") %>%
    et(amt=adjusted_doses[i], cmt=1, ii=12, addl=9) %>%
    et(0:240)
}

# 执行调整剂量的仿真
sim_results_adj <- list()
for (i in 1:length(weight_groups)) {
  sim_results_adj[[i]] <- sim_model$solve(sim_params_wt[[i]], adjusted_regimens[[i]])
  sim_results_adj[[i]]$weight_group <- weight_labels[i]
  sim_results_adj[[i]]$adjusted_dose <- adjusted_doses[i]
}

# 合并结果
sim_all_adj <- do.call(rbind, sim_results_adj)
sim_all_adj$weight_group <- factor(sim_all_adj$weight_group, levels=weight_labels)

# 绘制调整剂量后的浓度-时间曲线
p14 <- ggplot(sim_all_adj, aes(x=time, y=cp, group=interaction(id, weight_group), color=weight_group)) +
  geom_line(alpha=0.1) +
  stat_summary(aes(group=weight_group), fun=mean, geom="line", linewidth=1) +
  geom_hline(yintercept=c(target_min, target_max), linetype="dashed", color="red") +
  labs(title="不同体重组使用调整剂量的浓度-时间曲线", 
       x="时间 (小时)", y="血药浓度 (mg/L)", color="体重组") +
  theme_bw() +
  ylim(0, 30)

grid.arrange(p13, p14, ncol=2)

# 计算调整剂量后的Cmin和Cmax（稳态）
steady_state_adj <- sim_all_adj %>%
  filter(time >= 120) %>%
  group_by(id, weight_group, adjusted_dose) %>%
  summarize(Cmin = min(cp),
            Cmax = max(cp),
            Cavg = mean(cp))

# 计算调整剂量后达到目标范围的概率
prob_in_range_adj <- steady_state_adj %>%
  group_by(weight_group, adjusted_dose) %>%
  summarize(
    Prob_Cmin_Above_Min = mean(Cmin >= target_min) * 100,
    Prob_Cmax_Below_Max = mean(Cmax <= target_max) * 100,
    Prob_In_Range = mean(Cmin >= target_min & Cmax <= target_max) * 100
  )

# 打印调整剂量结果
print(prob_in_range_adj)

# 创建最终剂量建议表
dose_recommendations <- data.frame(
  Weight_Group = weight_labels,
  Weight = weight_groups,
  Standard_Dose = rep(optimal_dose, length(weight_groups)),
  Adjusted_Dose = adjusted_doses,
  Prob_In_Range_Standard = prob_in_range_wt$Prob_In_Range,
  Prob_In_Range_Adjusted = prob_in_range_adj$Prob_In_Range
)

print(dose_recommendations)

通过个体化剂量调整，我们可以为不同体重的患者提供定制的给药方案，提高治疗效果和安全性。

8. 结论与讨论

本案例展示了使用nlmixr2进行完整药代动力学分析的流程，包括：

数据准备与探索：检查数据结构和完整性，进行初步可视化
模型构建：基于药物特性和数据特点构建适当的模型
参数估计：使用不同方法估计模型参数，并选择最佳方法
模型诊断：全面评估模型适合度和预测能力
模型改进：通过添加协变量效应等方式改进模型
模型应用：使用模型进行剂量优化和个体化给药方案设计

通过这个案例，我们可以看到nlmixr2在药代动力学建模中的强大功能和灵活性。它不仅提供了多种参数估计方法，还集成了丰富的诊断工具和仿真功能，使研究者能够从数据中获取最大价值，并将模型应用于实际临床决策。

在实际应用中，模型的选择和评估应该基于多种因素，包括统计标准（如AIC、BIC）、诊断结果、生理合理性和实际应用需求。最终目标是建立一个既能准确描述数据又具有良好预测能力的模型，为药物开发和临床应用提供科学依据。

9. 参考文献

Fidler, M., Wilkins, J. J., Hooijmaijers, R., Post, T. M., Schoemaker, R., Trame, M. N., … & Wang, W. (2019). Nonlinear mixed-effects model development and simulation using nlmixr and related R open-source packages. CPT: Pharmacometrics & Systems Pharmacology, 8(9), 621-633.
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. Springer Science & Business Media.
Sheiner, L. B., & Beal, S. L. (1981). Evaluation of methods for estimating population pharmacokinetic parameters. II. Biexponential model and experimental pharmacokinetic data. Journal of Pharmacokinetics and Biopharmaceutics, 9(5), 635-651.
Upton, R. N., & Mould, D. R. (2014). Basic concepts in population modeling, simulation, and model-based drug development: part 3—introduction to pharmacodynamic modeling methods. CPT: Pharmacometrics & Systems Pharmacology, 3(1), 1-16.
Wang, W., Hallow, K. M., & James, D. A. (2016). A tutorial on RxODE: simulating differential equation pharmacometric models in R. CPT: Pharmacometrics & Systems Pharmacology, 5(1), 3-10.

总结与展望

nlmixr2的优势与局限性

通过本教程的学习，我们已经全面了解了nlmixr2在非线性混合效应建模中的应用。作为一个开源的R软件包，nlmixr2具有以下显著优势：

开源与可扩展性：作为开源软件，nlmixr2允许用户查看和修改源代码，根据特定需求进行定制。这种透明性和灵活性是商业软件所不具备的。
与R生态系统的无缝集成：nlmixr2可以直接利用R丰富的数据处理、统计分析和可视化功能，简化工作流程，提高效率。
多种估计方法：nlmixr2提供多种参数估计算法，包括FOCE、SAEM、FOCEI等，使用户可以根据数据特点和模型复杂性选择最合适的方法。
强大的仿真能力：通过与rxode2的集成，nlmixr2提供了高效的模型仿真功能，支持临床试验设计、剂量优化和特殊人群预测等应用。
活跃的社区支持：nlmixr2拥有活跃的开发者社区和用户群体，不断推出新功能和改进，并提供技术支持和问题解答。

然而，nlmixr2也存在一些局限性：

学习曲线：对于不熟悉R语言的用户，nlmixr2的学习曲线可能较陡峭，需要同时掌握R编程和药代动力学建模知识。
计算效率：在处理大规模数据集或复杂模型时，nlmixr2的计算效率可能不如一些专门优化的商业软件。
图形界面缺乏：与一些商业软件相比，nlmixr2缺乏直观的图形用户界面，主要依赖代码操作，可能不太适合偏好图形界面的用户。
文档和教程相对有限：虽然nlmixr2的文档在不断完善，但与成熟的商业软件相比，其教程和示例仍然相对有限。

未来发展趋势

随着药物开发和精准医疗的不断进步，非线性混合效应建模的重要性将继续增长。nlmixr2作为这一领域的重要工具，其未来发展趋势可能包括：

方法学创新：开发更高效、更稳健的参数估计算法，特别是针对复杂模型和稀疏数据的方法。
模型库扩展：建立更丰富的预定义模型库，涵盖常见的PK/PD模型结构，简化模型构建过程。
图形界面开发：开发更友好的图形用户界面，降低学习门槛，吸引更多非编程背景的用户。
与其他工具的集成：加强与其他药物开发工具的集成，如生理药代动力学（PBPK）模型、系统药理学模型等。
云计算支持：利用云计算资源加速计算密集型任务，提高大规模仿真和复杂模型拟合的效率。
机器学习结合：探索机器学习与传统NLME方法的结合，提高模型预测能力和个体化给药精度。

对读者的建议

对于希望在药代动力学/药效动力学建模领域发展的读者，我们提供以下建议：

打好基础：深入理解药代动力学/药效动力学的基本原理和非线性混合效应建模的统计基础。
实践为主：通过实际案例练习，从简单模型开始，逐步过渡到复杂模型，积累实战经验。
跨软件学习：不要局限于单一软件，尝试学习不同的建模工具（如NONMEM、Monolix、Phoenix NLME等），了解各自的优缺点。
关注方法学进展：定期阅读相关领域的最新研究文献，了解方法学创新和应用进展。
参与社区交流：加入nlmixr2用户社区，参与讨论，分享经验，共同进步。
跨学科合作：与临床医生、药学家、统计学家等专业人士合作，将建模结果应用于实际临床决策。

结语

非线性混合效应建模是连接药物开发和临床应用的桥梁，而nlmixr2作为一个强大的开源工具，为研究者提供了实现这一连接的有力支持。我们希望本教程能够帮助读者掌握nlmixr2的使用方法，并将其应用于自己的研究和实践中。

随着精准医疗时代的到来，个体化给药将成为临床实践的重要组成部分，而基于模型的方法将在其中发挥关键作用。通过不断学习和实践，我们可以更好地利用nlmixr2等工具，为患者提供更安全、更有效的药物治疗方案。

药代动力学/药效动力学建模是一门既需要科学知识又需要艺术直觉的学科，希望每位读者都能在这个领域找到自己的方向，并做出有意义的贡献。

祝您在非线性混合效应建模的探索之旅中取得成功！

参考文献

Almquist, J., Leander, J., & Jirstrand, M. (2015). Using sensitivity equations for computing gradients of the FOCE and FOCEI approximations to the population likelihood. Journal of Pharmacokinetics and Pharmacodynamics, 42(3), 191-209.
Bauer, R. J. (2019). NONMEM tutorial part I: Description of commands and options, with simple examples of population analysis. CPT: Pharmacometrics & Systems Pharmacology, 8(8), 525-537.
Beal, S. L., & Sheiner, L. B. (1982). Estimating population kinetics. Critical Reviews in Biomedical Engineering, 8(3), 195-222.
Comets, E., Brendel, K., & Mentré, F. (2008). Computing normalised prediction distribution errors to evaluate nonlinear mixed-effect models: the npde add-on package for R. Computer Methods and Programs in Biomedicine, 90(2), 154-166.
Duffull, S. B., Wright, D. F., & Winter, H. R. (2011). Interpreting population pharmacokinetic-pharmacodynamic analyses - a clinical viewpoint. British Journal of Clinical Pharmacology, 71(6), 807-814.
Fidler, M., Wilkins, J. J., Hooijmaijers, R., Post, T. M., Schoemaker, R., Trame, M. N., … & Wang, W. (2019). Nonlinear mixed-effects model development and simulation using nlmixr and related R open-source packages. CPT: Pharmacometrics & Systems Pharmacology, 8(9), 621-633.
Gastonguay, M. R. (2011). Full covariate models as an alternative to methods relying on statistical significance for inferences about covariate effects: a review of methodology and 42 case studies. PAGE 20, Abstr 2229.
Holford, N. (2005). The visual predictive check—superiority to standard diagnostic (Rorschach) plots. PAGE 14, Abstr 738.
Jonsson, E. N., & Karlsson, M. O. (1999). Xpose—an S-PLUS based population pharmacokinetic/pharmacodynamic model building aid for NONMEM. Computer Methods and Programs in Biomedicine, 58(1), 51-64.
Karlsson, M. O., & Savic, R. M. (2007). Diagnosing model diagnostics. Clinical Pharmacology & Therapeutics, 82(1), 17-20.
Kuhn, E., & Lavielle, M. (2005). Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics & Data Analysis, 49(4), 1020-1038.
Lavielle, M., & Mentré, F. (2007). Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software. Journal of Pharmacokinetics and Pharmacodynamics, 34(2), 229-249.
Nguyen, T. H., Mouksassi, M. S., Holford, N., Al-Huniti, N., Freedman, I., Hooker, A. C., … & Mentré, F. (2017). Model evaluation of continuous data pharmacometric models: metrics and graphics. CPT: Pharmacometrics & Systems Pharmacology, 6(2), 87-109.
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. Springer Science & Business Media.
Savic, R. M., & Karlsson, M. O. (2009). Importance of shrinkage in empirical Bayes estimates for diagnostics: problems and solutions. The AAPS Journal, 11(3), 558-569.
Sheiner, L. B., & Beal, S. L. (1981). Evaluation of methods for estimating population pharmacokinetic parameters. II. Biexponential model and experimental pharmacokinetic data. Journal of Pharmacokinetics and Biopharmaceutics, 9(5), 635-651.
Sheiner, L. B., Rosenberg, B., & Marathe, V. V. (1977). Estimation of population characteristics of pharmacokinetic parameters from routine clinical data. Journal of Pharmacokinetics and Biopharmaceutics, 5(5), 445-479.
Upton, R. N., & Mould, D. R. (2014). Basic concepts in population modeling, simulation, and model-based drug development: part 3—introduction to pharmacodynamic modeling methods. CPT: Pharmacometrics & Systems Pharmacology, 3(1), 1-16.
Wang, W., Hallow, K. M., & James, D. A. (2016). A tutorial on RxODE: simulating differential equation pharmacometric models in R. CPT: Pharmacometrics & Systems Pharmacology, 5(1), 3-10.
Zhang, L., Beal, S. L., & Sheiner, L. B. (2003). Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance. Journal of Pharmacokinetics and Pharmacodynamics, 30(6), 387-404.

引言#

nlmixr2安装和设置指南#

安装前的准备工作#

不同操作系统的安装步骤#

Windows系统#

macOS系统#

Linux系统（Ubuntu/Debian为例）#

验证安装#

安装开发版本（可选）#

使用R universe（推荐）#

使用remotes包#

安装支持包（可选）#

常见问题解决方案#

1. 编译器问题#

2. 依赖包安装失败#

3. 内存不足错误#

4. 路径问题#

5. 版本兼容性问题#

下一步#

非线性混合效应建模基本概念#

什么是非线性混合效应建模？#

非线性混合效应模型的数学表达#

常见的药代动力学模型#

一室模型#

二室模型#

口服给药模型#

药效动力学模型#

Emax模型（Hill方程）#

间接效应模型#

协变量模型#

nlmixr2中的模型表示#

非线性混合效应建模的优势#

非线性混合效应建模的应用#

实际应用示例：华法林PK/PD模型#

总结#

nlmixr2中的参数估计方法#

参数估计的基本原理#

nlmixr2中的估计方法概述#

FOCE和FOCEI方法#

原理#

nlmixr2中的使用#

控制参数#

优缺点#

SAEM方法#

原理#

nlmixr2中的使用#

控制参数#

优缺点#

nlme方法#

原理#

nlmixr2中的使用#

控制参数#

优缺点#

其他估计方法#

posthoc方法#

FO方法#

如何选择合适的估计方法#

实际应用示例#

参数估计的诊断和评估#

高级参数估计技巧#

1. 使用多种初始值#

2. 分步估计策略#

3. 使用不同的优化算法#

总结#

nlmixr2中的模型诊断技术#

模型诊断的目的#

nlmixr2中的基本诊断工具#

1. 模型拟合摘要#

2. 参数估计的不确定性#

3. 基本拟合图#

高级诊断技术#

1. 条件加权残差（CWRES）分析#

2. 归一化预测分布误差（NPDE）#

3. 视觉预测检验（VPC）#

4. 引导法（Bootstrap）#

5. 协变量筛选和模型比较#

使用xpose.nlmixr2进行高级诊断#

诊断结果的解释和模型改进#

1. 结构模型问题#

2. 统计模型问题#

引言

nlmixr2安装和设置指南

安装前的准备工作

不同操作系统的安装步骤

Windows系统

macOS系统

Linux系统（Ubuntu/Debian为例）

验证安装

安装开发版本（可选）

使用R universe（推荐）

使用remotes包

安装支持包（可选）

常见问题解决方案

1. 编译器问题

2. 依赖包安装失败

3. 内存不足错误

4. 路径问题

5. 版本兼容性问题

下一步

非线性混合效应建模基本概念

什么是非线性混合效应建模？

非线性混合效应模型的数学表达

常见的药代动力学模型

一室模型

二室模型

口服给药模型

药效动力学模型

Emax模型（Hill方程）

间接效应模型

协变量模型

nlmixr2中的模型表示

非线性混合效应建模的优势

非线性混合效应建模的应用

实际应用示例：华法林PK/PD模型

总结

nlmixr2中的参数估计方法

参数估计的基本原理

nlmixr2中的估计方法概述

FOCE和FOCEI方法

原理

nlmixr2中的使用

控制参数

优缺点

SAEM方法

原理

nlmixr2中的使用

控制参数

优缺点

nlme方法

原理

nlmixr2中的使用

控制参数

优缺点

其他估计方法

posthoc方法

FO方法

如何选择合适的估计方法

实际应用示例

参数估计的诊断和评估

高级参数估计技巧

1. 使用多种初始值

2. 分步估计策略

3. 使用不同的优化算法

总结

nlmixr2中的模型诊断技术

模型诊断的目的

nlmixr2中的基本诊断工具

1. 模型拟合摘要

2. 参数估计的不确定性

3. 基本拟合图

高级诊断技术

1. 条件加权残差（CWRES）分析

2. 归一化预测分布误差（NPDE）

3. 视觉预测检验（VPC）

4. 引导法（Bootstrap）

5. 协变量筛选和模型比较

使用xpose.nlmixr2进行高级诊断

诊断结果的解释和模型改进

1. 结构模型问题

2. 统计模型问题