1 引言

2016年底，針刺臨床研究論文《針刺治療慢性嚴(yán)重功能性便秘的隨機(jī)對(duì)照試驗(yàn)》見(jiàn)刊于內(nèi)科學(xué)年鑒，為針刺治療慢性便秘的療效提供了國(guó)際認(rèn)可的高質(zhì)量循證醫(yī)學(xué)證據(jù)。該研究共納入1075例患者，歷時(shí)3年完成，共有國(guó)內(nèi)15家三甲醫(yī)院參與。這項(xiàng)隨機(jī)對(duì)照研究通過(guò)中心分層區(qū)組隨機(jī)化的方法，設(shè)立電針組與假電針組兩組進(jìn)行平行對(duì)照，分別接受8周電針針刺穴位天樞、腹結(jié)、上巨虛與假電針淺刺雙側(cè)天樞旁、腹結(jié)旁、上巨虛旁非穴位點(diǎn)。結(jié)果發(fā)現(xiàn)8周電針治療可以安全有效地增加慢性嚴(yán)重便秘患者的完全自主排便次數(shù)，治療結(jié)束后療效可持續(xù)12周。

隨訪期的設(shè)置對(duì)于評(píng)價(jià)針刺的遠(yuǎn)期療效具有重要的意義，該研究的隨訪期為12周，但與之相矛盾的是隨著隨訪期長(zhǎng)度的增加，參與試驗(yàn)受試者的依從性則會(huì)下降，即脫落率與隨訪時(shí)間的長(zhǎng)度呈正相關(guān)關(guān)系。所以，如何利用恰當(dāng)?shù)慕y(tǒng)計(jì)分析方法處理缺失數(shù)據(jù)對(duì)于一項(xiàng)隨機(jī)對(duì)照試驗(yàn)十分重要。本文主要介紹內(nèi)容為在這項(xiàng)隨機(jī)對(duì)照試驗(yàn)中，針對(duì)缺失值處理的具體方法及基于SAS軟件的操作步驟進(jìn)行說(shuō)明。本文首先介紹如何使用SAS軟件對(duì)試驗(yàn)基線數(shù)據(jù)進(jìn)行統(tǒng)計(jì)描述。

2 研究項(xiàng)目的基線描述

2.1 問(wèn)題與數(shù)據(jù)

本研究共納入1 075例患者，人口統(tǒng)計(jì)學(xué)資料及基線特征按組別采用描述性統(tǒng)計(jì)方法展示，有關(guān)數(shù)據(jù)資料如下圖6-1所示。

2.2 數(shù)據(jù)結(jié)構(gòu)的分析

對(duì)研究資料圖6-1而言，試驗(yàn)因素是針灸治療(電針組和假電針組)和醫(yī)院(15家醫(yī)院)，具體變量信息如表6-1所示。

圖6-1 患者數(shù)據(jù)整理

2.3 分析的目的與統(tǒng)計(jì)分析方法的選擇

分析目的及方法：數(shù)據(jù)從excel導(dǎo)入SAS數(shù)據(jù)庫(kù)，按組別計(jì)算年齡、體重指數(shù)等連續(xù)性變量的平均值、中位數(shù)、四分位數(shù)、標(biāo)準(zhǔn)差、最小值和最大值；性別、民族、等分類(lèi)變量提供各類(lèi)別的例數(shù)和百分比。

2.4 SAS程序及說(shuō)明

proc import datafile='D:\test.xls'

表6-1 變量信息定義表

out=test dbms=excel replace;

getnames=yes;

run;

程序說(shuō)明：語(yǔ)句“proc import”表示調(diào)用import過(guò)程，進(jìn)行外部文件導(dǎo)入SAS過(guò)程。該語(yǔ)句中的“datafile=”規(guī)定要讀入外部文件的地址及名稱(chēng)；“out=”規(guī)定要輸出的SAS數(shù)據(jù)集，并命名為test；“dbms”選項(xiàng)規(guī)定外部數(shù)據(jù)文件格式的標(biāo)識(shí)名，excel文件可以統(tǒng)一設(shè)置為dbms=excel；replace規(guī)定替換已存在的文件，便于我們重復(fù)導(dǎo)入?！癵etnames=Yes”規(guī)定了導(dǎo)入數(shù)據(jù)第一行為變量名。

proc means data=test N NMISS MEAN STD MIN MAX Median Q1 Q3 CLM;

var BMI Course;

class group;

run;

程序說(shuō)明：以計(jì)算連續(xù)性變量BMI為例，語(yǔ)句“ proc means”表示調(diào)用means過(guò)程，計(jì)算變量的基本描述統(tǒng)計(jì)量。該語(yǔ)句中的“class”定義觀測(cè)組并分別計(jì)算觀測(cè)組的描述統(tǒng)計(jì)量；“VAR”規(guī)定要求計(jì)算簡(jiǎn)單描述統(tǒng)計(jì)量的數(shù)值變量和次序；“N NMISS MEAN STD MIN MAX Median Q1 Q3”規(guī)定了描述統(tǒng)計(jì)量關(guān)鍵詞分別是樣本量、含缺失值樣本數(shù)量、均數(shù)、標(biāo)準(zhǔn)差、最小值、最大值、中位數(shù)、1/4分位數(shù)和3/4分位數(shù)。

proc freq data=test;

tables sex*group/nopercent nocol chisq;

run;

程序說(shuō)明：該程序以分類(lèi)變量性別舉例，關(guān)鍵選項(xiàng)table語(yǔ)句表示輸出4個(gè)表，同時(shí)不顯示行百分比(nocol)、以及百分?jǐn)?shù)(nopercent)；chisq目的是輸出集中常用的卡方檢驗(yàn)統(tǒng)計(jì)量，列聯(lián)表為4個(gè)表時(shí)，還可以輸出Fisher精確檢驗(yàn)結(jié)果。

2.5 主要分析結(jié)果及解釋

圖6-2為在基線描述過(guò)程中proc import、proc means以及proc freq語(yǔ)句在SAS軟件中的基本格式。

圖6-3為以連續(xù)性變量BMI以及age為例，計(jì)算變量的基本描述統(tǒng)計(jì)量后的

輸出結(jié)果。

圖6-4為以分類(lèi)變量sex為例，計(jì)算變量的基本描述統(tǒng)計(jì)量后的輸出結(jié)果。

圖6-2 基線描述SAS程序代碼

圖6-3 兩組基線連續(xù)性指標(biāo)分析結(jié)果

圖6-4 兩組基線分類(lèi)指標(biāo)分析結(jié)果

3 主要指標(biāo)的敏感性分析

3.1 問(wèn)題與數(shù)據(jù)

在實(shí)際臨床試驗(yàn)中，常常會(huì)出現(xiàn)缺失數(shù)據(jù)。數(shù)據(jù)缺失有各種形式：譬如在臨床試驗(yàn)中患者中途退出試驗(yàn)研究，或者在某些時(shí)間點(diǎn)沒(méi)有進(jìn)行檢查，或者患者不情愿回答某些項(xiàng)目。這些缺失的原因也許與患者的病情有關(guān)，也許無(wú)關(guān)。

本研究對(duì)主要指標(biāo)的缺失值采用多重填補(bǔ)方法(Multiple imputation，MI)，該方法通常假設(shè)數(shù)據(jù)缺失機(jī)制為隨機(jī)缺失(Missing at random，MAR)。隨機(jī)缺失是目前較常見(jiàn)處理缺失值的方法，指反應(yīng)變量的缺失只依賴(lài)于已觀測(cè)到的反應(yīng)變量値，而與未觀測(cè)到的反應(yīng)變量値無(wú)關(guān)。然而，缺失數(shù)據(jù)無(wú)法被觀測(cè)，導(dǎo)致隨機(jī)缺失假設(shè)無(wú)法被證實(shí)。對(duì)于采用假定隨機(jī)缺失的臨床試驗(yàn)，國(guó)家研究委員會(huì)建議對(duì)違背MAR的統(tǒng)計(jì)假設(shè)進(jìn)行統(tǒng)計(jì)推斷，認(rèn)為主要指標(biāo)缺失機(jī)制的敏感性分析是統(tǒng)計(jì)報(bào)告中的必要組成部分。

本研究采用兩組等比例多中心隨機(jī)區(qū)組設(shè)計(jì)，治療組為接受電針試驗(yàn)組，對(duì)照組為接受假電針安慰劑組。變量Group為分組變量，Y0為主要指標(biāo)CSBM基線得分，Y1~Y4為治療期和隨訪期CSBM不同時(shí)間點(diǎn)測(cè)量值，有關(guān)資料如表6-2所示。

3.2 數(shù)據(jù)結(jié)構(gòu)分析

對(duì)研究資料表6-2而言，試驗(yàn)因素是針灸治療(電針組和假電針組)，試驗(yàn)效應(yīng)是治療效果(Y0-Y4)，具體變量信息如表6-3所示。

3.3 分析目的與統(tǒng)計(jì)分析方法選擇

據(jù)本研究統(tǒng)計(jì)分析計(jì)劃，假定缺失數(shù)據(jù)機(jī)制為隨機(jī)缺失，對(duì)研究主要指標(biāo)的缺失值處理采用了多重填補(bǔ)方法。本章節(jié)分析目的是，抽取378例患者數(shù)據(jù)，針對(duì)多重填補(bǔ)方式的隨機(jī)缺失假定進(jìn)行敏感性分析，采用模式混合模型(Pattern mixture model，PMM)(Little 1993；Molenberghs and Kenward 2007，pp.30，34–37)，該方法不僅適用于多種缺失機(jī)制的數(shù)據(jù)，且可在機(jī)制下得到更為準(zhǔn)確的參數(shù)估計(jì)值；并不要求一定知道缺失機(jī)制的具體分布形式，模型參數(shù)估計(jì)穩(wěn)健，該方法被國(guó)際醫(yī)學(xué)雜志如JAMA，內(nèi)科學(xué)年鑒等雜志的統(tǒng)計(jì)委員一致推薦(http:///aim/pages/AuthorInformationStatisticsOnly)。然而，該方法需要估計(jì)較多的特定模式參數(shù)，故一般形式下的PMM通常可識(shí)別性較低。為了增強(qiáng)模型的識(shí)別能力，需設(shè)定一些約束條件。根據(jù)本研究中觀察到的患者CSBM最大為10次，模型約束條件限定值選取10，該模型的基本原理及參數(shù)估計(jì)方法詳見(jiàn)。

表6-2 主要指標(biāo)數(shù)據(jù)結(jié)構(gòu)

表6-3 變量信息定義表

3.4 SAS程序及說(shuō)明

proc mi data=Mono3 seed=1423741;

class group;

var y0 y1 y2 y3 y4;

run;

程序說(shuō)明：首先檢查主要指標(biāo)缺失類(lèi)型。語(yǔ)句“proc mi”表示調(diào)用MI過(guò)程，該語(yǔ)句中的“data=test”是指分析的數(shù)據(jù)集名稱(chēng)為“test”；“seed=”是指種子數(shù)，可以隨機(jī)輸入；“Class”選項(xiàng)要求輸入分組變量；Var選項(xiàng)要求Y0~Y4進(jìn)行缺失值分析。

%midata( data=Mono2, smin=-2, smax=0, sinc=0.1, out=out1);

程序說(shuō)明：%midata是sas批處理程序，詳見(jiàn)附錄1。主要功能是對(duì)PMM模型進(jìn)行條件約束，并輸入結(jié)果。Smin選項(xiàng)要求約束條件最小值；smax選項(xiàng)要求約束條件最大值；sinc表示最小值到最大值的步長(zhǎng)；out選項(xiàng)輸入最終數(shù)據(jù)集為out1。

proc reg data=out1;

model y1= Trt y0;

by Shift _Imputation_;

ods output parameterestimates=regparms;

run;

程序說(shuō)明：語(yǔ)句“proc reg”表示調(diào)用reg過(guò)程，該語(yǔ)句中的“data=test”是指分析的數(shù)據(jù)集名稱(chēng)為“test”；Model選項(xiàng)中等號(hào)前放入因變量，等號(hào)后放入可能影響療效指標(biāo)的因素；by選項(xiàng)要求數(shù)據(jù)分析按照模型不同數(shù)據(jù)集及不同約束條件進(jìn)行分析，其中“shift”為PMM模型約束條件；“_Imputation_”為多重填補(bǔ)數(shù)據(jù)集標(biāo)識(shí)變量；“ods output parameterestimates=regparms”選項(xiàng)要求回歸模型結(jié)果輸出到同一數(shù)據(jù)集中，并命名該數(shù)據(jù)集為“regparms”。

proc mianalyze parms=regparms;

modeleffects Trt;

by Shift;

ods output parameterestimates=miparm1;

run;

程序說(shuō)明：語(yǔ)句“proc mianalyze”表示調(diào)用mianalyze過(guò)程，該語(yǔ)句中的parms選項(xiàng)制定數(shù)據(jù)來(lái)源集為regparms；modeleffects選項(xiàng)要求整合研究中組間結(jié)果，by選項(xiàng)要求數(shù)據(jù)分析按照不同約束條件進(jìn)行分析；“ods output parameterestimates=regparms”選項(xiàng)要求回歸模型結(jié)果輸出到同一數(shù)據(jù)集中，并命名該數(shù)據(jù)集為“miparm1”。

proc print label data=miparm1;

var Shift Probt;

format Probt 8.4;

run;

程序說(shuō)明：語(yǔ)句“ proc print”表示調(diào)用print過(guò)程，該語(yǔ)句中的“data=”要求調(diào)用數(shù)據(jù)集miparm1的數(shù)據(jù)，進(jìn)行后續(xù)分析；“var”選項(xiàng)要求對(duì)shift和Prrobt兩個(gè)變量進(jìn)行分析；“format”規(guī)定變量輸出格式，此處要求變量Probt小數(shù)點(diǎn)后4位。

3.5 主要分析結(jié)果及解釋

圖6-5顯示數(shù)據(jù)中數(shù)據(jù)缺失模式，各個(gè)變量意義：“Group”顯示原始數(shù)據(jù)的缺失類(lèi)型包括三種形式；“Freq”顯示每個(gè)類(lèi)型的的頻次；“Prercent”規(guī)定了每個(gè)類(lèi)型占樣本總數(shù)的百分比；COL1、COL2、COL3分別指示不同缺失類(lèi)型下，即將填補(bǔ)的三個(gè)變量平均數(shù)。

圖6-6顯示為在結(jié)果模式數(shù)據(jù)中不同shift下，每次填補(bǔ)中每組治療效應(yīng)(“估計(jì)”變量)以及可信區(qū)間(“下限”、“上限”變量)。

上表結(jié)果顯示在規(guī)定的shift(?10，10)的范圍內(nèi)，兩組間的差值及可信區(qū)間，結(jié)果顯示表明在規(guī)定的shift內(nèi)，兩組間依然存在統(tǒng)計(jì)學(xué)差異，P值均<0.05，結(jié)果穩(wěn)定。

圖6-5 數(shù)據(jù)缺失類(lèi)型結(jié)果

圖6-6 不同模式下兩組多重填補(bǔ)結(jié)果

表6-4 不同shift情況下兩組差值及可信區(qū)間

附錄1

PMM SAS macro code

/*----------------------------------------------------------------*/

/*--- Generate imputed data set for specified shift parameters ---*/

/*--- data= input data set ---*/

/*--- smin= min shift parameter ---*/

/*--- smax= max shift parameter ---*/

/*--- sinc= increment of the shift parameter ---*/

/*--- out= output imputed data set ---*/

/*----------------------------------------------------------------*/

%macro midata( data=, smin=, smax=, sinc=, out=);

data &out;

set _null_;

run;

/*------------ # of shift values ------------*/

%let ncase= %sysevalf( (&smax-&smin)/&sinc, ceil );

/*------- Imputed data for each shift -------*/

%do jc=0 %to &ncase;

%let sj= %sysevalf( &smin + &jc * &sinc);

proc mi data=&data seed=14823 nimpute=10 out=outmi;

class Trt;

monotone reg;

mnar adjust( y1 / shift=&sj adjustobs=(Trt='1') );

var Trt y0 y1;

run;

data outmi;

set outmi;

Shift= &sj;

run;

data &out;

set &out outmi;

run;

%end;

%mend midata;