如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料XV.用 SAS 軟件實(shí)現(xiàn)具有一個(gè)重復(fù)測(cè)量的兩因素和具有兩個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)定量資料的統(tǒng)計(jì)分析

王琪，胡良平，高輝

我們?cè)敿?xì)介紹了如何用 SAS 軟件實(shí)現(xiàn)具有一個(gè)重復(fù)測(cè)量的單因素設(shè)計(jì)定量資料的統(tǒng)計(jì)分析。本期，我們繼續(xù)探討生物醫(yī)藥研究中最常見(jiàn)的其他兩種重復(fù)測(cè)量設(shè)計(jì)，即具有一個(gè)重復(fù)測(cè)量的兩因素和具有兩個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)定量資料的統(tǒng)計(jì)分析的SAS 實(shí)現(xiàn)。

1 具有一個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)定量資料的SAS 實(shí)現(xiàn)

例1 選了 10只家兔觀察某藥物 A 對(duì)皮膚損傷情況。A有 2個(gè)水平：A1（高濃度），A2（低濃度）；為了排除實(shí)驗(yàn)部位對(duì)觀測(cè)結(jié)果的影響，藥物涂在兔的4個(gè)對(duì)稱(chēng)部位上，即部位因素 B 有 4個(gè)水平：B1（左前腿）、B2（右前腿）、B3（左后腿）、B4（右后腿）。10只家兔被完全隨機(jī)地均分到 A1、A2兩組中去，將藥物涂在任何一只兔 4條腿的對(duì)稱(chēng)部位上、且涂的面積相等，資料見(jiàn)表1，試進(jìn)行合適的統(tǒng)計(jì)分析。

表1 家兔皮膚損傷直徑

分析與SAS 實(shí)現(xiàn)：本例涉及“藥物濃度”和“部位”兩個(gè)實(shí)驗(yàn)因素，由于是在用藥后從同一只兔的不同部位重復(fù)獲得指標(biāo)“皮膚損傷直徑（mm）”的觀測(cè)值，所以，“部位”是一個(gè)重復(fù)測(cè)量因素，資料所對(duì)應(yīng)的實(shí)驗(yàn)設(shè)計(jì)類(lèi)型為具有一個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)。

具體的SAS 程序如下：

data sastjfx1; /*1*/do A=1 to 2;do rabbit=1 to 5;do B=1 to 4;input y @@; output;end; end; end;cards;16.520.718.319.516.719.017.918.416.518.516.417.618.022.519.220.217.520.018.719.918.221.219.418.918.520.519.618.3 proc mixed data=sastjfx1; /*5*/class A rabbit B;model y=A|B;repeated/type=AR(1) sub=rabbit(A);ods output fitstatistics=d;ods output dimensions=d1;run;%MACRO SHUJU(dataset,y); /*6*/data &dataset;set &dataset;&y=value;drop value;run;%MEND SHUJU;

19.722.521.320.521.523.722.421.920.722.721.621.8;run;proc mixed data=sastjfx1; /*2*/class A rabbit B;model y=A|B;repeated/type=VC sub=rabbit(A);ods output fitstatistics=a;ods output dimensions=a1;run;proc mixed data=sastjfx1; /*3*/class A rabbit B;model y=A|B;repeated/type=CS sub=rabbit(A);ods output fitstatistics=b;ods output dimensions=b1;run;proc mixed data=sastjfx1; /*4*/class A rabbit B;model y=A|B;repeated/type=UN sub=rabbit(A);ods output fitstatistics=c;ods output dimensions=c1;run;%SHUJU(a,VC) %SHUJU(b,CS)%SHUJU(c,UN) %SHUJU(d,AR1)%SHUJU(a1,VC) %SHUJU(b1,CS)%SHUJU(c1,UN) %SHUJU(d1,AR1)data e; /*7*/merge a b c d;run;data e1; /*8*/merge a1 b1 c1 d1;run;ods html;proc print data=e; /*9*/format _numeric_ 5.1;run;proc print data=e1; /*10*/run;ods html close;

程序說(shuō)明：SAS 程序中第 1 步為建立數(shù)據(jù)集；A 代表“藥物濃度”；B 代表“部位”；rabbit 代表“受試動(dòng)物個(gè)體”；y 代表觀測(cè)指標(biāo)“直徑”。第 2、3、4、5 步分別調(diào)用 mixed過(guò)程，采用 VC、CS、UN、AR(1) 四種協(xié)方差結(jié)構(gòu)模型對(duì)資料進(jìn)行方差分析。第 6 步為建立宏 SHUJU，以實(shí)現(xiàn)對(duì)數(shù)據(jù)集中已有變量 value的更名。第 7、8 步均用來(lái)實(shí)現(xiàn)對(duì)不同數(shù)據(jù)集的橫向合并。第 9、10 步均用來(lái)將數(shù)據(jù)集中的內(nèi)容輸出到 output 窗口中去。第 2、3、4、5 步所用語(yǔ)句基本相同，僅在“type=”后的選項(xiàng)不同，4個(gè)過(guò)程分別指定了 4 種協(xié)方差結(jié)構(gòu)模型。Mixed 過(guò)程中 repeated 語(yǔ)句用來(lái)規(guī)定個(gè)體的重復(fù)測(cè)量的協(xié)方差結(jié)構(gòu)，“/”后的sub（也可寫(xiě)為subject）用來(lái)指定數(shù)據(jù)集中的個(gè)體，若不含有分組因素，直接在“sub=”后面給出受試對(duì)象個(gè)體變量名稱(chēng)即可；若含有分組因素，則在“sub=”后面給出受試對(duì)象個(gè)體變量名稱(chēng)的同時(shí)，還需在后面加注“（）”，括號(hào)內(nèi)填入分組變量名稱(chēng)。在調(diào)用 mixed 過(guò)程進(jìn)行方差分析時(shí)，使用了兩個(gè) ods（output delivery system）語(yǔ)句，分別用來(lái)將模型擬合的有關(guān)信息（fitstatistics）和模型維度有關(guān)參數(shù)（dimensions）輸出。

SAS 輸出結(jié)果與結(jié)果解釋?zhuān)?/p>

Obs Descr VC CS UN AR11 –2 Res Log Likelihood 119.9 81.2 73.0 87.82 AIC (smaller is better) 121.9 85.2 93.0 91.83 AICC (smaller is better) 122.0 85.6 103.5 92.24 BIC (smaller is better) 122.2 85.8 96.0 92.4

Obs Descr VC CS UN AR11 Covariance Parameters 1 2 10 22 Columns in X 15 15 15 153 Columns in Z 0 0 0 04 Subjects 10 10 10 105 Max Obs Per Subject 4 4 4 4

這是上述程序中 ods 輸出的結(jié)果。首先給出了 4 種協(xié)方差結(jié)構(gòu)模型擬合本資料的有關(guān)情況，然后給出了協(xié)方差陣的有關(guān)信息（Covariance Parameters 表示模型中待估計(jì)的協(xié)方差結(jié)構(gòu)中參數(shù)的個(gè)數(shù)）。比較 4 種模型擬合資料情況的AIC、BIC 數(shù)值，可發(fā)現(xiàn) CS 與UN 兩種協(xié)方差結(jié)構(gòu)模型擬合資料較好。現(xiàn)比較這兩種協(xié)方差結(jié)構(gòu)模型擬合資料的效果之間的差異是否有統(tǒng)計(jì)學(xué)意義。

由ods 輸出結(jié)果的第二部分可知，q=2，q+v=10，所以v=8。而=15.51＞8.2，故P＞0.05。認(rèn)為可以選擇參數(shù)個(gè)數(shù)較少的協(xié)方差結(jié)構(gòu)模型，所以最后的結(jié)論應(yīng)按CS協(xié)方差結(jié)構(gòu)模型計(jì)算出來(lái)的結(jié)果來(lái)下。其假設(shè)檢驗(yàn)結(jié)果為：

Type 3 Tests of fixed effects

由上述結(jié)果可知：藥物濃度（A）、部位（B）、藥物濃度與部位的交互作用（A*B）均有統(tǒng)計(jì)學(xué)意義。

欲得出各均數(shù)之間兩兩比較的結(jié)果，可將原程序中所有過(guò)程步刪除，換上以下過(guò)程步：

程序運(yùn)行后，得如下結(jié)果：

Least squares means

Differences of least squares means

這是兩種藥物濃度的受試對(duì)象各部位皮膚損傷直徑均值之間兩兩比較的結(jié)果，首先給出了每種藥物濃度的受試對(duì)象各部位皮膚損傷直徑均值與0 比較的檢驗(yàn)結(jié)果，沒(méi)有任何實(shí)際意義。后面給出了兩兩比較的結(jié)果，A和B 列指定其中的某藥物濃度受試對(duì)象某個(gè)部位上的數(shù)據(jù)，_A和_B列也指定其中的某藥物濃度受試對(duì)象某個(gè)部位上的數(shù)據(jù)。

2 具有兩個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)定量資料的SAS 實(shí)現(xiàn)

例2 用5 條狗在背部做成同樣大小創(chuàng)傷 3個(gè)，分別用 3種藥物 A1、A2、A3治療，在用藥后第 2、4、6、8 天 4個(gè)時(shí)間點(diǎn)上（分別用 B1～B4表示）觀測(cè)創(chuàng)傷面積（cm2）。目的是觀察藥物對(duì)創(chuàng)傷面積影響的動(dòng)態(tài)變化情況，資料見(jiàn)表2，試進(jìn)行合適的統(tǒng)計(jì)分析。

表2 3 種藥物使用后不同時(shí)間點(diǎn)上觀測(cè)到的創(chuàng)傷面積

分析與SAS 實(shí)現(xiàn)：該資料涉及“藥物種類(lèi)”和“藥物作用時(shí)間”兩個(gè)實(shí)驗(yàn)因素，由于是在同一只狗身上重復(fù)使用3 種藥物（每種藥用在一個(gè)創(chuàng)傷部位），并在用藥后4個(gè)不同的時(shí)間點(diǎn)上重復(fù)測(cè)得每一處創(chuàng)傷的創(chuàng)傷面積，所以，“藥物種類(lèi)”和“藥物作用時(shí)間”都是重復(fù)測(cè)量因素，該資料所對(duì)應(yīng)的實(shí)驗(yàn)設(shè)計(jì)類(lèi)型為具有兩個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)。觀測(cè)指標(biāo)為“創(chuàng)傷面積（cm2）”，是定量指標(biāo)，故應(yīng)選用具有兩個(gè)重復(fù)測(cè)量的兩因素設(shè)計(jì)一元定量資料的方差分析來(lái)處理此資料。

具體的SAS 程序如下：

data sastjfx2; /*1*/do dog=1 to 5;do drug=1 to 3;do time=1 to 4;input y @@;output;end; end; end;cards;5.65.54.82.76.16.06.05.86.05.85.12.96.35.85.85.06.05.85.85.86.25.95.85.46.05.85.82.85.95.85.14.76.15.75.53.16.05.75.54.06.36.15.75.45.95.85.84.36.25.94.93.95.95.95.75.65.95.75.54.2;run;proc mixed data=sastjfx2; /*2*/class dog drug time;model y=drug|time;repeated/type=VC sub=dog;ods output fitstatistics=a;ods output dimensions=a1;run;proc mixed data=sastjfx2; /*3*/class dog drug time;model y=drug|time;repeated/type=CS sub=dog;ods output fitstatistics=b;ods output dimensions=b1;run;proc mixed data=sastjfx2; /*4*/class dog drug time;model y=drug|time;repeated/type=AR(1) sub=dog;ods output fitstatistics=c;ods output dimensions=c1;run;ods html;proc sql; /*5*/select a.Descr,a.VALUE as VC,b.VALUE as CS, c.VALUE as AR1 from a,b,c where a.Descr=b.Descr=c.Descr;quit;proc sql; /*6*/select a1.Descr,a1.VALUE as VC,b1.VALUE as CS, c1.VALUE as AR1 from a1,b1,c1 where a1.Descr=b1.Descr=c1.Descr;quit;ods html close;

程序說(shuō)明：SAS 程序中第 1 步為建立數(shù)據(jù)集；dog 代表“狗號(hào)”；drug 代表“藥物種類(lèi)”；time 代表“藥物作用時(shí)間”；y 代表觀測(cè)指標(biāo)“創(chuàng)傷面積”。第 2、3、4 步分別調(diào)用 mixed 過(guò)程，采用 VC、CS、AR(1) 三種協(xié)方差結(jié)構(gòu)模型對(duì)資料進(jìn)行方差分析。第 5、6 步將不同數(shù)據(jù)集的內(nèi)容進(jìn)行合并查詢(xún)，這與上一個(gè)程序中的宏 SHUJU 功能相同。由于采用 UN 協(xié)方差結(jié)構(gòu)模型對(duì)資料進(jìn)行方差分析時(shí)，中間計(jì)算無(wú)法形成正定矩陣，分析過(guò)程中斷。所以，上面沒(méi)有給出用 UN 協(xié)方差結(jié)構(gòu)模型對(duì)資料進(jìn)行方差分析的程序。而 SP（POW）要求重復(fù)測(cè)量因素均為定量變量，本資料中的藥物種類(lèi)是一個(gè)定性變量，所以也沒(méi)有采用 SP（POW）協(xié)方差結(jié)構(gòu)模型對(duì)資料進(jìn)行方差分析。

SAS 輸出結(jié)果與結(jié)果解釋?zhuān)?/p>

Description VC CS AR1–2 Res Log Likelihood 84.6 78.4 79.9 AIC (smaller is better) 86.6 82.4 83.9 AICC (smaller is better) 86.6 82.6 84.1 BIC (smaller is better) 86.2 81.6 83.1 Description VC CS AR1 Covariance Parameters 1 2 2 Columns in X 20 20 20 Columns in Z 0 0 0 Subjects 5 5 5 Max Obs Per Subject 12 12 12

這是上述程序中 ods 輸出的結(jié)果。首先給出了 3 種協(xié)方差結(jié)構(gòu)模型擬合本資料的有關(guān)情況，然后給出了協(xié)方差陣的有關(guān)信息。比較 3 種模型擬合資料情況的AIC、BIC 數(shù)值（越小越好），可發(fā)現(xiàn) CS 協(xié)方差結(jié)構(gòu)模型擬合資料較好，而 VC 協(xié)方差結(jié)構(gòu)模型擬合資料效果雖不如 CS 協(xié)方差結(jié)構(gòu)模型好，但其模型中參數(shù)個(gè)數(shù)要比 CS 協(xié)方差結(jié)構(gòu)模型少。現(xiàn)比較這兩種協(xié)方差陣模型擬合資料的效果之間的差異是否有統(tǒng)計(jì)學(xué)意義。

由 ods 輸出結(jié)果的第二部分可知，q=1，q+v=2，所以v=1。而=3.84＞6.2，故 P＞0.05。可認(rèn)為不適合用VC模型取代CS 模型，所以最后的結(jié)論應(yīng)按CS協(xié)方差結(jié)構(gòu)模型計(jì)算出來(lái)的結(jié)果來(lái)下。其假設(shè)檢驗(yàn)結(jié)果為：

Type 3 Tests of fixed effects

由上述結(jié)果可知：藥物種類(lèi)（drug）、藥物作用時(shí)間（time）以及兩者的交互作用（drug*time）均有統(tǒng)計(jì)學(xué)意義。即不同時(shí)間點(diǎn)上測(cè)得的家犬創(chuàng)傷面積均值之間的差異有統(tǒng)計(jì)學(xué)意義，不同藥物對(duì)創(chuàng)傷面積均值影響之間的差異有統(tǒng)計(jì)學(xué)意義。

欲得出各均數(shù)之間兩兩比較的結(jié)果，可將原程序中所有過(guò)程步刪除，換上以下過(guò)程步：

此步運(yùn)行結(jié)果與上述 type=CS 過(guò)程步計(jì)算結(jié)果基本相同，僅多了多個(gè)均數(shù)兩兩比較的結(jié)果，輸出結(jié)果所占篇幅較多，此處從略。

需要說(shuō)明的是：有些重復(fù)測(cè)量設(shè)計(jì)定量資料中含有協(xié)變量，如實(shí)驗(yàn)中對(duì)每一個(gè)個(gè)體在不同的時(shí)間點(diǎn)上進(jìn)行了重復(fù)測(cè)量，其中在“服藥前”觀測(cè)了 1 次，在“服藥后”觀測(cè)了 k（k ≥ 2）次，這屬于帶有協(xié)變量的重復(fù)測(cè)量設(shè)計(jì)問(wèn)題。對(duì)這類(lèi)資料，最合適的分析方法是以服藥前觀測(cè)的結(jié)果作為“基礎(chǔ)值”（即作為“協(xié)變量”），采用帶有協(xié)變量的重復(fù)測(cè)量設(shè)計(jì)資料定量資料的協(xié)方差分析。當(dāng)然，對(duì)于具有一個(gè)重復(fù)測(cè)量的單因素設(shè)計(jì)定量資料，由于再無(wú)其他實(shí)驗(yàn)分組因素，故即便含有“基礎(chǔ)值”，也不需要采用協(xié)方差分析。關(guān)于帶有協(xié)變量的重復(fù)測(cè)量設(shè)計(jì)定量資料的統(tǒng)計(jì)分析問(wèn)題，請(qǐng)讀者參考相關(guān)文獻(xiàn)。

[1]Hu LP.Application of statistical triple-type theory in the experiment design.Beijing: People’s Military Medical Press, 2006:107-120.(in Chinese)胡良平.統(tǒng)計(jì)學(xué)三型理論在實(shí)驗(yàn)設(shè)計(jì)中的應(yīng)用.北京:人民軍醫(yī)出版社, 2006:107-120.

[2]Hu LP.Scientific research design and statistical analysis of cardiovascular disease.Beijing: People’s Military Medical Press,2010:93-111.(in Chinese)胡良平.心血管病科研設(shè)計(jì)與統(tǒng)計(jì)分析.北京:人民軍醫(yī)出版社,2010:93-111.