RCT动态随机化，省下10万不香吗？
前言
随机分为随机抽样与随机分组，临床RCT说的都是随机分组(randomization)。随机化目的在于试验组间的非研究因素（covariate, cov）达到均衡。传统随机化方法都是事先制作好所有固定的随机号。动态随机化法(dynamic randomization, DR)，则是一边招募一边动态调整对象随机号。从随机化方法进化论的角度来说，是简单随机→整群随机→区组随机→分层区组随机→动态随机。前3种都没有考虑cov带来的影响，也就是当样本量不大的时候可能age sex等因素会组间不均衡。分层区组能解决cov均衡性问题，但是要求样本量大，比如age(老中青)sex(男女)，32组合，需要6倍普通的区组随机样本量。就多中心研究来说样本量也不是无穷尽的，且样本量增长会大大拖延RCT招募执行周期。动态随机法则可以很好解决小样本量cov不均衡问题，例如将100人随机分2组，最小化动态随机可以同时保障20多个cov均衡性。是不是很神奇？
原理
动态随机的常见方法有偏性掷币法、瓮法和最小化法等。偏性掷币法和瓮法主要是为了保障组间例数相近；而最小化随机分组的核心目的是保障分组后组间的某些影响因素是均衡的，应用最广泛。最小化法与1975年Pocock and Simon提出，三个参数确定病例分组：因素不平衡函数D，总体不平衡函数G，最优分配概率P。区分cov重要性的情况则要设定权重W。P=1，新病例直接分配到目标组，结果易于预测会导致选择偏倚。P=a ∈（0,1）新病例一定概率a分配到目标组，ICH指南推荐。最小化动态随机过程就是，根据前面已入组个体的情况进行动态分配随机号保证cov变量组间均衡。动态随机分组，研究对象被分入某组的概率不是固定不变的，是根据前面已入组个体的情况进行动态调整的。具体流程为：1、确定需要平衡的影响因素、各自的权重和目标组分配概率；2、第一个研究对象完全随机分组；3、从第二个研究对象开始，计算该研究对象被分组特定组后，两组间的预后因素差异大小；4、按照差异最小化的原则，将该研究对象按照分配概率进行随机分组。例如，某RCT最小化随机设计。treat(A、B)两种术式的近期疗效比较研究，cov(年龄、疾病分期和病理类型)是影响outcome(疗效好坏)的重要因素。需要所考虑的cov在treat组间均衡，才能获得无偏估计的treat→outcome效应。本例以第17个患者为例，说明最小化动态随机后应该入组A术式。本例D值计算由Pocock和Simon提出了四种计算方法，包括极差法、方差法、最大限值法和符号法，其中最常用的是极差法（Djk=|Xi1jk-Xi2jk|）和方差法（Djk=Var(Xi1jk,Xi2jk,…,XiNjk)），本例使用极差法计算不平衡函数D，cov权重设为123，最优分配概率P取0.8。，G(A)<G(B)，则A组为目标组。将新病例按概率0．8分配到A组，概率0.2分配到B组。
动态随机系统
大型临床试验采用中央随机化系统（IWRS），是一个基于多中心的受试者随机系统。通常包括受试者随机化和药品随机化，除了可以支持常规的静态随机算法（如：简单随机、区组随机、分层随机）外，还支持动态随机算法（如：最小化法、投币法、瓮法等）。其中云数医药、dacima、MACT等系统都包涵动态随机，但是作为小研究团队很难支付他们狮子大开口动辄5-10万的随机系统。其实有一些免费的动态随机软件，同样可以实现动态随机。例如按钮交互界面的MinimPy2、开源软件OxMAR，笔者测试下来体验不佳，但好在不用编程。笔者结合网络资料整理了R 与 SAS 代码。笔者认为，只要不是必须要求实时按秒计反馈分组信息，多中心RCT基于微信群或者邮件发送新患者信息，完全可以在1分钟内解决动态分组同步问题。省下10万不香吗，谁的钱也不是大风吹来的。再不济，用区组随机、分层区组随机都是可以的。
MinimPy2软件实现最小化随机
SAS代码
/* 编号：【ID】001 002 …
因素1：f1；记录日期：date（作为生成随机数的种子数）；因素2：f2；因素8：f8*/
/*@@@@@@@手动输入 新入组病人cov特征@@@@@@@*/
data patient;
input ID $ f1 Date mmddyy10. f2 f3 f4 f5 f6 f7 f8;
cards;
001 1 12-02-2017 1 1 2 1 2 1 1
002 1 01-11-2018 1 2 2 1 2 2 2
003 3 02-05-2018 1 1 1 1 2 1 1
004 1 03-02-2018 1 2 2 1 2 1 1
005 3 04-19-2018 1 2 1 2 1 1 2
006 1 05-08-2018 1 2 2 2 2 1 2
007 1 05-25-2018 1 2 1 2 2 1 1
008 1 06-15-2018 2 1 2 2 2 1 1
009 1 07-13-2018 1 1 2 1 2 2 1
010 1 09-24-2018 1 1 2 2 2 1 2
;
run;
data dataset;
set patient;
i+1;
run;
/*以上是录入数据，最新一行可视作新入组受试个体*/
proc sql; /*建立空表*/
create table test
(ID char(8), /*病人编号*/
f1 num(3),
f2 num(3),
f3 num(3),
f4 num(3),
f5 num(3),
f6 num(3),
f7 num(3),
f8 num(3),/*考虑的因素*/
date num(10) ,
group_ char(3), /*分组*/
DiffA num(3), /*假设新个体纳入A组时，两组在新个体相同水平上的差异*/
DiffB num(3), /*假设新个体纳入B组时，两组在新个体相同水平上的差异*/
P num(3)
);
quit;
%macro DR(id,f1,f2,f3,f4,f5,f6,f7,f8,date); /*输入编号，非处理因素f1-f8，入组时间*/
proc sql NOPRINT;
insert into test /*假设受试个体纳入A组时，A组与新个体在各因素的水平相同的频数*/
set id="&id" ,f1= %sysevalf(&f1+0),f2= %sysevalf(&f2+0),f3= %sysevalf(&f3+0),f4= %sysevalf(&f4+0),f5= %sysevalf(&f5+0),f6= %sysevalf(&f6+0),
f7= %sysevalf(&f7+0),f8= %sysevalf(&f8+0),date=%sysevalf(&date+0),group_='A';
select count(f1) into :AAcount_f1 from test where f1= %sysevalf(&f1) and group_='A';
select count(f2) into :AAcount_f2 from test where f2= %sysevalf(&f2) and group_='A';
select count(f3) into :AAcount_f3 from test where f3= %sysevalf(&f3) and group_='A';
select count(f4) into :AAcount_f4 from test where f4= %sysevalf(&f4) and group_='A';
select count(f5) into :AAcount_f5 from test where f5= %sysevalf(&f5) and group_='A';
select count(f6) into :AAcount_f6 from test where f6= %sysevalf(&f6) and group_='A';
select count(f7) into :AAcount_f7 from test where f7= %sysevalf(&f7) and group_='A';
select count(f8) into :AAcount_f8 from test where f8= %sysevalf(&f8) and group_='A';
/*假设受试个体纳入A组时，B组与新个体在各因素的水平相同的频数*/
select count(f1) into :ABcount_f1 from test where f1= %sysevalf(&f1) and group_='B';
select count(f2) into :ABcount_f2 from test where f2= %sysevalf(&f2) and group_='B';
select count(f3) into :ABcount_f3 from test where f3= %sysevalf(&f3) and group_='B';
select count(f4) into :ABcount_f4 from test where f4= %sysevalf(&f4) and group_='B';
select count(f5) into :ABcount_f5 from test where f5= %sysevalf(&f5) and group_='B';
select count(f6) into :ABcount_f6 from test where f6= %sysevalf(&f6) and group_='B';
select count(f7) into :ABcount_f7 from test where f7= %sysevalf(&f7) and group_='B';
select count(f8) into :ABcount_f8 from test where f8= %sysevalf(&f8) and group_='B';
update test
set group_= 'B' where id="&id"; /*假设受试个体纳入B组时，A组与新个体在各因素的水平相同的频数*/
select count(f1) into :BAcount_f1 from test where f1= %sysevalf(&f1) and group_='A';
select count(f2) into :BAcount_f2 from test where f2= %sysevalf(&f2) and group_='A';
select count(f3) into :BAcount_f3 from test where f3= %sysevalf(&f3) and group_='A';
select count(f4) into :BAcount_f4 from test where f4= %sysevalf(&f4) and group_='A';
select count(f5) into :BAcount_f5 from test where f5= %sysevalf(&f5) and group_='A';
select count(f6) into :BAcount_f6 from test where f6= %sysevalf(&f6) and group_='A';
select count(f7) into :BAcount_f7 from test where f7= %sysevalf(&f7) and group_='A';
select count(f8) into :BAcount_f8 from test where f8= %sysevalf(&f8) and group_='A';
/*假设受试个体纳入B组时，B组与新个体在各因素的水平相同的频数*/
select count(f1) into :BBcount_f1 from test where f1= %sysevalf(&f1) and group_='B';
select count(f2) into :BBcount_f2 from test where f2= %sysevalf(&f2) and group_='B';
select count(f3) into :BBcount_f3 from test where f3= %sysevalf(&f3) and group_='B';
select count(f4) into :BBcount_f4 from test where f4= %sysevalf(&f4) and group_='B';
select count(f5) into :BBcount_f5 from test where f5= %sysevalf(&f5) and group_='B';
select count(f6) into :BBcount_f6 from test where f6= %sysevalf(&f6) and group_='B';
select count(f7) into :BBcount_f7 from test where f7= %sysevalf(&f7) and group_='B';
select count(f8) into :BBcount_f8 from test where f8= %sysevalf(&f8) and group_='B';
quit;
proc sql noprint;
select count(*) into :count_A from test where id^="&id" and group_="A";
select count(*) into :count_B from test where id^="&id" and group_="B";
quit;
/*在Test数据集中入组对应信息，1对应A，2对应B(对照组)*/
data test;
set test;
format date DDMMYY10.;
if id="&id" then
do;
count_A=%sysevalf(&count_A+0); count_B=%sysevalf(&count_B+0);
date= %sysevalf(&date);
seed=%sysevalf(&date*&id);
/*假设受试个体纳入A组时，A组与新个体在各因素的水平相同的频数总和*/
DiffA=%sysevalf(%sysfunc(abs(&AAcount_f1-&ABcount_f1))+%sysfunc(abs(&AAcount_f2-&ABcount_f2))+%sysfunc(abs(&AAcount_f3-&ABcount_f3))
+%sysfunc(abs(&AAcount_f4-&ABcount_f4))+%sysfunc(abs(&AAcount_f5-&ABcount_f5))+%sysfunc(abs(&AAcount_f6-&ABcount_f6))
+%sysfunc(abs(&AAcount_f7-&ABcount_f7))+%sysfunc(abs(&AAcount_f8-&ABcount_f8)));
DiffB=%sysevalf(%sysfunc(abs(&BAcount_f1-&BBcount_f1))+%sysfunc(abs(&BAcount_f2-&BBcount_f2))+%sysfunc(abs(&BAcount_f3-&BBcount_f3))
+%sysfunc(abs(&BAcount_f4-&BBcount_f4))+%sysfunc(abs(&BAcount_f5-&BBcount_f5))+%sysfunc(abs(&BAcount_f6-&BBcount_f6))
+%sysfunc(abs(&BAcount_f7-&BBcount_f7))+%sysfunc(abs(&BAcount_f8-&BBcount_f8)));
/*假设受试个体纳入A组时，B组与新个体在各因素的水平相同的频数总和*/
if DiffA=DiffB then do; /*差异相等时，按相等概率进入A组或B组*/
P=uniform(%sysevalf(&date*&id));
if count_A=count_B then do;
if P<0.5 then group_='A' ; else group_='B';end;
if count_A<count_B then do;
if P<0.8 then group_='A' ; else group_='B';end;
if count_A>count_B then do;
if P<0.8 then group_='B' ; else group_='A';end;
end;
if DiffA<DiffB then /*A组差异较小时，按0.8概率进入A组，0.2概率进入B组*/
do;
P=uniform(%sysevalf(&date*&id));
if P<=0.8 then group_='A' ; else group_='B';
end;
if DiffA>DiffB then /*B组差异较小时，按0.8概率进入B组，0.2概率进入A组*/
do;
P=uniform(%sysevalf(&date*&id));
if P<=0.8 then group_='B' ; else group_='A';
end;
end;
run;
%mend;
/*@@@@@@@需要设置 nobs= ，第几个人@@@@@@@*/
%let nobs=10;
%macro test;
%do t=1 %to &nobs;
proc sql noprint;
select id,f1,f2,f3,f4,f5,f6,f7,f8,date into :id, :f1,:f2,:f3,:f4,:f5,:f6,:f7,:f8,:date from dataset where i=%sysevalf(&t);
quit;
%let id=&id; /*需要在宏test里面用%let赋值转换成局部宏参数，另一种方法用%global声明亦可*/
%let F_1=%sysevalf(&f1);
%let F_2=%sysevalf(&f2);
%let F_3=%sysevalf(&f3);
%let F_4=%sysevalf(&f4);
%let F_5=%sysevalf(&f5);
%let F_6=%sysevalf(&f6);
%let F_7=%sysevalf(&f7);
%let F_8=%sysevalf(&f8);
%let seed=%sysevalf(&date);
%DR(&id,&f_1,&f_2,&f_3,&f_4,&f_5,&f_6,&f_7,&f_8,&seed);
%end;
run;
%mend;
/*@@@@@@@直接运行@@@@@@*/
%test;
R代码
setwd("D:/写书/随机化/动态随机")
pacman::p_load(Minirand)
ntrt <- 2 #分组数AB
nsample <- 100 #计算好样本数
trtseq <- c(1, 2) #分组序列A=1,B=2
ratio <- c(1, 1) #分组比例1:1
library(readxl)
#第1个病人随机A或B
#定义f1-f8 labelcolnames(covmat1) = c("Gender", "Age", "Hypertension", "Use of Antibiotics")
covmat1 <- read_excel("8因素动态随机.xlsx", sheet = "1")
##ID f1 f2 f3 f4 f5 f6 f7 f8
#第2个病人入组A术式，Djk最小*/
covmat2 <- read_excel("8因素动态随机.xlsx", sheet = "2")
covmat2s <-covmat2[1:2,2:9]
#--------此处手动增加 第 n个病人--------------
covmat5 <- read_excel("8因素动态随机.xlsx", sheet = "5")
covmat5s <- covmat10[1:5,2:9]
covmat10 <- read_excel("8因素动态随机.xlsx", sheet = "10")
covmat10s <- covmat10[1:10,2:9]
covwt <- c(1/8, 1/8, 1/8, 1/8,1/8, 1/8, 1/8, 1/8) #equal weights;cov协变量权重一样
#假设多中心RCT预期样本量为100人
rm(res)
res <- rep(100, nsample)
#第一个病人按照seed，随机分A=1或B=2;最小化极差法动态随机化
set.seed(2021)
res[1] = sample(trtseq, 1, replace = TRUE, prob = ratio/sum(ratio))
# 第2人
for (j in 2:nsample){
res[j] <- Minirand(covmat=covmat2s, j, covwt=covwt, ratio=ratio,ntrt=ntrt, trtseq=trtseq, method="SD", result=res, p = 1)
}
trt1 <- res[1:2]
head(trt1,2)
#p=0.9 会产生问题，新患者导致前面分组结果有比例改变。90%概率 2 1；10%概率 2 2 。
# 第5人
for (j in 2:nsample){
res[j] <- Minirand(covmat=covmat5s, j, covwt=covwt, ratio=ratio,ntrt=ntrt, trtseq=trtseq, method="SD", result=res, p =1)
}
trt1 <- res[1:5]
head(trt1,5)
# p=0.9 会产生多种可能2 1 2 2 1；2 1 1 1 2
# 第10人
for (j in 2:nsample){
res[j] <- Minirand(covmat=covmat10s, j, covwt=covwt, ratio=ratio,ntrt=ntrt, trtseq=trtseq, method="SD", result=res, p = 1)
}
#trt1 为治疗分组AB, cov均衡但是未显示P值。
trt1 <- res[1:10]
head(trt1,10)
#2 1 1 2 2 1 1 2 2 1 与SAS 2 1 1 2 2 1 1 2 1 2 稍微有点不同
#显示分组频数
covmat10s <- as.matrix(covmat10s)
balance1 <- randbalance(trt1, covmat10s, ntrt, trtseq)
balance1
totimbal(trt = trt1, covmat = covmat10s, covwt = covwt, ratio = ratio, ntrt = ntrt, trtseq = trtseq, method = "SD")
references:
https://www.tandfonline.com/doi/abs/10.1080/03610918.2019.1619765?journalCode=lssp20 https://blog.csdn.net/weixin_41065383/article/details/83113624 赵丽君, 李宏田, 段蜚藩,等. 最小化随机分组方法介绍及其SAS实现[J]. 中国生育健康杂志, 2018, v.29(02):5-9. Yan Cui, Huaien Bu, Hongwu Wang, Shizhong Liao, "MACT: A Manageable Minimization Allocation System", Computational and Mathematical Methods in Medicine, vol. 2014, Article ID 645064, 8 pages, 2014. https://doi.org/10.1155/2014/645064 https://osdn.net/projects/minimpy2/downloads/73910/MinimPy2.zip/ https://sourceforge.net/projects/oxmar/files/ https://www.youtube.com/watch?v=D1zQRwdVPk4
原创不易，转载请说明来自本公众号。