临床试验非劣效试验的样本量模拟计算SAS宏

这是张高魁, 夏结来, 姚晨《相对率的非劣效性试验检验效能及样本量的模拟计算方法及SAS实现》的程序
例: 一种足癣标准治疗方案, 疗程结束时的真菌学清除率为80%。现提出一种相对简单的治疗方案, 为考察对比新方案与标准治疗方案的疗效, 拟设计一个非劣效性试验, 临床提出新方案的疗效达到标准方案清除率的90% 作为临床非劣效的相对标准, 即: 当新方案与标准方案在疗程末的真菌学相对清除率( 新方案P标准方案) 的95% 可信区间下限不低于0. 9时, 即可判断其临床非劣效于标准方案。问题: ( 1) 、病例为100 对, 检验的效能是多少? ( 2) 、要使效能达到0. 8, 需要多少病例?
问题（1）解决的SAS宏如下，运行没有问题。

1. %macro RR(delta,n1,n2,distr1,distr2,essai,report=1);
   
2. %global pow;                        
   
3. data _null_;         
   
4.      call symput("seed ",date());
   
5.   run;
   
6.   Data depart;
   
7.     do numessai=1 to &essai;
   
8.        do pat=1 to %eval(&n1+&n2);
   
9.           if pat<=&n1 then do;
   
10.           Distr=rantbl(&seed,&distr1);
    
11.           group=1;
    
12.           end ;
    
13.         if pat> %eval(&n1) and pat <= %eval(&n1+&n2) then do;
    
14.            Distr=rantbl(&seed,&distr2);
    
15.            group=2;
    
16.         end;
    
17.       output;
    
18.       end;
    
19.      end;
    
20.    Run;
    
21. /**对每组随机数进行CMH检验，并计算相对率；*/
    
22. proc freq data=depart(where=(group in (1,2))) noprint;
    
23.    by numessai;
    
24.    tables group * distr /cmh2 scores=ridit;
    
25.    output out=TRT_AB N l_LGRRC1;
    
26. run;
    
27. /** 计算达到非劣效标准的次数；*/
    
28. Data power1;
    
29.    set TRT_AB;
    
30.    if l_lgrrc1^=.;
    
31.    power=l_lgrrc1>=δ
    
32.    scenario="TRT-A/TRT-B";
    
33. Run;
    
34. proc univariate data=power1 noprint;
    
35.   by scenario;
    
36.   var power;
    
37.   output out =power2 n=n sum=rep;
    
38. run;
    
39. /** 计算效能，并准备报告；*/
    
40. Data power3;
    
41.   length distr1 distr2 $ 30 scenario $ 20;
    
42.   merge power2;
    
43.   by scenario;
    
44.   lcl=100-100*Betainv(0.975,n-rep+1,rep)+0;
    
45.   ucl=100*Betainv(0.975,rep+1,n-rep)+0;
    
46.   lC="["|| compress (put(lcl,best4.))||"%;"|| compress (put(ucl,best4.))||"%"||"]";
    
47.   distr1=compress(symget("distr1")*100||"%");
    
48.   distr2=compress(symget("distr2")*100||"%");
    
49.   n1=symget("n1");
    
50.   n2=symget("n2");
    
51.   delta=symget("delta");
    
52.   simul=symget("essai");
    
53.   power=100*rep/n;
    
54.   call symput('pow',put(power,best4.1));
    
55. Run;
    
56. /**当报告开关为'开'时，报告模拟结果；*/
    
57. %if  %eval(&report)=1 %then %do;
    
58. Title2"以相对率为指标时检验效能的计算结果";
    
59. Proc report data= power3 headline headskip nowindows split ="@" ls=120;
    
60.   column("--" scenario simul delta power lc  ("Expected sample size" "--" n1 n2)  ("Expected rates" "--" distr1 distr2));
    
61. define scenario /order format=$18. left width=18 "Scenario";
    
62. define simul/display format=$5.center width=12 "N of @ simulations";
    
63. define delta /display format=$5.center width=6"Delta";
    
64. define power/display format=4.1 center width=6"Power";
    
65. define lC/display format=$12.center width=12"95% cl of @ Power";
    
66. define n1/display format=$10.center width=12 "TRT-A";
    
67. define n2/display format=$10.center width=12 "TRT-B";
    
68. define distr1/display format=$12.center width=12 "TRT-A";
    
69. define distr2/display format=$12.center width=12 "TRT-B";
    
70.    compute after;
    
71.      line @31 108*"-";
    
72.    endcomp;
    
73.    footnote "seed=&seed";
    
74. %end;
    
75. Run;
    
76. 
    
77. %mend RR;
    
78. 
    
79. %RR(0.9, 100, 100, 0.80, 0.80, 1000);

复制代码

Appreciate it.
我看了原文,方法比较容易理解。遗憾的是没有找到其中几个公式的出处。原来的程序，可能因为中英文乱码的原因，很难一一纠查。索性自己编写。根据作者的思路，我差不多照搬了第一个程序%RR。不同的是我直接用BINOMIAL数据，尤其对于尺寸大的数据，效率应该有所改善。并且抑制输出。原来的第二个程序%PSIZE看起来比较复杂些，没有细读。条条大路通罗马。我使用二分 算法，效率虽然不高，但是简单。
http://blog.sina.com.cn/s/blog_a3a926360101fp1q.html
京剧

现在问题是问题（2）的SAS宏，我安照文章编写，最后运行出错，程序如下：

1. %macro psize(delta,distr1,distr2,base=100,ratio=1,power=80,essai=1000,from=%STR(),to=%STR());
   
2.   options nodate nonumber ps=20;
   
3.   title;
   
4.   %if &from ne %then %goto fromto;
   
5.   %let n1=%eval(&base);
   
6.   %let n2=%eval(&n1*&ratio);
   
7.   %rr(&delta,&n1,&n2,&distr1,&distr2,&essai);
   
8.   %do %while(&pow < %eval(&power));
   
9.     %let step=%sysfunc(max(1,%sysfunc(int(%eval(&power)-&pow))));
   
10.     %let essai=%sysfunc(int(10000/%eval(&step)));
    
11.     %let n1=%eval(&n1+&step);
    
12.     %let n2=%eval(&n1+&ratio);
    
13.     %rr(&delta,&n1,&n2,&distr1,&distr2,&essai);
    
14.     %if &pow>= %eval(&power) %then %goto exit;
    
15.   %end;
    
16.   %do %while (&pow>%eval(&power));
    
17.     %let step=%sysfunc(max(1,%sysfunc(int(&pow-%eval(&power)))));
    
18.     %let essai= %sysfunc (int(10000/%eval (&step)));
    
19.     %let n1= %eval(&n1-&step);
    
20.     %let n2= %eval(&n1*&ratio);
    
21.     %rr(&delta,&n1, &n2,&distr1,&distr2,&essai);
    
22.     %if &pow <= %eval(&power) %then  %goto exit;
    
23.     %end;
    
24.   %fromto:
    
25.      %do t= %eval(&from-1) %to %eval(&to-1);
    
26.       %let n1=%eval(&t+1);
    
27.       %let n2=%eval(&n1*&ratio);
    
28.       %rr(&delta,&n1, &n2,&distr1,&distr2,&essai);
    
29.     %end;
    
30.   %exit:
    
31.   run;
    
32.   %mend psize;
    
33.   %psize (0.9,0.8,0.8);

复制代码

出错提示如下：

1. NOTE: 由调用宏“PSIZE”生成行。
   
2. 3       %rr(&delta,&n1,&n2,&distr1,&distr2,&essai);
   
3.         -
   
4.         180
   
5. WARNING: 没有解析宏 RR 的调用。
   
6. 
   
7. ERROR 180-322: 语句无效或未按正确顺序使用。
   
8. 
   
9. WARNING: 没有解析符号引用 POW。
   
10. ERROR: 在需要数值操作数的 %EVAL 函数或 %IF 条件中发现字符操作数。条件是: &pow < %eval(&power)
    
11. ERROR: %DO %WHILE 循环中的条件 &pow < %eval(&power) 生成无效或缺失值 ，宏将终止执行。
    
12. ERROR: 宏 PSIZE 将终止执行。

复制代码

哪哪位帮忙修正下？谢谢！！

我觉得换个思路：
解决作者提出的问题，如果用sas提供的过程power，差不多就5-10行程序，而且提供的信息更全面。
针对以上的宏：在参数未知的情况下，一味使用%EVAL（）来估计和判断，不出错误只是运气。
京剧

jingju11 发表于 2013-9-15 00:00
我觉得换个思路：
解决作者提出的问题，如果用sas提供的过程power，差不多就5-10行程序，而且提供的信息更 ...

PROC POWER? 愿闻其详？

把一楼，31行 改为

power=L_lgrrc1>=&delta;

后 一楼和二楼的宏都可以运行. 只是二楼的宏要手动停止， 还有一楼 25行处有一个 warning:
NOTE: Line generated by the invoked macro "RR".
5      distr /cmh2 scores=ridit;    output out=TRT_AB N L_LGRRC1; run;
                                                        --------
                                                        1

WARNING 1-322: Assuming the symbol LGRRC1 was misspelled as L_LGRRC1.

NOTE: There were 410160 observations read from the data set WORK.DEPART.
      WHERE group in (1, 2);

yongyitian 发表于 2013-9-15 11:41
把一楼，31行 改为

power=L_lgrrc1>=&delta;

谢谢回复。

其实是改成output out=TRT_AB N lgrrc1;  
(lgrrc1=Logit common relative risk, column 1 .)

正确的宏程序如下：

1. %macro RR(delta,n1,n2,distr1,distr2,essai,report=1);
   
2. %global pow;                        
   
3. data _null_;         
   
4.      call symput("seed ",date());
   
5.   run;
   
6.   Data depart;
   
7.     do numessai=1 to &essai;
   
8.        do pat=1 to %eval(&n1+&n2);
   
9.           if pat<=&n1 then do;
   
10.           Distr=rantbl(&seed,&distr1);
    
11.           group=1;
    
12.           end ;
    
13.         if pat> %eval(&n1) and pat <= %eval(&n1+&n2) then do;
    
14.            Distr=rantbl(&seed,&distr2);
    
15.            group=2;
    
16.         end;
    
17.       output;
    
18.       end;
    
19.      end;
    
20.    Run;
    
21. /**对每组随机数进行CMH检验，并计算相对率；*/
    
22. proc freq data=depart(where=(group in (1,2))) noprint;
    
23.    by numessai;
    
24.    tables group * distr /cmh2 scores=ridit;
    
25.    output out=TRT_AB N lgrrc1;
    
26. run;
    
27. /** 计算达到非劣效标准的次数；*/
    
28. Data power1;
    
29.    set TRT_AB;
    
30.    if l_lgrrc1^=.;
    
31.    power=l_lgrrc1>=δ
    
32.    scenario="TRT-A/TRT-B";
    
33. Run;
    
34. proc univariate data=power1 noprint;
    
35.   by scenario;
    
36.   var power;
    
37.   output out =power2 n=n sum=rep;
    
38. run;
    
39. /** 计算效能，并准备报告；*/
    
40. Data power3;
    
41.   length distr1 distr2 $ 30 scenario $ 20;
    
42.   merge power2;
    
43.   by scenario;
    
44.   lcl=100-100*Betainv(0.975,n-rep+1,rep)+0;
    
45.   ucl=100*Betainv(0.975,rep+1,n-rep)+0;
    
46.   lC="["|| compress (put(lcl,best4.))||"%;"|| compress (put(ucl,best4.))||"%"||"]";
    
47.   distr1=compress(symget("distr1")*100||"%");
    
48.   distr2=compress(symget("distr2")*100||"%");
    
49.   n1=symget("n1");
    
50.   n2=symget("n2");
    
51.   delta=symget("delta");
    
52.   simul=symget("essai");
    
53.   power=100*rep/n;
    
54.   call symput('pow',put(power,best4.1));
    
55. Run;
    
56. /**当报告开关为'开'时，报告模拟结果；*/
    
57. %if  %eval(&report)=1 %then %do;
    
58. Title2"以相对率为指标时检验效能的计算结果";
    
59. Proc report data= power3 headline headskip nowindows split ="@" ls=120;
    
60.   column("--" scenario simul delta power lc  ("Expected sample size" "--" n1 n2)  ("Expected rates" "--" distr1 distr2));
    
61. define scenario /order format=$18. left width=18 "Scenario";
    
62. define simul/display format=$5.center width=12 "N of @ simulations";
    
63. define delta /display format=$5.center width=6"Delta";
    
64. define power/display format=4.1 center width=6"Power";
    
65. define lC/display format=$12.center width=12"95% cl of @ Power";
    
66. define n1/display format=$10.center width=12 "TRT-A";
    
67. define n2/display format=$10.center width=12 "TRT-B";
    
68. define distr1/display format=$12.center width=12 "TRT-A";
    
69. define distr2/display format=$12.center width=12 "TRT-B";
    
70.    compute after;
    
71.      line @31 108*"-";
    
72.    endcomp;
    
73.    footnote "seed=&seed";
    
74. %end;
    
75. Run;
    
76. 
    
77. %mend RR;
    
78. 
    
79. %RR(0.9, 100, 100, 0.80, 0.80, 1000);
    
80. 
    
81. 
    
82. %macro psize(delta,distr1,distr2,base=100,ratio=1,power=80,essai=1000,from=%STR(),to=%STR());
    
83.   options nodate nonumber ps=20;
    
84.   title;
    
85.   %if &from ne %then %goto fromto;
    
86.   %let n1=%eval(&base);
    
87.   %let n2=%eval(&n1*&ratio);
    
88.   %rr(&delta,&n1,&n2,&distr1,&distr2,&essai);
    
89.   %do %while(&pow < %eval(&power));
    
90.     %let step=%sysfunc(max(1,%sysfunc(int(%eval(&power)-&pow))));
    
91.     %let essai=%sysfunc(int(10000/%eval(&step)));
    
92.     %let n1=%eval(&n1+&step);
    
93.     %let n2=%eval(&n1+&ratio);
    
94.     %rr(&delta,&n1,&n2,&distr1,&distr2,&essai);
    
95.     %if &pow>= %eval(&power) %then %goto exit;
    
96.   %end;
    
97.   %do %while (&pow>%eval(&power));
    
98.     %let step=%sysfunc(max(1,%sysfunc(int(&pow-%eval(&power)))));
    
99.     %let essai= %sysfunc (int(10000/%eval (&step)));
    
100.     %let n1= %eval(&n1-&step);
     
101.     %let n2= %eval(&n1*&ratio);
     
102.     %rr(&delta,&n1, &n2,&distr1,&distr2,&essai);
     
103.     %if &pow <= %eval(&power) %then  %goto exit;
     
104.     %end;
     
105.   %fromto:
     
106.      %do t= %eval(&from-1) %to %eval(&to-1);
     
107.       %let n1=%eval(&t+1);
     
108.       %let n2=%eval(&n1*&ratio);
     
109.       %rr(&delta,&n1, &n2,&distr1,&distr2,&essai);
     
110.     %end;
     
111.   %exit:
     
112.   run;
     
113.   %mend psize;
     
114.   %psize (0.9,0.8,0.8);

复制代码

其实，在实际工作中，此例子是是比较少见的，此例假设如下：

1. 一种足癣标准治疗方案,疗程结束时的真菌学清除率为80%。现提出一种相对简单的治疗方案,为考察对比新方案与标准治疗方案的疗效,拟设计一个非劣效性试验,临床提出新方案的疗效达到标准方案清除率的90%作为临床非劣效的相对标准,即:当新方案与标准方案在疗程末的真菌学相对清除率(新方案Π标准方案)的95%可信区间下限不低于0.9时,即可判断其临床非劣效于标准方案。问题:(1)、病例为100对,检验的效能是多少?(2)、要使效能达到0.8,需要多少病例?

复制代码

问题出现在“清除率的95%可信区间下限0.9”这里，在实际的临床试验中，不会知道这个下限的，仅知道非劣效界值，此值通常在5%~15%之间，不会超过15%，当然，界值约小需要样本量越多，比如某临床试验定了非劣效界值是15%，试验组及对照组的估计率是80%，这里的率95%可信区间下限不能等于80%-15%=65%。


如果本例中的“清除率的95%可信区间下限0.9”改为非劣效界值是15%，疗程结束时的真菌学清除率为80%，通过PASS计算得出每组样本量是112例。


当然，也可以通过公式计算，公式如下。

Rock2000 发表于 2013-9-15 18:07
其实，在实际工作中，此例子是是比较少见的，此例假设如下：

问题出现在“清除率的95%可信区间下限0.9” ...

楼上的大哥 ，PASS是买的？

webgu 发表于 2013-9-15 20:52
楼上的大哥 ，PASS是买的？

下载的 http://www.ncss.com/download/pass/updates/pass11/