如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料VIII. 用SAS 軟件實現(xiàn)拉丁方設(shè)計定量資料的統(tǒng)計分析

陶麗新，胡良平，郭晉

?

如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料VIII. 用SAS 軟件實現(xiàn)拉丁方設(shè)計定量資料的統(tǒng)計分析

陶麗新，胡良平，郭晉

100850 北京，軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心

編者按

生物統(tǒng)計學(xué)是生物學(xué)領(lǐng)域科學(xué)研究和實際工作中必不可少的工具，在分子生物學(xué)迅速發(fā)展的今天，生物統(tǒng)計學(xué)更顯示出了它的重要性。實驗設(shè)計與數(shù)據(jù)統(tǒng)計分析是現(xiàn)代生物學(xué)的基石，是生物學(xué)研究者檢驗假說、尋找模式、建立生物學(xué)理論的有利工具，也是生物學(xué)研究者探索微觀和宏觀生物世界的必備基礎(chǔ)知識。對于每天甚至是每時每刻涌現(xiàn)的大量的、以天文數(shù)字計量的分子遺傳數(shù)據(jù)，必須借助統(tǒng)計學(xué)知識加以分析處理，才能從中獲得有意義的信息?！吧锒鄻有詳?shù)據(jù)分析”是開展生物多樣性研究的一個重要方面，數(shù)據(jù)分析能力的高低極大地影響著我們對各種生態(tài)學(xué)現(xiàn)象認識的深度和廣度?，F(xiàn)在，電子計算機的普及使得生物統(tǒng)計分析過程大大簡化，生物統(tǒng)計分析軟件包的普及將生物統(tǒng)計學(xué)從統(tǒng)計學(xué)家的書本里解放了出來，簡化了生物統(tǒng)計分析過程，使之成為生物學(xué)研究者的常用工具。本刊特邀軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心主任胡良平教授，以“如何用 SAS 軟件正確分析生物醫(yī)學(xué)科研資料”為題，撰寫系列統(tǒng)計學(xué)講座，希望該系列講座能對生物醫(yī)學(xué)科研工作者有所幫助。

1 拉丁方設(shè)計的概念及應(yīng)用場合

拉丁方設(shè)計是當(dāng)實驗中只考察 1 個水平的實驗因素，但同時又涉及 2 個都具有水平的區(qū)組因素（即重要的非實驗因素），且它們之間的交互作用無統(tǒng)計學(xué)意義時所采用的實驗設(shè)計類型。

正交拉丁方設(shè)計是在當(dāng)實驗中只考察 1 個水平的實驗因素，但同時又涉及3 個都具有水平的區(qū)組因素，且它們之間的交互作用無統(tǒng)計學(xué)意義時所采用的實驗設(shè)計類型，又叫作希臘拉丁方設(shè)計。

很顯然，拉丁方設(shè)計和希臘拉丁方設(shè)計是根據(jù)實驗所涉及到的區(qū)組因素的個數(shù)來區(qū)分的，有時候也將這 2 種實驗設(shè)計類型分別稱為三因素拉丁方設(shè)計和四因素拉丁方設(shè)計。

拉丁方設(shè)計適用于下列場合：若每個區(qū)組中的個數(shù)據(jù)取自條件相同或相近的個受試對象，那么實驗因素對觀測指標(biāo)不會產(chǎn)生“攜帶效應(yīng)”的影響，往往是以“較大的樣本含量為代價的”。但是如果每個區(qū)組中的個數(shù)據(jù)重復(fù)測自相同的受試對象，雖然可以達到減少受試對象個數(shù)的目的，但是由于存在“攜帶效應(yīng)”（前一種處理對觀測指標(biāo)造成的影響在受試對象接受下一個處理時仍然存在），使其一般不適用于干預(yù)性實驗研究，除非能夠保證實驗因素各水平對觀測指標(biāo)的影響是短暫可逆的，即在設(shè)定的時間間隔內(nèi)，前一種處理對觀測指標(biāo)的影響會消除，觀測指標(biāo)的取值會恢復(fù)到最原始的水平的觀察性研究中應(yīng)用拉丁方設(shè)計效果較好；拉丁方設(shè)計和正交拉丁方設(shè)計都要求實驗所涉及的3 個或 4 個實驗因素之間不存在交互作用或者交互作用可以忽略不計且各因素的水平數(shù)相同（以實驗因素的水平數(shù)為基準(zhǔn)）[1-2]。這2 類拉丁方設(shè)計的構(gòu)造方法和列表格式參見下面的實例。

2 拉丁方設(shè)計一元定量資料方差分析的 SAS 實現(xiàn)

例 1 受試對象為 7 名志愿受試的健康女大學(xué)生，年齡在 19 ~ 21 歲，均有多次暈車史。正式試驗前進行的運動病敏感試驗證明她們確實為運動病敏感者，運動病耐力得分為（14 ± 7.4）分。藥物及分組：安慰劑為淀粉片（100 mg/片·人）（A）；暈海寧（50 mg/片·人）（B）；山莨菪堿（10 mg/片·人）（C）；腦益嗪組（25 mg/片·人）（D）；暈可平（30 ml/人）（E），主要成分為黛赭石、車前子等；生姜合劑（30 ml/人）（F），主要成分為生姜、明天麻等；復(fù)方中西藥制劑（G），每人使用上述生姜合劑 15 ml 加山莨菪堿 2.5 mg、腦益嗪 6.25 mg。藥物均口服。

實驗方法：①運動病誘發(fā)試驗：被試者坐在電動轉(zhuǎn)椅上，頭部固定，然后按 A 程式、B 程式順序進行反復(fù)急停刺激（采用“Modified Graybiel Method”[3]）。A 程式：被試者閉眼，作順時針旋轉(zhuǎn)，以 10o·s-2加速到 100o·s-2（10 s 內(nèi)），恒速旋轉(zhuǎn) 20 s，然后在 0.5 s 內(nèi)急停，休息30 s，再重復(fù)進行，共 20 次；B 程式：被試者開眼，作逆時針旋轉(zhuǎn)，參數(shù)同上，亦在 0.5 s 內(nèi)急停，休息 30 s 后再重復(fù)進行，共20 次。以被試者感到惡心為終止點。耐受上述刺激 A 程式每次計 1 分，B 程式每次計 2 分。②被試者于飯后 1 h 進入實驗室，按順序給藥。服藥 1 h 后進行運動病誘發(fā)試驗，計算被試者每次的運動病耐力得分。試驗于每周六或周日進行，連續(xù) 7 周完成。設(shè) 7 種藥物編號為A、B、C、D、E、F、G。7 名志愿者編號分別為1、2、3、4、5、6、7，如表 1 所示。

表1 不同藥物、不同志愿者、不同用藥順序?qū)\動耐力得分的影響

分析與 SAS 實現(xiàn)：該研究資料中，樣本含量為 7，試驗因素為“藥物種類”，它有 7個水平；同時還涉及 2 個區(qū)組因素分別為“受試者的個體差異”和“用藥次序”，分別有 7 個水平，并且因素之間的交互作用可以忽略不計，定量的觀測指標(biāo)為“運動耐力得分”，因此該資料應(yīng)該為 7 × 7 拉丁方設(shè)計一元定量資料。

對于該類資料進行統(tǒng)計分析，首先檢查資料是否滿足方差分析的前提條件，若滿足則采用 7 × 7 拉丁方設(shè)計一元定量資料的方差分析；若不滿足，則可以采用適當(dāng)?shù)淖兞孔儞Q方法使變換后的資料滿足方差分析的前提條件，然后再對變換后的資料采用 7 × 7 拉丁方設(shè)計一元定量資料的方差分析。需要注意的是，要處理好“攜帶效應(yīng)”的問題，因為試驗因素是藥物種類，因此應(yīng)該確保前一次用藥對后一次用藥的影響是微弱的，最好是在間隔一段時間后，定量指標(biāo)的取值能基本恢復(fù)到原先的水平，原文試驗每次用藥間隔 1 周時間。

具體的 SAS 程序（注：檢查資料是否具備參數(shù)檢驗前提條件的程序從略，下同）：

data a; do sub=1 to 7; do ord=1 to 7; input med $ value@@; output; end; end;cards;A 33 B 38 C 36 D 50 E 39 F 48 G 41B 35 E 39 A 32 G 36 F 42 D 26 C 38C 37 F 51 G 43 B 32 D 33 A 30 E 34D 40 G 38 E 38 F 45 C 34 B 35 A 33E 31 D 28 B 37 C 31 A 31 G 39 F 40F 44 C 32 D 43 A 32 G 44 E 32 B 33G 40 A 31 F 37 E 33 B 33 C 33 D 39;run;ods html;proc anova data=a; class sub ord med; model value=sub ord med; means med/snk;run;quit;ods html close;

程序說明：數(shù)據(jù)步建立數(shù)據(jù)集 a，輸入變量名和變量值，“sub”、“ord”、“med”和“value”分別代表志愿者編號、用藥順序、藥物種類和運動耐力得分。過程步調(diào)用適于分析均衡數(shù)據(jù)的方差分析過程，其過程名為 anova；class 語句指出在 anova 過程中要使用的分類變量；model 語句用來規(guī)定因變量和自變量效應(yīng)；means 語句用來計算藥物種類所對應(yīng)的運動耐力得分的均值，snk 選項表示對 7 種藥物使用后的運動耐力得分的平均值進行兩兩比較，quit 語句用來結(jié)束 anova 過程，因為 anova 過程為交互式的，run 語句不能結(jié)束該過程，它只是告訴 anova 去執(zhí)行 run 語句之前的語句。

SAS 輸出結(jié)果與結(jié)果解釋：

SourceDFAnova SSMean squareF valuePr > F sub6205.7142857 34.28571431.980.1002 ord6 41.7142857 6.95238100.400.8722 med6696.8571429116.14285716.700.0001

從以上輸出結(jié)果中可以看出，不同的受試者和不同的用藥順序分別對運動耐力得分均值影響之間的差異均無統(tǒng)計學(xué)意義，而不同的藥物種類對運動耐力得分均值影響之間的差異則有統(tǒng)計學(xué)意義（= 6.70，= 0.0001 < 0.05）。

此處藥物因素全部水平下均值之間的兩兩比較結(jié)果從略，依下面簡化后模型輸出的結(jié)果為準(zhǔn)。

從第一部分輸出結(jié)果中可以看出，不同受試對象和不同的用藥順序?qū)\動耐力得分均值的影響之間的差異無統(tǒng)計學(xué)意義，可以在 SAS 程序的 class 和 model 語句中將 sub 和 ord 刪除，這樣可以增大誤差項的自由度，使結(jié)果更加穩(wěn)定可靠。重新運行程序，得到檢驗結(jié)果如下：

SourceDFAnova SSMean squareF valuePr > F med6696.8571429116.14285716.36< 0.0001

Means with the same letter are not significantly different. SNK groupingMeanNmed A43.8577F BA40.1437G BC37.0007D BC35.1437E BC34.7147B BC34.4297C C31.7147A

從以上結(jié)果中可以看出，生姜合劑組與安慰劑組、暈海寧組、山莨菪堿組、腦益嗪組、暈可平組以及安慰劑組與復(fù)方中西藥制劑組的運動耐力得分均值之間的差異具有統(tǒng)計學(xué)意義，其他兩兩組之間的運動耐力得分均值之間的差異沒有統(tǒng)計學(xué)意義。生姜合劑組的運動耐力得分均值明顯大于安慰劑組、暈海寧組、山莨菪堿組、腦益嗪組、暈可平組的運動耐力得分均值；復(fù)方中西藥制劑組的運動耐力得分均值明顯大于安慰劑組的運動耐力得分均值。

例 2 在某項研究中，將4 種不同劑量的某藥液分別注射于 4 個受試對象，每個受試對象先后接受此藥物的不同的劑量各 1 次（即每個受試對象接受藥物 4 次），資料如表 2 所示。

表 2 不同的給藥次序、不同的受試對象、不同的藥物劑量對血糖升高值的影響

分析與 SAS 實現(xiàn)：該資料中樣本量為 4，試驗中涉及 1 個試驗因素，即“藥物劑量”，它有 4 個水平；同時還涉及 2 個區(qū)組因素：“給藥次序”和“受試者編號”，均有 4 個水平，且交互作用可以忽略不計，定量觀測指標(biāo)為“血糖的升高值（mg/dl）”，因此該資料為 4 × 4 拉丁方設(shè)計一元定量資料。

先檢查資料是否滿足參數(shù)檢驗的前提條件，若滿足，則采用 4 × 4 拉丁方設(shè)計一元定量資料的方差分析；若不滿足，則采用適當(dāng)?shù)淖兞孔儞Q方法使變換后的資料滿足方差分析的前提條件，然后再對變換后的資料采用 4 × 4 拉丁方設(shè)計一元定量資料的方差分析。

具體的 SAS 程序：

data b;do ord=1 to 4; do sub=1 to 4; input med $ value@@;output; end; end;cards;C 42 B 96 A 53 D 110B 50 A 31 D 78 C 55A 50 D 64 C 55 B 70D 98 C 41 B 79 A 49;ods html;proc glm;class ord sub med;model value=ord sub med;means med/snk;run;quit;ods html close;

程序說明：該例資料的 SAS 程序與例 1 的程序唯一不同之處是：過程步調(diào)用了 glm過程，glm 過程主要用于對非均衡數(shù)據(jù)進行方差分析，均衡數(shù)據(jù)也可以使用該過程，此時，可比 anova 過程消耗更多的計算機資源。對該資料而言，2 個過程的 SAS 輸出結(jié)果是基本相同的。

SAS 輸出結(jié)果與結(jié)果解釋：

SourceDFType III SSMean squareF valuePr > F ord31049.187500 349.7291671.530.2996 sub3 423.687500 141.2291670.620.6277 med34913.1875001637.7291677.180.0207

從上面結(jié)果可以看出，給藥次序和受試者的個體差異對血糖升高值的影響之間的差別沒有統(tǒng)計學(xué)意義，藥物劑量對血糖升高值的影響之間的差異具有統(tǒng)計學(xué)意義（= 7.18，= 0.0207 < 0.05）。

此處藥物劑量全部水平下均值之間的兩兩比較結(jié)果從略，依下面簡化后模型輸出的結(jié)果為準(zhǔn)。

由于給藥次序和受試者編號對血糖升高值之間的差別沒有統(tǒng)計學(xué)意義，可以將 SAS 程序 class 和 model 語句中的 ord 和 sub 刪除，這樣可以使誤差項的自由度增大，使結(jié)論更加可靠。重新編寫程序并發(fā)送給系統(tǒng)執(zhí)行，得到下面的檢驗結(jié)果：

SourceDFType III SSMean squareF valuePr > F med34913.1875001637.7291676.920.0059

Means with the same letter are not significantly different. SNK groupingMeanNmed A87.504D BA73.754B B48.254C B45.754A

從上面的結(jié)果中可以看出，該藥物 D 劑量與 C 劑量對血糖升高值的影響之間的差別具有統(tǒng)計學(xué)意義，D 劑量組的血糖升高值明顯高于 C 劑量組的血糖升高值；該藥物 D 劑量與 A 劑量對血糖升高值的影響之間的差別具有統(tǒng)計學(xué)意義，D 劑量組的血糖升高值明顯高于 A 劑量組的血糖升高值；其他任何 2 種劑量對血糖升高值的影響之間的差別沒有統(tǒng)計學(xué)意義。

3 正交拉丁方設(shè)計一元定量資料方差分析的 SAS 實現(xiàn)

例 3 在應(yīng)用重組人腫瘤壞死因子（rh-TNF）促進傷口愈合的實驗研究中，欲比較其 5 種不同劑量的效應(yīng)，同時要觀察動物個體間、不同用藥部位間、給藥順序間的差別。在該研究中，每隔 1 日更換內(nèi)層敷料，滴藥 1 次，記錄傷口愈合時間：

表3中，在行間隨機分配 5 只不同家兔，列間每只家兔背部切除位置不同的全層皮膚的 5 塊 2 cm × 2 cm 創(chuàng)面，希臘字母為每次用藥的 5 種先后順序，拉丁字母代表 5 種不同用藥滴數(shù)，即：A（0）、B（1）、C（2）、D（3）、E（4）。

分析與 SAS 實現(xiàn)：在該研究資料中，樣本量為 5，試驗因素為“重組人腫瘤壞死因子的劑量”，有 5 個水平；同時涉及到 3 個區(qū)組因素，分別為“家兔的個體差異”、“用藥順序”和“用藥部位”，分別有5 個水平，并且因素之間的交互作用可以忽略不計，定量的觀測指標(biāo)為“家兔創(chuàng)面的愈合時間（d）”，因此該資料為 5 × 5 的正交拉丁方設(shè)計一元定量資料。

當(dāng)資料滿足方差分析的前提條件時，使用 5 × 5 正交拉丁方設(shè)計一元定量資料的方差分析；若不滿足，則采用適當(dāng)?shù)淖兞孔儞Q方法使變換后的資料滿足方差分析的前提條件，然后再對變換后的資料采用 5 × 5 正交拉丁方設(shè)計一元定量資料的方差分析。

表 3 家兔創(chuàng)面不同用藥量的愈合時間（d）

具體的 SAS 程序：

data c; do sub=1 to 5; do buwei=1 to 5; input jiliang $ ord $ time@@;output; end; end;cards;A α 26 B β 21 C γ 20 D δ 17 E ε 18B δ 20 C ε 19 D α 18 E β 17 A γ 25C β 19 D γ 20 E δ 18 A ε 24 B α 20D ε 18 E α 18 A β 24 B γ 21 C δ 19E γ 17 A δ 24 B ε 21 C α 18 D β 17;run;ods html;proc anova; class sub buwei ord jiliang; model time=sub buwei ord jiliang; means jiliang/snk;run;quit;ods html close;

SAS 輸出結(jié)果與結(jié)果解釋：

SourceDFAnova SSMean squareF valuePr > F sub4 2.9600000 0.74000001.370.3257 buwei4 2.9600000 0.74000001.370.3257 ord4 3.3600000 0.84000001.560.2753 jiliang4161.360000040.340000074.70< 0.0001

從以上結(jié)果可以看出，家兔的個體差異、不同的用藥部位、給藥順序?qū)谟蠒r間造成的影響之間的差異沒有統(tǒng)計學(xué)意義，而藥物的不同劑量對傷口愈合時間造成的影響之間的差異具有統(tǒng)計學(xué)意義（= 74.70，< 0.0001）。

此處用藥滴數(shù)全部水平下均值之間的兩兩比較結(jié)果從略，依下面簡化后模型輸出的結(jié)果為準(zhǔn)。

由于家兔的個體差異、不同的用藥部位、給藥順序?qū)彝脛?chuàng)面?zhèn)谟蠒r間造成的影響之間的差異沒有統(tǒng)計學(xué)意義，可以將 SAS 程序 class 和 model 語句中的 sub、buwei、order刪除，重新運行程序，得到檢驗結(jié)果如下：

SourceDFAnova SSMean squareF valuePr > F jiliang4161.360000040.340000059.32< 0.0001

Means with the same letter are not significantly different. SNK groupingMeanNjiliang A24.60005A B20.60005B C19.00005C DC18.00005D D17.60005E

從上述結(jié)果中可以看出，用藥滴數(shù)為“2（C）”與“3（D）”之間、“3（D）”與“4（E）”之間在家兔創(chuàng)面愈合時間上的差異沒有統(tǒng)計學(xué)意義，其他兩兩之間比較的差異都具有統(tǒng)計學(xué)意義。用藥滴數(shù)為“0（A）”組的創(chuàng)面愈合天數(shù)明顯長于其他各組；用藥滴數(shù)為“1（B）”組也明顯長于“2（C）”以上各組；用藥滴數(shù)為“2（C）”組也長于“4（E）”組。從該實驗可以看出用藥劑量為隔日 3 滴的重組人腫瘤壞死因子有助于加快創(chuàng)面愈合。

[1] Hu LP. Application of statistics triple-type theory in the experiment design. Beijing: People’s Military Medical Press, 2006:71-76. (in Chinese)

胡良平. 統(tǒng)計學(xué)三型理論在實驗設(shè)計中的應(yīng)用. 北京: 人民軍醫(yī)出版社, 2006:71-76.

[2] Hu LP. The practical course in the statistical analysis for windows SAS, version 6.12 & 8.0. Beijing: Press of Military Medical Sciences, 2001:203-206. (in Chinese)

胡良平. Windows SAS 6.12 & 8.0實用統(tǒng)計分析教程. 北京: 軍事醫(yī)學(xué)科學(xué)出版社, 2001:203-206.

[3] Graybiel A, Lackner JR. A sudden-stop vestibulovisual test for rapid assessment of motion sickness manifestations. Aviat Space Environ Med, 1980, 51(1):21-23.

胡良平，Email：lphu812@sina.com

10.3969/cmba.j.issn.1673-713X.2010.02.016