SAS编程：按SOC和PT类别汇总AE的受试者发生率

临床试验项目中，安全性分析会对AE受试者发生率按试验组进行汇总，最近手动写了这类表QC侧程序，基于此捋一捋这类表的输出过程。

AE 受试者发生率

这篇文章针对表中各层级频数汇总，单独处理后汇总结果。在SAS编程-Table：层级拼接法输出AE SOC、PT的受试者发生率中，采用的是层级拼接法处理这类Table。

表中说明

统计量介绍

输出一张表，要先理解各个统计量的含义，程序代码的实现是理解含义之后自然产物。受试者AE发生率的Table涉及到４类统计量的计算：(对应图中数字标记)

1. 各试验组BigN的计算，即个发生率的分母；
2. 各试验组发生AE的受试者的频数及发生率；
3. 各试验组发生具体SOC的受试者的频数及发生率；
4. 各试验组发生具体PT的受试者的频数及发生率；

首先，需要明确的是，受试者的发生率是“人数比人数”，所以我们所有的计数都是基于“数人数”的前提。

对于第１类统计量BigN，人数就是简单地计数进入安全性分析集中的各组受试者，直观地数人头。剩下３类与第１类不同，分别计数发生AE的人数、发生具体SOC的AE的人数、发生具体PT的AE的人数。这会遇到一个受试者发生多个AE的情况，这时，受试者也只会被计数一次。在编程上，这里就需要对“重复受试者”去重。

其次，关于输出的排序，一般是按照汇总组或试验组的SOC、PT频数进行降序排序，即频数由大到小排序。

对于，频数相同的情况，如果TFLShell中没有明确要求，我们需要与统计师沟通确认。我做的这张表是先按试验组SOC、PT的频数降序排序后，再按SOC、PT的首字母排序。

变量介绍

CRF收集的不良事件名称，会保存的AETERM (Reported Term for the Adverse Event)变量中。根据外部字典(eg, MedDRA)对AETERM进行编码，标准化名称保存在AEDECOD (Dictionary-Derived Term)变量中，也就是我们常说的PT(Preferred Terrm)。

一个不良事件可能涉及多个身体系统器官，其中一个是主要的器官，主要器官信息就保存在AESOC (Primary System Organ Class)变量中。

在临床试验分析中，申办方可能特别关注AE对某个器官的影响，只需要将它纳入分析而不关注其他，这个信息保存在AEBODSYS (Body System or Organ Class)变量中。

纳入分析的系统器官并不一定是主要器官，只不过大多数情况两者是一致的。所以，我们常发现，AEBODSYS和AESOC 这两个变量虽然有着不同的含义，但两个变量值是相同的。

这张表SOC、PT对应所使用的变量是AEBODSYS、AEDECOD 。

具体编程过程

在输出TFL的过程中，我一般习惯分几个大步骤进行，以QC侧举例：

1. Create formats for output;
2. Get analysis data
3. Calculate statistics
4. Create dataset for QC
5. Compare

AE Table QC 流程图

下面具体描述一下这张表中的编程过程。

***1. Create formats for output;

为了程序功能显示简洁，我习惯把TFL编程中的Format设置，集中放到程序的开头。这样做，除了简洁之外，需要对Format进行更新时，也方便定位。

需要注意的是，对于包含统计量信息的Format，要在计算统计量之后在生成。相应的，位置也应该放到统计量计算之后。

频数的计数有很多种，我倾向于使用Means过程步，主要因为可以使用preloadfmt选项。SAS默认不输出分类计数为０的那一条记录，当遇到这种情况，我们需要补充计数为０的类别。使用preloadfmt选项，提前确定好分类变量的Format，在某一分类计数为０时，也可以输出该类的汇总结果，不必要做补全计数的操作。

演示代码如下：

***1. Create formats for output;
**Format for displaying all values of trt01an;
proc format;
  value trtn;
    1 = 1
    2 = 2
  ;
run;

***2. Get analysis data

生成AE的受试者发生率的表，涉及到２个数据集ADSL、ADAE。

获取数据集时，要注意数据集记录的筛选条件。一般在SAP中会指出，AE的分析属于安全性分析，这里的条件是saffl="Y"，通常指至少服用一次试验用药的数据集。

同时，安全性人群的分析，试验分组基于受试者实际用药情况，而不是计划用药的情况。具体到分组变量的选择是trtxxan(n), 而不是trtxxp(n)。

由于Means过程步只能对数值变量，进行分析，所以程序中还会新建一个变量(flag = 1)用于Means过程步的计数。

需要指出的是，这里获取分析数据，可以直接根据简单条件筛选出分析所需要的记录。如果筛选记录涉及统计量的计算，这一部分获取数据集的程序，我也会放到计算统计量程序之后。这一点跟前面设置Format类似。

ADAE选取时，一般需要对缺失值的SOC、PT进行处理。

演示代码如下：

***2. Get analysis data;
**2.1 Data for BigN;
data adsl;
  set adam.adsl;
  where saffl = "Y" and trt01an in (0,1);

   flag = 1; /*Flag for count*/
run;

**2.2 Data from adae;
data adae;
  set adae;
  where saffl = "Y" and trt01an in (0,1) and trtemfl = "Y";
  
  if aebodsys = "" then aebodsys = "_Missing System Organ Class";
  if aedecod = "" then aedecod = "_Missing Preferred Term";

  flag = 1; /*Flag for count*/
run;

这里获取的ADAE数据集，包含的“重复受试者”的信息，在计数AE、具体SOC、具体PT发生的人数时，还需要去重。这是这类表的一个重点，读者可以根据发生率的含义以及以下程序进行理解。

*Data for subjects count;
proc sort data = adae out = adae_subj nodupkey;
  by usubjid trt01an;
run;

*Data for SOC count;
proc sort data = adae out = adae_soc nodupkey;
  by usubjid trt01an aebodsys;
run;

*Data for PT count;
proc sort data = adae out = adae_pt nodupkey;
  by usubjid trt01an aebodsys aedecod;
run;

***3. Calculate statistics;

**3.1 Derive BigN and save them to macro vars

关于计算统计量，首先是BigN的计算。为方便方便后续调用，一般将BigN保存在宏变量中。常用的计数过程步为Means和Freq，我习惯使用Means过程步。

宏变量的生成，一般有两种方法。一是，SQL过程步中的INTO子句；另一个是，Data步中的call symput语句。赋值BigN时，我常使用后者，这里的原因，有空再另写文章说明，这里就不作详细解释。

宏变量命名时，可以将分组信息保留在名称中，便于调用时识别。通常，我也会在这一程序之前的Commets中，加入TRTN与TRT的对应关系的说明。

对于生成的宏变量以及取值，我会从SAS字典中读取保存到固定数据集中，方便编程过程中回看(*Check Bign)。

**3.1 Derive BigN and save them to macro vars;
*1: PLACEBO;
*2: TREATMENT;

proc means data = adsl nway completeTypes;
  format trt01an trtn.;
  class trt01an / preloadfmt order = data;
  
  var flag;
  output n = bign nmiss = nmiss out = Bing;
run;

data _null_;
  set Bign;
  call symput("N_"||strip(put(trt01an, best.)), strip(put(bign, best.)) );
run;

*Check Bign;
proc sql noprint;
  create table Bigncheck as 
    select *
    from dictionary.macros
    where name like "N_%";
quit;

**3.2 Count subject;

计数发生AE的人数时，需要使用按usubjid trt01an去重的数据集。选项completeTypes、preloadfmt以及format语句使得汇总结果包含所有trt01pn的选项（对于计算BigN，汇总结果中分组变量一般都种类齐全，我是出于习惯加上这些语句）。

输出结果后，需要将trt01an分组的纵向排列结果转置为横向，与TFLShell对应。

**3.2 Count subject;
proc means data = adae_subj  nway completeTypes format trt01pn trtn.;
  format trt01an trtn.;
  class trt01an / preloadfmt order = data;
  var flag;
  output n = subjn nmiss = nmiss out = subjn1;
run;

proc transpose data = subjn1 out = subjn prefix=_;
  id trt01an;
  var subjn;
run;

**3.3 Count SOC;

计数发生具体SOC的AE的人数，需要使用按usubjid trt01an aebodsys去重的数据集。汇总结果转置时，需要保留SOC的横向信息，即by aebodsys。

**3.3 Count SOC;
proc means data = adae_soc  nway completeTypes;
  class aebodsys;

  format trt01an trtn.;
  class trt01an / preloadfmt order = data;
  var flag;
  output n = socn nmiss = nmiss out = socn1;
run;

proc transpose data = socn1 out = socn2 prefix=_;
  by aebodsys;
  id trt01an;
  var subjn;
run;

由于表格需要按SOC频数降序、字母升序排列的顺序，在汇总结果输出后，需要获取SOC的排列顺序(by descending _2 aebodsys;)。

*Sort SOC by descending order of freq;
proc sort data = socn2 out = socn3;
  by descending _2 aebodsys;
run;

*Get the order var for SOC;
data socn;
  set socn3;
  soc_ord = _n_;

  proc sort;
    by aebodsys;
run;

**3.4 Count SOC＊PT;

计数发生具体PT的AE的人数，需要使用按usubjid trt01an aebodsys aedecod去重的数据集。汇总结果转置时，需要保留SOC和PT的横向信息，即by aebodsys aedecod。

计数PT的过程相比于SOC，是要复杂一点。需要保留PT对应的SOC信息，方便与SOC计数结果拼接时，能与对应的SOC计数排序划分在一组。

但是，这里计数时，不能使用preloadfmt、completeTypes选项。首先，数据集中SOC、PT出现的种类未知，无法提前定义Format，就不能使用preloadfmt选项。其次，一个PT对应一个SOC，completeTypes选项会补全PT与数据集中其他SOC相对应的分组情况，与我们想要的结果不符。

基于以上考虑，需要补全所有aebodsys aedecod trt01an的分组信息，我采用的是Dummy数据集的方法。

**3.4 Count SOC*PT;
proc means data = adae_pt nway;
  class aebodsys aedecod trt01an;
  var flag;
  output  n = ptn nmiss = nmiss out = ptn1;
run;

*Get all SOC and PT values for dummy dataset;
proc sort data = ptn1 out = socpt (keep = aebodsys aedecod) nodupkey;
  by aebodsys aedecod;
run;

*Create a dummy dataset for trt01an;
data dummy;
  set socpt;
  
  trt01an = 1; output;
  trt01an = 2; output;
run;

*Get SOC*PT count;
data ptn2;
  merge ptn1 dummy;
  by aebodsys aedecod trt01an;
  
  if ptn = . then ptn = 0;

  keep aebodsys aedecod trt01an ptn;
run;

proc transpose data = ptn2 out = ptn3 prefix=_;
  by aebodsys aedecod;
  id trt01an;
  var ptn;
run;

与计数SOC一样，获取汇总计数结果后，需要再获取PT的频数降序、字母排序升序的信息(by aebodsys descending _2 aedecod;)。PT的排序涉及两部分，PT所属SOC的排序以及PT在SOC内的排序。

SOC的排序在上一步（3.3 Count SOC）产生，后续可以通过aebodsys作为关键变量拼接获取，所以这一步只需要获取每个SOC内PT的排列顺序。

*Sort PT by descending order of freq in each SOC;
proc sort data = ptn3 out = ptn4;
  by aebodsys descending _2 aedecod;
run;

*Get the order var for PT;
data ptn;
  set ptn4;
  pt_ord = _n_;

  proc sort;
    by aebodsys;
run;

***4. Create dataset for QC;

在上一步中，表格所需要的统计量全都计算完毕，在进行比较之前，需要根据TFLShell以及读取RTF的结果，对统计量信息进行整合。

**4.1 Create dataset for header;

第一步，是进行header的信息的处理。我在Review其他人代码时发现，不少人做Table的QC并没有在程序中比对Header的信息，这是不规范的。

Header中的对应的信息有两类，一是具体的Label显示内容，二是BigN的大小。如果单单只靠人工去判断这两类内容，是很容易出差错的。

参考的代码如下：

**4.1 Create dataset for header;
data header;
  row_num = 0;

  length c1 - c3 $200;
  c1 = "System Organ Class Preferred Terrm";
  c2 = "Placebo (N = %sysfunc(strip(&N_1.))) n (%)";
  c3 = "Treatment (N = %sysfunc(strip(&N_2.))) n (%)";
run;

**4.2 Combine counting datasets above;

计数完成后，需要除以BigN计算发生率。考虑到，每个计数输出数据集的结构相同，统一计算发生率比较高效简洁，所以先将这些数据集纵向拼接在一起后，进行计算。

同时，也需要划分一些分组变量，为后续的排序做准备。

*Combine ;
data final1;
  set subjn(in = a) socn(in = b drop = soc_ord) ptn(in = c);

  length c1 - c3 $200;
  
  if a then do;
    sec = 1;
    c1 = "Number of subjects reporting treatment-emergent adverse events";
  end;

  if b then do;
    sec = 2;
    c1 = aebodsys;
  end;

  if c then do;
    sec = 2;
    c1 = aedecod;
  end;

  c2 = strip(put(_1, 8.)) || " (" || strip(put(_1/&N_1.*100,8.1)) || ") ";
  c3 = strip(put(_2, 8.)) || " (" || strip(put(_2/&N_2.*100,8.1)) || ") ";

  proc sort;
    by aebodsys;
run;

整张表的内容排序是先按SOC频数降序排序、再按照PT频数降序排序。所以需要将之前在计数时生成的排序变量拼接到输出数据集中。

一个SOC对应一个排序，一个SOC分类下可能有多个PT，我们需要为PT拼接对应的SOC顺序(soc_ord)。为避免变量覆盖，在上一步的程序中，已经把数据集中的soc_ord变量删除。PT也有对应的顺序(pt_ord)，排序时要加上这两个变量。

*Get SOC order;
data final2;
  merge final1(in = a) socn(in = b keep = aebodsys soc_ord);
  by aebodsys;
  
  if a and b or sec = 1;
  
  proc sort;
    by sec soc_ord pt_ord;
run;

我们公司宏在读取RTF时，默认是会添加读取内容的行号的。Header的行号为０，表格中的行号从１开始计数，QC数据集中需要增加行号的信息。

*Get row number;
data final;
  set final2;

  row_num = _n_;

  keep row_num c1-c3;
run;

最后，将Header信息与表格的主体信息拼接。

*Combine header and results;
data qc;
  set header final;
run;

***5. Compare;

各家公司可能都有自己一套QC的宏程序，我觉得自家公司宏程序输出的QC结果数据集有些臃肿，在完全Pass QC之前，我习惯使用简单的Compare过程步进行比对调试。

***5. Compare;
proc compare base = readrtf comp = qc out=df outbase outcomp outnoequal outdif;
run;

结语

以上是我做AE受试者发生率表的全部过程，读者可以对照流程图细细捋一遍，加深印象。

编程前，最好理解这张表的统计量的具体含义；编程中，要留意受试者去重以及SOC、PT的排序。

如果读者对此有疑惑，可以多看几遍相关文字描述和代码示例；也可以copy对应代码，到自己的项目环境中进行实际运行和调试。

以上内容也展示了我出TFL的思维框架和编程习惯，希望对读者有帮助。

感谢阅读， 欢迎关注：SAS茶谈！
若有疑问，欢迎评论交流！