专栏名称: 生信菜鸟团

生信菜鸟团荣誉归来，让所有想分析生物信息学数据的小伙伴找到归属，你值得拥有！

群体药代动力学-Nlmixr2学习

生信菜鸟团 · 公众号 · 生物 · 2024-12-22 21:32

正文

nlmixr文献附件中的示例和代码《Nonlinear Mixed‐Effects Model Development and Simulation Using nlmixr and Related R Open‐Source Packages》；

如果英文状态不太好理解，也可以看一下这篇毕业论文的详解《面向药动药效学的非线性混合自动机器学习模型研究》

欢迎交流~

背景

nlmixr是一个 免费的开源R包 ，用于拟合非线性药代动力学（PK）、药效学（PD）、联合PK-PD和定量系统药理学混合效应模型。目前，nlmixr能够拟合传统的房室PK模型以及使用常微分方程实现的更复杂的模型。作者开发这个软件包的主要目的是成为商业软件工具（如 NONMEM、Monolix 和 Phoenix NLME）的可靠免费替代品。

方法

本教程演示了四个步骤的模型构建原则：（i）如何使用nlmixr开发具有一阶吸收和线性消除的两室PK模型，（ii）如何使用nlmixr.xpose评估模型性能，（iii）如何在具有不确定性（如临床试验模拟中）和无不确定性（如VPC评估）的情况下模拟来自开发模型的新数据，（iv）使用shinyMixR遵循以项目为中心的工作流程进行模型开发和评估。

案例研究：一个简单的两室模型

为了说明nlmixr的工作原理，使用RxODE模拟一个简单的PK数据集，然后使用nlmixr进行拟合。讨论了nlmixr的模型语法、拟合模型以及查看和编辑优化结果。

数据结构

RxODE用于模拟40名受试者的该基础模型的数据，以及体重和性别对清除率和体重对中央分布容积的额外影响。 nlmixr 中的数据变量主要就是这几列比较重要，患者ID，时间，时间对应的DV观测值一般指血药浓度，AMT给药方式和剂量，EVID事件发生

下面生成测试数据 dat

rm(list = ls())
library(gridExtra)
library(rxode2)
library(MASS)
library(data.table)
library(ggplot2)
library(nlmixr2)
library(xpose)
library(xpose.nlmixr2)
library(shinyMixR)


#load('win-tutorial.Rdata')
## ----simulate_data, cache=T, echo=F, include=F---------------------------

mod <- RxODE({
  k10 = CL/V2
  k12 = Q/V2
  k21 = Q/V3
  d/dt(depot) =-KA*depot;
  d/dt(centr) = KA*depot - k10*centr - k12*centr + k21*peri;
  d/dt(peri)  =                        k12*centr - k21*peri;
  C2 = centr/V2;
  C3 = peri/V3;
  cp = C2 #*(1+cp.err)
})

theta <- c(TKA=1.05, TCL=0.121, TV2=1.939,
           TQ=0.282, TV3=5.65)

#The values for Cl1/F, Cl2/F, V1/F, V2/F and Ka were 121 ml/h, 282 ml/h, 1,939 ml, 5,650 ml and 1.05 h-1, respectively

omegaCor <- matrix(c(1,    0.5, 0.25, 0.1,  0,
                     0.5,  1,   0.5,  0.1,  0,
                     0.25, 0




    
.5, 1,    0.1,  0,
                     0.1,  0.1, 0.1,  1,    0,
                     0,    0,   0,    0,    1), dimnames=list(NULL,c("eta.CL","eta.V2","eta.V3", "eta.Q", "eta.KA")), nrow=5)

iiv.sd <- c(0.25, 0.25, 0.25, 0.3, 0.3) ## SDs of model parameters

iiv <- iiv.sd %*% t(iiv.sd)
omega <- iiv * omegaCor  # covariance matrix

sigma <- diag(1)*0.1
dimnames(sigma) <- list(NULL, c("cp.err"))

set.seed(740727)

mv <- mvrnorm(40, rep(0, dim(omega)[1]), omega) # Sample from covariance matrix

#Combine population parameters with IIV
params.all <-
  data.table(
    "ID" = seq(1:40),
    "CL" = theta['TCL'] * exp(mv[, 1]),
    "V2" = theta['TV2'] * exp(mv[, 2]),
    "V3" = theta['TV3'] * exp(mv[, 3]),
    "Q"  = theta['TQ']  * exp(mv[, 4]),
    "KA" = theta['TKA'] * exp(mv[, 5]),
    "WT" = round(rnorm(40,70,15)),
    "SEX" = rbinom(n = 40, prob = 0.5, size =1)
  )
#set the doses (looping through the 4 doses)
params.all[, AMT := 1200]

params.all$CL <- params.all$CL * (params.all$WT/70)^0.75
params.all$V2 <- params.all$V2 * (1 - 0.2 * params.all$SEX)


s = lapply(1:40, function(i) {
  #selects the parameters associated with the subject to be simulated
  params <- params.all[i]
  #creates an eventTable with 7 doses every 24 hours
  ev <- eventTable()
  ev$add.dosing(
    dose = params$AMT,
    nbr.doses = 28,
    dosing.to = 1,
    dosing.interval = 24,
    rate = NULL,
    start.time = 0
  )

  smp <- c(round(runif(1, 0, 1),3),
           round(runif(1, 1, 3),3),
           round(runif(1, 3, 6), 3),
           round(runif(1, 6, 12), 3),
           round(runif(1, 18




    
, 23.9), 3),
           round(runif(1, 168, 169),3),
           round(runif(1, 169, 171),3),
           round(runif(1, 171, 180),3),
           round(runif(1, 188, 191.9),3))

  ev$add.sampling(smp)

  x <- as.data.table(mod$run(params, ev, seed=740727))
  x$rv <- rnorm(nrow(x), 0, 0.075)
  x$DV <- round(x$cp * (1 + x$rv),1)
  x$ID <- i

  x[, names(params) := params]
})

sim <- as.data.table(do.call("rbind", s))

setnames(sim, "time", "TIME")

Dose <- expand.grid(TIME = seq(0, 7 * 24, 24), ID = params.all$ID, DV=0)
Dose <- data.table(merge(Dose, params.all, by = "ID"))
Dose[, EVID := 101]

sim[, EVID := 0]
sim[, AMT := 0]

# sim <- subset(sim, DV<10)
# sim <- subset(sim, DV>0.015)

# dilute by 50% at random
# sim <- sim[sample(nrow(sim), nrow(sim)*0.75), ]


sim <- sim[,c("ID","TIME","DV","WT","SEX","AMT","EVID")]
dat <- rbind(sim, Dose[,c("ID","TIME","DV","WT","SEX","AMT","EVID")])
setkey(dat, ID, TIME)

write.csv(dat, file="examplomycin.csv", quote = F, na = ".", row.names = F)

dat <- read.csv("examplomycin.csv", na.strings=c("."))

模型拟合基础

在金标准NONMEM的基础上开发的R包的nlmixr可以提供大多数在定量药理学和相关领域工作的分析师都能立即熟悉的模型规范语法。

整体模型结构

nlmixr的模型规范是使用包含 ini 和 model 的函数编写的，首先看一个非常简单的没有协变量的单房室模型。

ini模块

ini 指定初始条件，包括初始估计值，以及支持这些条件的算法的边界（目前，nlme和SAEM方法不支持参数边界）。 ini 类似于NONMEM中的 OMEGA和$SIGMA块。描述种群模型参数和残差模型参数的基本方法是使用以下表达式之一：

parameterName = initial estimate
parameterName = c(lower boundary, initial estimate)
parameterName = c(lower boundary, initial estimate, upper boundary)

与NONMEM不同，计算结果为数字的简单表达式可用于定义初始条件和边界，并且nlmixr允许指定无限参数：

parameterName = c(-Inf, log(70), 7)

估计与参数的典型总体值（类似于NONMEM中的“ETA”参数）的多变量正态个体偏差，以及这些偏差的方差/协方差矩阵（“OMEGA”矩阵）。这些个体间随机效应参数采用初始估计值，在nlmixr中使用波浪号（~）指定。在

指定吸收速率常数（k _a ）中可变性的初始估计值的示例是：

eta.ka ~ 0.1

加法运算符“+”可用于指定联合分布随机效应的方差/协方差矩阵，表达式的右侧以下三角矩阵形式指定初始估计值。CL 和 V 的示例如下所示：

eta.cl + eta.v ~ c(0.1,0.005, 0.1)

model

model 模块指定模型，类似于NONMEM中的 PRED和$ERROR模块。

模型是根据 ini 中定义的参数和数据中的协变量定义的。在指定微分方程或模型方程之前定义各个模型参数。例如，要根据population参数和individual参数定义单个间隙，可以使用：

Cl = exp(lCl+eta.Cl)

定义各个模型参数后，可以使用方程或RxODE代码直接定义模型。指定方程的RxODE方法基于使用 “d/dt” 符号编写微分方程的Leibnitz形式。要定义状态“depot”的ODE及其与“cent”的交互，可以使用以下方法定义状态值的浓度：

d/dt(depot) = -ka*depot
d/dt(cent) = ka*depot - Cl/V*cent
cp = cent/V

状态的初始条件可以通过定义state（0）来定义：

depot(0) = 0

残差模型可以使用加法误差、比例误差或加法和比例误差的组合来定义。对于加性误差，无法解释的变异性参数估计为标准差（SD），对于比例误差，估计为分数变异系数。add(parameter)；prop(parameter)；add(parameter1) + prop(parameter2)

在上面的示例中可以将中央室浓度的加性误差指定为：

cp ~ add(add.err)

模型可以用ODE表示，并且在单室、双室和三室线性PK模型变体的情况下用推注或输注给药表示为求解系统，可能与先前的吸收室相结合。使用求解的系统进行估计通常更快。指定求解的系统后， nlmixr根据 ini 模块中定义的参数名称推导区室模型的类型。要将整个ODE系统和加性残差模型折叠为一个简单的语句，可以编写以下内容：

linCmt() ~ add(add.err)

在本例中，浓度自动计算为中央室中的药物量除以中央体积。

拟合模型实操

使用模拟数据集开始尝试模型拟合，

run1：从nlme算法的单房室拟合开始：

## ----1cpt-nlme, cache=T, warning=F, message=F, results="hide", echo=F----

###这里是选择的加性残差，如果选择比例残差呢？比例参差的数值应该怎么设置呢？
#ka是吸收速率常数，k10是中心室的消除速率常数，k12是周边室的消除速率常数。
# tq 中央室到周边室的转移速率常数
model.1cpt.cf <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(1.5)
    tv  <- log(3)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v ~ 0.1
    add.err <- 0.01
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.err)
  })
}

fit.1cpt.cf.nlme <- nlmixr2(model.1cpt.cf, dat, est="nlme", table=tableControl(cwres=TRUE, npde=TRUE))

## ----1cpt-nlme-results, echo=F-------------------------------------------
print(fit.1cpt.cf.nlme)

该模型被过度拟合化，shrink%都红色99.9%，23.6%也略大，加性残差更是太大了。

run2：既然加性残差不合适，把加性换成比例残差，算法依旧是nlme，结果无法运行。尝试无果，放弃继续下一步。

run3：比例残差+saem算法的单房室

## ----1cpt-saem, cache=T, warning=F, message=F, results="hide", echo=F----

model.1cpt.ode <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(1.5




    
)
    tv  <- log(3)
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v ~ 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(center) = ka * depot - cl / v * center
    cp = center / v
    cp ~ prop(prop.err)
  })
}

fit.1cpt.ode.saem <- nlmixr2(model.1cpt.ode, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE), seed=740727)


## ----1cpt-saem-results, echo=F-------------------------------------------

print(fit.1cpt.ode.saem)

比例残差结果OK，SD参数还是太大，再换第三种算法

run4：比例残差+focei算法的单房室

## ----1cpt-focei, cache=T, warning=F, message=F, results="hide", echo=F----
model.1cpt.focei <- function() {
  ini({
    tka <- log(1.5)
    tcl <- log(1.5)
    tv  <- log(3)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v ~ 0.1
    prop.sd <- 0.01
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ prop(prop.sd)
  })
}

fit.1cpt.cf.focei <- nlmixr2(model.1cpt.focei, dat, est="focei", table=tableControl(cwres=TRUE, npde=TRUE))


## ----1cpt-focei-results, echo=F-------------------------------------------
print(fit.1cpt.cf.focei)

明显结果不好

run5：加性残差+saem算法的二房室

## ----2cpt, cache=T, warning=F, message=F, results="hide", echo=F---------

model.2cpt.ode <- function() {
  ini({
    tka <- log(1.05)
    tcl <- log(0.121)
    tv2  <- log(1.939)
    tv3  <- log(5.65)
    tq   <- log(0.282)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.q ~ 0.1
    add.err <- 75
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q  <- exp(tq + eta.q)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q/v3 * periph - q/v2 * center
    d/dt(periph) = q/v2 * center - q/v3 * periph
    cp = center / v2
    cp ~ add(add.err)
  })
}

fit.2cpt.ode.saem <- nlmixr2(model.2cpt.ode, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE), seed=740727)
print(fit.2cpt.ode.saem)

加性残差的结果一直拟合的不好，其他参数SD也不太行。SD一般得为正的偏小一点的数值。

run6：比例残差+saem算法的二房室

## ----2cptp, cache=T, warning=F, message=F, results="hide", echo=F--------




    


model.2cptp.ode <- function() {
  ini({
    tka <- log(1.05)
    tcl <- log(0.121)
    tv2  <- log(1.939)
    tv3  <- log(5.65)
    tq   <- log(0.282)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.q ~ 0.1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q  <- exp(tq + eta.q)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q/v3 * periph - q/v2 * center
    d/dt(periph) = q/v2 * center - q/v3 * periph
    cp = center / v2
    cp ~ prop(prop.err)
  })
}

fit.2cptp.ode.saem <- nlmixr2(model.2cptp.ode, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE), seed=740727)
print(fit.2cptp.ode.saem)

但SAEM算法的二房室结果也一般

run7：比例残差+saem算法的三房室

## ----3cpt, cache=T, warning=F, message=F, results="hide", echo=F---------

model.3cpt.ode <- function() {
  ini({
    tka <- log(1.42)
    tcl <- log(0.044)
    tv2  <- log(2)
    tv3  <- log(10)
    tv4  <- log(10)
    tq2   <- log(0.5)
    tq3   <- log(0.5)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.v4 ~ 0.1
    eta.q2 ~ 0.1
    eta.q3 ~ 0.1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    v4 <- exp(tv4 + eta.v4)
    q2  <- exp(tq2 + eta.q2)
    q3  <- exp(tq3 + eta.q3)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q2/v3 * periph1 - q2/v2 * center + q3/v4 * periph2 - q3/v2 * center
    d/dt(periph1) = q2/v2 * center - q2/v3 * periph1
    d/dt(periph2) = q3/v2 * center - q3/v4 * periph2
    cp = center / v2
    cp ~ prop(prop.err)
  })
}

fit.3cpt.ode.saem <- nlmixr2(model.3cpt.ode, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE), seed=740727)
print(fit.3cpt.ode.saem)

三房室结果也不太好

先总体看一下不加协变量时哪个模型相对较好，后续在这个基础上添加协变量。看起来目前结果相对而言较好一点的是第三个，OFV较小，AIC和BIC均较小，代表模型拟合的相较还行。于是后续选用这个模型（比例残差+SAEM）

这里发现了一个问题，虽然我设置了seed，每次运行结果都有一些差异。

run8：比例残差+saem算法的二房室+协变量体重

## ----2cptwt, cache=T, warning=F, message=F, results="hide", echo=F-------

dat$lnWT <- log(dat$WT/70)
###加上了体重的协变量，在ini中添加了wteff <- 0.35，在model中添加了cl = exp(tcl + wteff*lnWT + eta.cl)中的+ wteff*lnWT
model.2




    
cpt.ode.wtcl <- function() {
  ini({
    tka <- log(1.14)
    tcl <- log(0.0190)
    tv2  <- log(2.12)
    tv3  <- log(20.4)
    tq   <- log(0.383)
    wteff <- 0.35
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v2 ~ 1
    eta.v3 ~ 1
    eta.q ~ 1
    prop.err <- 0.075
  })
  model({
    ka = exp(tka + eta.ka)
    cl = exp(tcl + wteff*lnWT + eta.cl)
    v2 = exp(tv2 + eta.v2)
    v3 = exp(tv3 + eta.v3)
    q  = exp(tq + eta.q)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q/v3 * periph - q/v2 * center
    d/dt(periph) = q/v2 * center - q/v3 * periph
    cp = center / v2
    cp ~ prop(prop.err)
  })
}

fit.2cpt.ode.wtcl.saem <- nlmixr2(model.2cpt.ode.wtcl, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE))
print(fit.2cpt.ode.wtcl.saem)

run9：比例残差+saem算法的二房室+协变量性别

## ----2cptsex, cache=T, warning=F, message=F, results="hide", echo=F------
###加上了性别的协变量，在ini中添加了sexeff <- -0.2，在model中添加了v2 = exp(tv2 + sexeff*(SEX) + eta.v2)中的+ sexeff*(SEX)
#体重是在cl中加入的，性别是在v中加入的
model.2cpt.ode.sexv2 <- function() {
  ini({
    tka <- log(1.14)
    tcl <- log(0.0190)
    tv2  <- log(2.12)
    tv3  <- log(20.4)
    tq   <- log(0.383)
    sexeff <- -0.2
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v2 ~ 1
    eta.v3 ~ 1
    eta.q ~ 1
    prop.err <- 0.075
  })
  model({
    ka = exp(tka + eta.ka)
    cl = exp(tcl + eta.cl)
    v2 = exp(tv2 + sexeff*(SEX) + eta.v2)
    v3 = exp(tv3 + eta.v3)
    q  = exp(tq + eta.q)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q/v3 * periph - q/v2 * center
    d/dt(periph) = q/v2 * center - q/v3 * periph
    cp = center / v2
    cp ~ prop(prop.err)
  })
}

fit.2cpt.ode.sexv2.saem <- nlmixr2(model.2cpt.ode.sexv2, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE))
print(fit.2cpt.ode.sexv2.saem)

run10：比例残差+saem算法的二房室+协变量性别&体重

## ----2cptwtsex, cache=T, warning=F, message=F, results="hide", echo=F----
###同时加上了体重和性别的协变量，在cl和v中
model.2cpt.ode.wtcl.sexv2 <- function() {
  ini({
    tka <- log(1.14)
    tcl <- log(0.0190)
    tv2  <- log(2.12)
    tv3  <- log(20.4)
    tq   <- log(0.383)
    wteff  <- 0.35
    sexeff <- -0.2
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v2 ~ 1
    eta.v3 ~ 1
    eta.q ~ 1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + wteff*lnWT + eta.cl)
    v2 <- exp(tv2 + sexeff*(SEX) + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q  <- exp(tq + eta.q)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v2 * center + q/v3 * periph - q/v2 * center
    d/dt(periph) = q/v2 * center - q/v3 * periph
    cp = center / v2
    cp ~ prop(prop.err)
  })
}

## results
fit.2cpt.ode.wtcl.sexv2.saem <- nlmixr2(model.2cpt.ode.wtcl.sexv2, dat, est="saem", table=tableControl(cwres=TRUE, npde=TRUE))
print(fit.2cpt.ode.wtcl.sexv2.saem)

看一下各个模型的AIC和BIC