最近有机会review了一下导师Andrew Hooker主导开发的群体药代动力学实验设计的一种方法（optimal design），就写出来与大家讨论，因为不是很正式的形式，所以行文会有些松散和随意，有问题我们可以评论区讨论，文中如果有缺漏我也会尽快修补。（ ref: Introduction to PopED )
在设计群体药代动实验，或者更一般的药代动力学实验时，我们大部分人都还是依据经验，好一些的会考虑到药物的半衰期等去调整一下实验方案，并没有一个公认的框架让我们去有意识优化实验方案。
首先是实验目的，通常来说，药代动力学实验的数据主要是用来1、进行统计检验；2、拟合模型得到参数。对于前一种目的来说，基于MCMC（蒙特卡洛马尔科夫链）的仿真模拟应用比较广泛，但对于后者，理论并不多。这篇文章主要介绍第二种，也是药代动力学实验更常见的一种类型。
对于以估算模型参数为导向的实验，很明显不同的设计可以对应不同的“信息量”，这里的信息量不是IT专业的词汇，而是指用不同实验方案得到的数据去估算模型参数时，参数估算结果的精度不同，更直接一点就是计算得到的参数RSE大小不同；很显然，能够使参数估算RSE尽可能小的实验方案更“优秀”一些，以下图1为例，我们直觉上就会在圆点标记的地方设置采血点（观察），因为那些点处在导数变化较剧烈的地方，携带了更多曲线的信息。
如何定量这种“直觉”就是这个optimal design的核心了，因为受众的关系（其实是不想打公式），我就用最简单直白的语言简单地介绍下如何定量实验设计在估算参数时的“信息量”的。在统计学中有一个Cramer-Rao不等式，粗略地说的就是类似于群体药代动力学模型这种过程估算的参数的不确定度（方差-协方差矩阵）有一个下限，它等于模型的Fisher information matrix（FIM)的逆矩阵。于是乎，FIM的逆矩阵就可以一定程度上地去反应PPK模型参数估算的精度了（实际上是一个FIM的函数），后面这个FIM的函数我就简称为J。这一段过于枯燥（其实是公式太多），我就点到为止，哪位有兴趣我们可以单聊。
根据上述的基于FIM的定量方法，我们可以得到实验设计和J的一一对应的关系，这里要注意是”实验设计“，也就是说你在实验得到数据之前就可以计算J的大小，这也是PopED优化实验设计的基础，因为一旦可以将实验设计和信息量函数J一一对应，就意味着我们可以优化实验设计找到那个J函数最小的实验设计方案，此即为优化。
另外，在实验设计时我们还需要考虑”经济“问题，实验经费永远都是有限的，比如像图2设置的采血点，虽然比图1多了不少，但实际上在J函数的值上并没有什么变化，也就是花大钱办了小事，这也是优化的方向。
以上说了这么多，回到实验设计的框架问题，首先是第一原则：以最经济的方案得到最适合拟合模型的数据，它通常包含以下要素：1、先验知识也就是实验预设的模型；2、分组数量；3、每组样本量；4、采血时间点；5、其他与实验相关的参数（比如剂量，给药间隔等等）。这就是药代动力学实验设计的框架，对里面各个要素的优化就是对实验设计的优化。
我们来一一分析以上的要素，首先是先验知识和预设模型，可能有人会疑惑，我就是准备做实验来建立模型的，怎么先就有了个模型。其实这里的模型指的不仅仅是复杂房室模型，可以非常简单直觉，比如你仅仅是觉得这观察值会随着时间变大也是一种模型，正比例模型嘛。这也引出了使用实验设计软件的一大注意事项：软件永远不会知道你不知道的事情。所以不要指望任何软件里会有一个大红按钮，你闭着眼睛一按就给你吐出一套优化的实验方案，这种事梦里才有。接下来是2到4点，这个但凡做过实验的都接触过，药代动实验不可缺少的部分。而第5点要说明一下，其他设计软件不清楚，但是对于接下来要讲解的PopED软件来说，这个相关参数需要分成2类：一个是连续型的比如剂量、给药间隔；另一个是离散的，比如给药次数之类的。因为PopED的优化引擎的关系，这两类参数是分开计算的，所以需要在设计之前就将需要优化的要素类型分析清楚。
综合以上，如果有一句话总结optimal design的原理，那就是将实验设计框架中的5个要素与一个反应数据信息量的函数J进行一一对应，并以此为基础对实验设计框架中的各个要素进行优化，找到能够使J函数达到极值的那个实验设计方案。
原理部分就到此为止，让我们来看看软件的实操。
这个软件包在R中是免费的，依然是由Andrew Hooker在维护，所以只需要在R中使用命令：install.package("PopED")就可以安装，然后library("PopED")加载。
PopED的使用主要分为3步：1、预设模型的搭建；2、参数的初始化和残差加入；3、PopED数据的生成。
以下我直接使用PopED的实例代码（可以在包内的帮助文件里找到），然后以注释形式解释它们的用途，但R语言相关的问题就不解释了，请自行搜索。
第一步：模型搭建，这个实例里用的是多次给药的口服吸收一房室模型例子，再次强调，你在设计实验之前必须对于这个实验有一个清晰的认识，并落实到模型的层面。如果实验结果得到新信息可以改善自己最初的模型，所以可以走文献初探→小规模实验验证设计→大规模实验的阶梯式实验设计方法。
ff <- function(model_switch,xt,parameters,poped.db){
#model_switch是用于帮助切换不同模型的变量；xt即为采血的时间点；parameters为预设
#模型的参数数组；poped.db是PopED数据，用于连接PopED内的各个函数。
with(as.list(parameters),{
y=xt
N = floor(xt/TAU)+1 #给药次数的计算
y=(DOSE*Favail/V)*(KA/(KA - CL/V)) *
(exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
#以上的一大串代码是该模型的解析解，其中DOSE为剂量、Favail为生物利用度、TAU为给 #药间隔，公式可以搜索得到。如果你的模型比较复杂是由一系列微分方程构成的，也可以 #用ODE的形式完成该函数
return(list( y=y,poped.db=poped.db))
})
}
函数建立完后，预设模型的骨架就搭好了。
第二步：参数初始化以及残差加入
sfg <- function(x,a,bpop,b,bocc){
#此处就是定义预设模型中的参数的步骤，这里要注意的是其中的x和a，我上文提到过在#PopED中，由于优化引擎的关系，与实验相关的参数如剂量、给药间隔、给药次数等会分 #为两类进行分析，一类是连续型的输入a里面，离散型的输入x里面，此例子中没有离散型 #的参数，故x为缺省状态；除此之外，bpop为群体药代动力学模型中的典型值；b为参数在#群体内分布的方差；bocc为参数的场合间变异。
#此函数内另外一点需要注意的就是，如果设计的实验不涉及群体效应，可以在此处缺省b的#输入项，此时即为普通药代动实验
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU=a[2])
return( parameters )
}
接下来是残差的加入，实际上如果残差模型简单，比如加和性和比例性再或者二者的混合，这种模型可以用PopED自带的内建函数即可，具体参阅PopED内的feps.add.prop等函数。
自定义残差的话，可以参考以下加和性和比例性混合的例子：
feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <-do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = y*(1+epsi[,1])+epsi[,2]
return(list( y= y,poped.db =poped.db ))
}
第三步就是生成PopED数据了，这个数据是准备工作的核心，也是准备工作的最后一步。
poped.db <- create.poped.database(
ff_fun=“ff”, #输入预设模型
fg_fun=“sfg”, #输入初始化的参数
fError_fun=“feps”, #输入残差模型，若简单残差可用feps.add.prop等内置模型
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
#输入初值，文献或者经验或者小规模实验 这里再次强调，做实验设计的是你而不是软件，#它不知道你不知道的东西，尝试在实验初期就充分的了解自己的实验，给出合理的初值
notfixed_bpop=c(1,1,1,0), #是否针对该参数优化实验，实验优化的目的是更好的估算参数，#但这个过程中具体优化哪几个参数的估算在此处进行设置，这里就不对Favail进行优化
d=c(V=0.09,KA=0.09,CL=0.25^2), #输入群体方差
sigma=c(0.04,5e-6), #输入模型输入残差方差
notfixedsigma=c(1,0), #与notfixed_bop类似
m=2, #样本分组数
groupsize=20, #每组样本量
xt=c( 1,2,8,240,245), #采集时间点
minxt=c(0,0,0,240,240), #采集点时间下限，优化时间的时候可以根据需要设置各个时间点#的上下限
maxxt=c(10,10,10,248,248), #采集点时间上限
bUseGrouped_xt=T, #两组使用是否使用同样的x采集点方案，同样也可以设置别的，具体 #可以参考PopED中bUseGrouped_系列的函数。
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), #其他实验参数，这里是剂量和给药间隔
maxa=c(DOSE=200,TAU=24), #其他实验参数上限
mina=c(DOSE=0,TAU=24))#其他实验参数下限
到此，实验设计的准备工作就做完了，可以使用PopED的各种功能了。
虚拟数据作图：
一般在实验之前，我们都会自己脑内生成一张目标的曲线图一类的东西，这个工作现在PopED可以帮助你完成了，使用的函数为:
plot_model_prediction(poped.db, model_num_points = 300)
此处为不带群体效应的曲线。
plot_model_prediction(poped.db, IPRED=T, DV=T, separate.groups=T, model_num_points = 300)
这是更复杂的带有群体效应和残差的图（不光滑的浅灰色区域标注的是残差效应），用这个图可以帮助你避免设置观察值负值或者低于定量下限的时间点。该函数的更多功能可以自行探索。
RSE预测：
这个属于PopED的核心功能了，就是计算该实验设计对应的J函数的值以及预测参数估算的RSE下限。这里提一句，NONMEM在7.5版本里也加入了该功能，在$DESIGN中，在得到了参数的final estimates后可以迅速计算uncertainty而不需要调用$COVARIANCE模块，省时省力，不足点是稳定性不够，不一定每次都能计算出结果，但不妨一试。以个人经验，如果$DESIGN失败的话，$COVARIANCE步骤同样失败的可能性超过一半，因此它也可以作为一个快速诊断工具使用。
RSE预测函数：evaluate_design(poped.db)
得到主要三部分结果，一个是基于J函数的OFV值，它和别的软件不同，OFV越大越好。第二是硬核玩家去使用的FIM矩阵。第三就是 大家比较关心的RSE了，它表示如果你按照这套实验方案去做实验，在设定初值的基础上，各个参数估算的RSE，也是大家比较关心的核心指标了。
实验优化：
最后这个函数就是PopED最重要的函数，用于优化实验参数，它会调用优化引擎来优化你设定好的实验参数。
函数：output <- poped_optim(poped.db, opt_xt=TRUE)
此处”opt_“为一系列输入，此处的xt为采集时间点优化，其它优化可以自行探索。这里需要说明一点，类似于每组人数、采集时间点总数之类的参数，在没有特殊状况的条件下都是多多益善的，所以通常不会优化出什么好的结果，有特殊需求的可以修改以上的函数去优化。另外就是，各个样本的采集时间点个数优化不包括在R的包中，需要去MATLAB平台的PopED使用。
结果如下同样分为三个部分：1、优化方案；2、OFV和效率比较；3、RSE比较。有两点需要说明一下，一个是优化方案，这个结果里面出现了2个10的采集点，这说明4个采集点和5个采集点的数据在估算参数上区别不大，PopED建议放弃一个。第二就是Efficiency，他比较的是两个方案的效率差异，因为考虑了模型中的参数数量，比OFV更准确一些。
特殊优化：
这里仅举两个例子：一个是采集时间点受限的情况，比如说采血点最好不要设置在凌晨或者饭点的时候，再或者仪器只能再有限的时间里开机等等，这里用的是输入一个采集时间点池，然后让PopED在其中选择最佳时间点。
poped.db.discrete <- create.poped.database(poped.db,discrete_xt =list(c(0:10,240:248))) output_discrete <- poped_optim(poped.db.discrete, opt_xt=TRUE)
discrete_xt=list(...)就是采集时间点池的输入过程，至于结果的解读，和上述函数一致。
另一个例子就是对其他实验参数的优化，比如剂量。因为实际中，如果经费不受限制，很多类似的量同样是多多益善无法优化的，所以这个优化过程需要自行定制，也是最考验大家实验设计功力的地方。
这里的例子提出的特殊要求是，在240小时的时间点，两组样本的谷浓度需要接近0.2和0.35，以此优化给药剂量。
自行定制的过程表现在自定义优化函数：
crit_fcn <- function(poped.db,...){
pred_df <- model_prediction(poped.db)
sum((pred_df[pred_df["Time"]==240,"PRED"] -c(0.2,0.35))^2)
#此处用的是最小平方和原理去接近目标谷浓度
}
然后再使用优化函数：
output_cost <- poped_optim(poped.db, opt_a = TRUE, opt_xt = FALSE, ofv_fun=crit_fcn, maximize = FALSE)
这里opt_xt和maximize都是false，表示这是在时间点优化完成的基础上再优化给药剂量，若不这么做可能优化引擎会出现问题。
以上只是PopED非常粗略的入门，希望它给你带来一些实验设计上的启发，然后还是强调一下，真正做实验设计的不是电脑而是你的人脑，你想的越细实验设计的会越完善，在此基础上PopED会是一个很趁手的工具。
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
