近期小编整理了一下生物等效性(BE)试验中PK分析部分的常规处理程序。于是就来分享一下这部分Winnonlin上的操作以及对应的利用SAS实现的程序。
BQL处理
在拿到样本检测数据后（浓度数据通常以Excel外部数据形式进行传输），将外部数据转化成SAS数据集，同时会对浓度数据中低于定量下限的BQL数据进行处理，根据方案中指定的规则进行BQL数据的替换，通常处理规则如下图。
BQL处理规则（摘自某方案）
上面截图为最常规的填补方式，也不乏有其他的花式填补。本文以最常规的填补方式为例。
winnonlin中的操作
在winnonlin中，选择数据选择>send to>Data Management>BQL。
选择BQL
在Main中设置SORT（受试者、周期/采血）、Time（采血时间点）、Concentration(浓度)、Carry（保留需要的变量）。
设置Main
设置填补规则，下图所示，设置Tmax前填补0，Tmax后填补ND。如果BQL无差别替换的话，将下图中的几个ND都替换成０即可。
设置替换规则
SAS中的代码
在SAS中，BQL的替换也是很简单的，下面来看看代码。
代码语言：javascript复制/****************************************************************************************************************************************
宏名称 : BQL
目的 : BQL填补
参数说明 :
inds 输入数据集
outds 输出数据集
con 浓度变量
time 时间变量(时间点,数值型)
treat 序列(T / R)
subjid 受试者编号
Keyword 关键词
bql1 填补方式1生成的变量
bql2 填补方式2生成的变量
YN 针对Tmax 之后出现的BQL
0：不考虑Tmax之后的连续BQL（统一按照缺失或者0处理）
1：考虑Tmax之后出现的连续BQL（连续出现的BQL之后的可测量结果填补ND，自身也ND）
Tmaxafter 限制Tmax之后出现的BQL
填补规则 :
第一种填补：将BQL全部替换成0
第二种填补：将Tmax之前的BQL替换成0，Tmax之后的BQL替换成ND
%bql(inds=adpc,outds=final,subjid=RANDNUM,
treat=TRTAL,con=AVAL,time=APCTIM,Keyword=BQL,bql1=BQL01,bql2=BQL02,YN=0);
________________________________________________________________________________________________________________________
__________________________________________________________________________________________________________________________
版本 日期 修改人 修改描述
--- ----------- ------------ ----------------------------------------------------------------------------------
1.0 2021.04.06 Setup 创建
****************************************************************************************************************************************/
%macro bql(inds=,outds=,subjid=,treat=,con=,time=,Keyword=BQL,bql1=AVAL01,bql2=AVAL02,YN=0,Tmaxafter=1);
data _temp_01;
set &inds.;
_temp_con=input(&con.,??best.);
run;
*计算Tmax时间点;
proc sql ;
create table _temp_02(where=(Cmax=_temp_con)) as
select distinct &subjid.,&treat.,&time.,_temp_con,max(_temp_con) as Cmax
from _temp_01 as a
group by &subjid.,&treat.;
create table _temp_03 as
select distinct a.*,Cmax,b.&time. as Tmax label="Tmax"
from _temp_01 as a
left join _temp_02 as b
on a.&subjid. =b.&subjid.
and a.&treat. =b.&treat. ;
quit;
proc sort data=_temp_03 out=_temp_03 sortseq=linguistic(numeric_collation=on);by &subjid. &treat. &time. ;quit;
data _temp_04;
set _temp_03;
by &subjid. &treat. &time. ;
retain missfl 0;
if first.&treat. then do;
*限制Tmax之后出现的连续的BQL;
if missing(_temp_con) %if &Tmaxafter. eq 1 %then %do;and &time.<=Tmax %end ; then missfl=0;
end;
else do;
if missing(_temp_con) %if &Tmaxafter. eq 1 %then %do;and &time.>Tmax %end ; then missfl=missfl 1;
%if &Tmaxafter. eq 1 %then %do;
else if missing(_temp_con) and &time.<=Tmax then missfl=0;
%end ;
else if ^missing(_temp_con) and missfl>=2 %if &Tmaxafter. eq 1 %then %do;and &time.>Tmax %end ; then missfl=missfl 1;
else if ^missing(_temp_con) and missfl<2 %if &Tmaxafter. eq 1 %then %do;and &time.>Tmax %end ; then missfl=0;
end;
run;
proc sort data=_temp_04 out=_temp_04 sortseq=linguistic(numeric_collation=on);by &subjid. &treat. descending &time. ;quit;
data _temp_04;
set _temp_04;
by &subjid. &treat. descending &time.;
missfllag1=lag(missfl);
if first.&treat. then call missing(missfllag1);
run;
proc sort data=_temp_04 out=_temp_05 sortseq=linguistic(numeric_collation=on);by &subjid. &treat. &time. ;quit;
data &outds.;
set _temp_05;
/*第一种填补方式：将所有的BQL当做0处理*/
&bql1.=&con.;
if index(upcase(&con.),"&Keyword.") then &bql1.="0";
/*考虑Tmax后连续2条BQL的情况--连续的俩条BQL之后的数据填补ND*/
%if &YN. eq 1 %then %do;
if missfl>=2 or (missfl=1 and missfllag1=2) then &bql1.="ND";
%end;
/*将Tmax之前的BQL替换成0，Tmax之后的BQL替换成ND*/
&bql2.=&con.;
if index(upcase(&con.),"&Keyword.") and &time.<Tmax then &bql2.="0";
else if index(upcase(&con.),"&Keyword.") and &time.>Tmax then &bql2.="ND";
%if &YN. eq 1 %then %do;
if missfl>=2 or (missfl=1 and missfllag1=2) then &bql2.="ND";
%end;
drop Cmax Tmax _temp_con missfllag1 missfl;
run;
proc delete data=work._temp_01 _temp_02 _temp_03 _temp_04 _temp_05 ;quit;
%mend;在BQL填补完成后，就可以在winnonlin中进行PK参数的计算，这里小编就不截图演示了，PK参数的计算小编也没有用SAS实现过，在winnonlin中点点点就可以了。
多因素方差分析
在生物等效性评价中，通常会对Cmax、AUClast、AUCinf_obs进行对数转化后，进行多因素方差分析。
方差分析样表
winnonlin中的操作
如下图，在计算出PK参数后，按照下图所示的点选，然后选择对应的的变量，就能方差分析、单向双侧T检验、置信区间都计算出来，关于Tmax非参数检验将下图所选的改成Crossover即可。这里就不展示winnonlin中具体的操作。
选择Bioequivalence
SAS中的代码
在生物等效性分析中，一般采用Proc Mixed过程对数据进行分析。此处以常规的两制剂、单次给药、双周期、双交叉试验为例。下面以Cmax为例进行计算。
代码语言：javascript复制data _rungrpds;
set adam.adpp;
if besfl="是" and UPCASE(param)="CMAX".;
_temp_aval=Log(aval);
run;
%let subject=RANDNUM;*受试者变量;
%let trtseq=trtseq;*给药顺序;
%let period=period;*周期;
%let trtal=trtal;*制剂;
ods exclude all;
ods output type3=_temp_01;
proc mixed data=_rungrpds method=type3;
class &subject. &trtseq. &period. &trtal.;
model _temp_aval=&trtseq. &period. &trtal. ;
random &subject.(&trtseq.);
Lsmeans &trtal.;
run;
ods exclude none;
data _final;
length _param _top $200.;
set _temp_01;
_param="ln(Cmax)";
if upcase(Source)=upcase("&TRTSEQ.") then do _top="序列";SEQ=1;end;
else if index(Source,')') then do _top="受试者(序列)";SEQ=2;end;
else if upcase(Source)=upcase("&TRTAL.") then do _top="制剂";SEQ=3;end;
else if upcase(Source)=upcase("&PERIOD.") then do _top="周期";SEQ=4;end;
else do ; _top="误差";SEQ=5;end;
label _top="试验因素" _param="药代动力参数" ;
keep _top _param DF SS MS FValue ProbF SEQ ;
run;
proc sort data=_final out= _final sortseq=linguistic(numeric_collation=on);by SEQ;quit;如上代码，即可实现Cmax相关的多因素方差分析。
输出结果(试验数据不展示了)
双向单侧T检验
下面来看一下第二个常规表格，双向单侧T检验，如下图所示。
双向单侧T检验样表
双向单侧T检验分析的实现也是较为简单的，采用PROC MIXED计算出参数，然后利用公式进行计算。采用此公式计算出来的结果与winnonlin结果数值是一致的（T2的正负除外）。
参考另外一篇文献，将T2公式稍改一下，即可得出与winnonlin一致的结果。
SAS中的代码
代码语言：javascript复制data _rungrpds;
set adam.adpp;
if besfl="是" and UPCASE(param)="CMAX".;
_temp_aval=Log(aval);
run;
%let subject=RANDNUM;*受试者变量;
%let trtseq=trtseq;*给药顺序;
%let period=period;*周期;
%let trtal=trtal;*制剂;
ods exclude all;
proc mixed data=_rungrpds;
class &subject. &trtseq. &period. &trtal.;
model _temp_aval=&trtseq. &period. &trtal. &subject.*&trtseq. /s ddfm=satterth ;
random intercept/subject=&subject. ;
lsmeans &trtal./diff=control('R') alpha=0.1 cl;
ods output diffs=_temp_diffs lsmeans=_temp_Lsmeans;
run;
ods exclude none;
proc sql UNDO_POLICY=NONE;
create table _temp_01 as
select distinct a.*,b.Estimate as Tmean,c.Estimate as Rmean
from _temp_diffs as a
left join _temp_Lsmeans(where=(&trtal.='T')) as b on 1=1
left join _temp_Lsmeans(where=(&trtal.='R')) as c on 1=1
;
quit;
data _temp_01;
set _temp_01;
/*Tost 值计算*/
t1_Tost=(input(Tmean,??best.)-input(Rmean,??best.)-log(0.80))/StdErr;
t1_TostP =1-probt(abs(t1_Tost),DF);
/*t2_Tost=(log(1.25)-(input(Tmean,??best.)-input(Rmean,??best.)))/StdErr;*书上公式，但算出来符号与winNonlin不一致;*/
t2_Tost=((input(Tmean,??best.)-input(Rmean,??best.))-log(1.25))/StdErr;
t2_TostP =1-probt(abs(t2_Tost),DF);
run;
置信区间法
常规表格三，生物等效性中的置信区间的计算，如下图所示。
置信区间样表
这部分的计算依然采用的是PROC MIXED过程步。这个表也是主要评价是否等效的一个关键表，下面来看一看这个表实现的过程。
SAS中的代码
代码语言：javascript复制
data _rungrpds;
set adam.adpp;
if besfl="是" and UPCASE(param)="CMAX";
_temp_aval=Log(aval);
run;
%let subject=RANDNUM;*受试者变量;
%let trtseq=trtseq;*给药顺序;
%let period=period;*周期;
%let trtal=trtal;*制剂;
ods exclude all;
proc mixed data=_rungrpds;
class &subject. &trtseq. &period. &trtal.;
model _temp_aval=&trtseq. &period. &trtal. ;
random &subject.(&trtseq.);
lsmeans &trtal./diff=control('R') alpha=0.1 cl;
ods output covparms=_temp_Covparms tests3=_temp_type diffs=_temp_diffs lsmeans=_temp_Lsmeans;
run;
ods exclude none;
/*T/R 几何均值*/
data _temp_01;
set _temp_Lsmeans;
Geomean=compress(put(round(exp(Estimate),0.001),12.3));
Mean=Estimate;
keep &trtal. Geomean Mean ;
run;
proc transpose data=_temp_01 out=_temp_01 prefix=AVAL;
by &trtal.;
var Geomean Mean;
run;
data _temp_01;
set _temp_01;
_nname_=strip(&trtal.)||strip(_NAME_);
run;
proc transpose data=_temp_01 out=_temp_01(drop=_name_) ;
var AVAL1;
id _nname_;
run;
/*T/R 例数*/
proc sql;
create table _temp_02 as
select distinct &trtal.,cats(count(distinct &subject. )) as n
from _rungrpds
group by &trtal.
order by &trtal. desc;
quit;
proc transpose data=_temp_02 out=_temp_02(drop=_name_) prefix=N;
var n;
id &trtal.;
run;
data _temp_03;
set _temp_diffs;
ratio=exp(Estimate);
Estimate2=exp(Estimate)*100;
Lower1=exp(Lower);
Upper1=exp(Upper);
keep Estimate2 ratio Lower1 Upper1 ;
run;
data _temp_Covparms;
set _temp_Covparms;
AVAL=sqrt(exp(Estimate)-1);
if CovParm="Residual" then do;n_name_="CV01";n_label_="个体内变异";end;
else do;n_name_="CV02";n_label_="个体间变异";end;
run;
proc transpose data=_temp_Covparms out=_temp_04 ;
var AVAL;
id n_name_ ;
idlabel n_label_ ;
run;
proc sql noprint;
select count(distinct &subject.) into :count trimmed from _rungrpds;
select ratio into :ratio trimmed from _temp_03;
select CV01 into :cv trimmed from _temp_04;
quit;
ods output output=_temp_05;
proc power;
twosamplemeans test=equiv_ratio
alpha = 0.05
lower = 0.8
upper = 1.25
meanratio = &ratio.
cv = &cv.
npergroup =&count.
power =.;
run;
ods exclude none;
data _final;
merge _temp_01-_temp_05;
run;非参数检验
针对Tmax的分析，通常是采用配对非参数的方法进行差异性检验（指导原则）。winnonlin上的crossover出来的结果应该不是采用这个方法做的，本来试着去还原一下winnonlin的方法，但是尝试几次，始终结果不一致，就放弃。
指导原则上的要求
SAS中的代码
代码语言：javascript复制%let subject=RANDNUM;*受试者变量;
%let trtseq=trtseq;*给药顺序;
%let period=period;*周期;
%let trtal=trtal;*制剂;
data _rungrpds;
set adam.adpp;
if besfl="是" and UPCASE(param)="TMAX";
run;
proc sort data=_rungrpds out=_rungrpds ;by &subject. ;;quit;
proc transpose data=_rungrpds out=Temp01 ;
by &subject. ;
var AVAL;
id &trtal. ;
run;
data Temp01;
set Temp01;
DIF=R-T;
run;
ods exclude none;
ODS TRACE ON;
ODS output TestsForLocation=temp2;
proc univariate DATA=Temp01;
var DIF;
run;
ODS TRACE OFF;
ods exclude all;
总结
上面就是生物等效性分析中的常规表格，还有几张常规表格，基本就类似于清单一样直接输出，当然还有几张图，譬如药时曲线等，作图可见小编历史推文SAS-时药曲线的绘制、SAS-时药曲线的绘制（完）、SAS-临床试验程序绘图合集（一），也都是很简单的操作。另外小编这里也开通一项功能，承接临床试验数据分析（包括但不限于PK参数计算及分析等）。
