如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料 X. 用SAS軟件實(shí)現(xiàn)嵌套設(shè)計(jì)和裂區(qū)設(shè)計(jì)定量資料的統(tǒng)計(jì)分析

陶麗新，胡良平

如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料 X. 用SAS軟件實(shí)現(xiàn)嵌套設(shè)計(jì)和裂區(qū)設(shè)計(jì)定量資料的統(tǒng)計(jì)分析

陶麗新，胡良平

編者按

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

1 嵌套設(shè)計(jì)及其資料的統(tǒng)計(jì)分析

1.1 嵌套設(shè)計(jì)的概念及應(yīng)用場合

試驗(yàn)中涉及兩個(gè)或多個(gè)試驗(yàn)因素，且依據(jù)專業(yè)知識可以認(rèn)為各試驗(yàn)因素對觀測指標(biāo)的影響有主次之分，主要因素各水平下嵌套著次要因素，次要因素各水平下又嵌套著更次要的因素，這樣的試驗(yàn)設(shè)計(jì)稱為嵌套設(shè)計(jì)。此類設(shè)計(jì)有兩種情形：第一種情形是，受試對象本身具有分組再分組的各種分組因素，處理（即最終的試驗(yàn)條件）是各因素各水平的全面組合，且因素之間在專業(yè)上有主次之分（如年齡與性別對心室射血時(shí)間的影響，性別的影響大于年齡）；第二種情形是，受試對象本身并非具有分組再分組的各種分組因素，處理（即最終的試驗(yàn)條件）不是各因素各水平的全面組合，而是各因素按其隸屬關(guān)系系統(tǒng)分組，且因素之間在專業(yè)上有主次之分（如研究不同代次不同家庭成年男性的身高資料，不同家庭之間的差別大于同一個(gè)家庭內(nèi)部不同代次之間的差別）。

嵌套設(shè)計(jì)主要應(yīng)用于試驗(yàn)因素對觀測指標(biāo)的影響有主次之分的試驗(yàn)研究中，試驗(yàn)因素之間的主次關(guān)系要有專業(yè)依據(jù)，不能憑空想象而定[1-2]。嵌套設(shè)計(jì)的構(gòu)造方法和列表格式參見例 1。

1.2 嵌套設(shè)計(jì)定量資料方差分析的 SAS 實(shí)現(xiàn)

例 1在4塊相同的試驗(yàn)田中考察肥力和小麥品種對小麥產(chǎn)量的影響，并且假定品種和肥力對產(chǎn)量的影響大小順序?yàn)槠贩N > 肥力，資料如表1所示。試分析小麥品種和肥力對小麥產(chǎn)量的影響有無統(tǒng)計(jì)學(xué)意義。

分析與 SAS 實(shí)現(xiàn)：該資料涉及兩個(gè)試驗(yàn)因素，分別為品種與肥力，各有3個(gè)水平，兩因素各水平全面組合，各種組合條件下進(jìn)行了 4 次獨(dú)立重復(fù)試驗(yàn)；又由于品種與肥力對產(chǎn)量的影響大小順序?yàn)槠贩N > 肥力，所以該資料為系統(tǒng)分組（或嵌套）設(shè)計(jì)。

表1 不同小麥品種和肥力對小麥產(chǎn)量（kg/ha）的影響

具體的 SAS 程序如下：

data sastjfx1;do r=1 to 4;do kind=1 to 3;do feili=1 to 3;input y@@;output;end; end; end;cards;6300 5700 4800 5400 4275 3900 4650 4800 4200 5850 5250 5700 5250 4650 3900 5850 5100 4425 5550 5550 4500 5175 3450 3300 4500 4800 3900 5900 5550 5000 5275 4125 3700 5000 4900 4175;proc sort;by kind feili r;run;ods html;proc nested;class kind feili;var y;run;proc glm;class kind feili;model y=kind feili(kind)/SS1;test H=kind E=feili(kind);means kind/snk e=feili(kind);lsmeans feili(kind)/tdiff pdiff;run;quit;ods html close;

作者單位：100850 北京，軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心

程序說明：數(shù)據(jù)步建立數(shù)據(jù)集 sastjfx1，輸入變量名和變量值，“r”、“kind”、“feili”、“y”分別代表試驗(yàn)田編號、小麥品種、肥力和小麥產(chǎn)量。sort 過程是對數(shù)據(jù)集中的變量按照主次重新排序；nested 過程用于對變量已經(jīng)按照主次排序且來自嵌套設(shè)計(jì)的平衡資料的方差分析；第 3 個(gè)過程為glm 過程，glm 過程的適用面很寬，資料可以是不平衡的，也可以不排序，僅取決于 model 語句和 test 語句的寫法，“model y = kind feili (kind)/SS1”表示因素肥力嵌套在品種里，“SS1”表示第 I 型離均差平方和（計(jì)算結(jié)果與變量在model 語句中的順序有關(guān)），不能選用 SS3；“H = kind E =feili (kind)”表示在考察品種的效應(yīng)時(shí)，以模型中肥力的效應(yīng)作為誤差項(xiàng)。程序中同時(shí)使用了 nested 過程和 glm 過程對資料進(jìn)行分析。

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

The NESTED procedure

上述結(jié)果表明：F = 2.71，P = 0.1453；F = 11.06，P <0.0001，說明就該資料而言，肥力對小麥產(chǎn)量的影響具有統(tǒng)計(jì)學(xué)意義，但是還沒有理由說明小麥品種對小麥產(chǎn)量的影響具有統(tǒng)計(jì)學(xué)意義。

The GLM procedure dependent variable: y

Source DF Type I SS Mean square F value Pr > F Kind 2 7727638.889 3863819.444 29.93 < 0.0001 Feili (kind) 6 8563750.000 1427291.667 11.06 < 0.0001 Tests of hypotheses using the type I MS for feili (kind) as an error term Source DF Type I SS Mean square F value Pr > F Kind 2 7727638.889 3863819.444 2.71 0.1453

上述結(jié)果表明：F = 15.77，P < 0.0001，說明總模型具有統(tǒng)計(jì)學(xué)意義；模型的分項(xiàng)部分只有對肥力分析的結(jié)果是正確的，其 F = 11.06，P < 0.0001，說明肥力對小麥產(chǎn)量的影響具有統(tǒng)計(jì)學(xué)意義；由 test 語句給出的結(jié)果，才是對主要因素（小麥品種）的合理檢驗(yàn)，其 F = 2.71，P = 0.1453，說明小麥品種對小麥產(chǎn)量的影響沒有統(tǒng)計(jì)學(xué)意義。

Least squares means for effect feili (kind)t for H0: LSMean (i) = LSMean (j) / Pr > |t|Dependent variable: y

以上是 9 種組合情況的小麥產(chǎn)量均值之間兩兩比較的結(jié)果，小麥品種為 A 時(shí)，由 t = –3.54241，P = 0.0015 說明高肥與低肥比較小麥產(chǎn)量之間的差別有統(tǒng)計(jì)學(xué)意義；由t = –1.52521，P = 0.1388 說明高肥與中肥比較小麥產(chǎn)量之間的差別沒有統(tǒng)計(jì)學(xué)意義；由 t = –2.01721，P = 0.0537 說明中肥與低肥比較小麥產(chǎn)量之間的差別沒有統(tǒng)計(jì)學(xué)意義。小麥品種為 B 時(shí)，由 t = –4.52642，P = 0.0001 說明高肥與中肥比較小麥產(chǎn)量之間的差別有統(tǒng)計(jì)學(xué)意義；由 t = –6.19922，P < 0.0001 說明高肥與低肥比較小麥產(chǎn)量之間的差別有統(tǒng)計(jì)學(xué)意義；由 t = –1.67281，P = 0.1059 說明中肥與低肥比較小麥產(chǎn)量之間的差別沒有統(tǒng)計(jì)學(xué)意義。小麥品種為 C時(shí)，由 t = –0.3936，P = 0.6970 說明高肥與中肥比較小麥產(chǎn)量之間的差別沒有統(tǒng)計(jì)學(xué)意義；由 t = –3.24721，P = 0.0031說明高肥與低肥比較小麥產(chǎn)量之間的差別有統(tǒng)計(jì)學(xué)意義；由t = –2.85361，P = 0.0082 說明中肥與低肥比較小麥產(chǎn)量之間的差別有統(tǒng)計(jì)學(xué)意義。

The GLM procedure

Student-newman-keuls test for y

上述結(jié)果表明：不同小麥品種之間小麥產(chǎn)量的差別沒有統(tǒng)計(jì)學(xué)意義。而小麥品種為 A時(shí)，高肥的小麥產(chǎn)量高于低肥的小麥產(chǎn)量；小麥品種為 B 時(shí)，高肥的小麥產(chǎn)量高于中肥和低肥的小麥產(chǎn)量；小麥品種為 C 時(shí)，高肥的小麥產(chǎn)量高于低肥的小麥產(chǎn)量，中肥的小麥產(chǎn)量高于低肥的小麥產(chǎn)量。

2 裂區(qū)設(shè)計(jì)及其資料的統(tǒng)計(jì)分析

2.1 裂區(qū)設(shè)計(jì)的概念及應(yīng)用場合

試驗(yàn)因素分階段進(jìn)入試驗(yàn)過程，通常由先進(jìn)入試驗(yàn)的試驗(yàn)因素（設(shè)為 A）構(gòu)成單因素多水平設(shè)計(jì)或由先進(jìn)入試驗(yàn)的試驗(yàn)因素（設(shè)為 A）與區(qū)組因素（設(shè)為 B）構(gòu)造出含 m次獨(dú)立重復(fù)試驗(yàn)的隨機(jī)區(qū)組設(shè)計(jì)；再把接受因素 A 各水平處理或接受因素 A 與因素 B 各組合水平處理的 m 個(gè)受試對象隨機(jī)地分配給在第二階段進(jìn)入試驗(yàn)的試驗(yàn)因素 C的 m 個(gè)水平，這樣安排試驗(yàn)因素的方法稱為裂區(qū)設(shè)計(jì)或分割設(shè)計(jì)。

結(jié)合實(shí)際問題，當(dāng)試驗(yàn)研究過程自然形成 2 個(gè)或多階段（有時(shí)稱為工序），各階段涉及的試驗(yàn)因素彼此不同，但需要等整個(gè)試驗(yàn)過程結(jié)束后，才能觀測定量指標(biāo)的結(jié)果，就需要用到此設(shè)計(jì)[1-2]。裂區(qū)設(shè)計(jì)的構(gòu)造方法和列表格式參見例 2。

2.2 裂區(qū)設(shè)計(jì)定量資料方差分析的 SAS 實(shí)現(xiàn)

例 2在 3 塊不同的試驗(yàn)田中進(jìn)行玉米的施肥試驗(yàn)。在苗期對每塊試驗(yàn)田使用不同的施氮量，即正常施用苗肥（78 kg/hm2）和減少 1/3 苗肥（52.5 kg/hm2）；在 5 葉期，采用灌根法每公頃用清水 600 kg，分別加入不同用量的農(nóng)夫樂，充分?jǐn)嚢韬螅种陝┝繚补嘣谟衩赘?。資料如表2所示。試分析不同的施氮量和不同的農(nóng)夫樂用量對玉米產(chǎn)量的影響。

表2 不同施氮量和不同農(nóng)夫樂施用量對玉米產(chǎn)量的影響

分析與 SAS 實(shí)現(xiàn)：該資料中，試驗(yàn)因素“不同施氮量”和區(qū)組因素在苗期先出現(xiàn)在試驗(yàn)中，每個(gè)區(qū)組被分成兩個(gè)一級單位；當(dāng)玉米生長到 5 葉期時(shí)，每個(gè)一級單位再被分為3 個(gè)二級單位，分別施加 3 種不同用量的農(nóng)夫樂，定量的觀測指標(biāo)為玉米產(chǎn)量，兩個(gè)試驗(yàn)因素分兩個(gè)階段進(jìn)入試驗(yàn)過程，所以該資料為裂區(qū)（或分割）設(shè)計(jì)。

具體的 SAS 程序如下：

data sastjfx2;do a=1 to 2;do b=1 to 3;do s=1 to 3;input x @@;output;end;end;end;cards;11.98 11.97 11.93 10.77 10.75 10.75 9.90 9.80 9.60 11.72 11.71 11.70 10.74 10.71 10.73 9.78 9.75 9.76;ods html;proc glm;class a s b;model x=a s a*s b a*b/ss3;test h=a s e=a*s;means a*b;lsmeans a*b/tdiff pdiff;run;quit;ods html close;

程序說明：數(shù)據(jù)步建立數(shù)據(jù)集 sastjfx2，輸入變量名和變量值，“a”、“b”、“s”、“x”分別代表不同施氮量、農(nóng)夫樂的不同施用量、試驗(yàn)田編號和玉米產(chǎn)量。過程步調(diào)用 glm過程，“test h = a s e = a*s;”表示在分析因素 a 和區(qū)組因素s 時(shí)，以它們的交互作用的均方作為誤差的均方（注意：若出現(xiàn)各組例數(shù)不等或某些試驗(yàn)條件下未做試驗(yàn)，應(yīng)使用 glm過程來實(shí)現(xiàn)計(jì)算）。

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

The GLM procedure dependent variable: x

以上結(jié)果表明：F = 459.74，P < 0.0001，說明總的模型具有統(tǒng)計(jì)學(xué)意義。

R-square Coeff var Root MSE x mean 0.998070 0.519693 0.056026 10.78056 Source DF Anova SS Mean square F value Pr > F a 1 0.04013889 0.04013889 12.79 0.0072 s 2 0.01471111 0.00735556 2.34 0.1581 a*s 2 0.00964444 0.00482222 1.54 0.2725 b 2 12.86831111 6.43415556 2049.82 < 0.0001 a*b 2 0.05497778 0.02748889 8.76 0.0097

以上結(jié)果表明：對 a * s 的檢驗(yàn)，F(xiàn) = 1.54，P = 0.2725，說明不同施氮量與隨機(jī)區(qū)組 s 之間的交互作用沒有統(tǒng)計(jì)學(xué)意義；對 b 的檢驗(yàn)，F(xiàn) = 2049.82，P < 0.0001，說明農(nóng)夫樂不同施用量對玉米產(chǎn)量的影響具有統(tǒng)計(jì)學(xué)意義；對 a * b 的檢驗(yàn)，F(xiàn) = 8.76，P = 0.0097，說明不同的施氮量和農(nóng)夫樂的不同施用量之間的交互作用具有統(tǒng)計(jì)學(xué)意義。

Tests of hypotheses using the anova MS for a*s as an error term Source DF Anova SS Mean square F value Pr > F a 1 0.04013889 0.04013889 8.32 0.1021 s 2 0.01471111 0.00735556 1.53 0.3960

以上結(jié)果表明：F = 8.32，P = 0.1021；F = 1.53，P =0.3960，說明不同的施氮量和不同的試驗(yàn)田對玉米產(chǎn)量的影響之間的差別沒有統(tǒng)計(jì)學(xué)意義。

Level of a Level of b N x Mean Std Dev 1 1 3 11.9600000 0.02645751 1 2 3 10.7566667 0.01154701 1 3 3 9.7666667 0.15275252 2 1 3 11.7100000 0.01000000 2 2 3 10.7266667 0.01527525 2 3 3 9.7633333 0.01527525

Least squares means for effect a*b t for H0: LSmean (i) = LSmean (j) / Pr > |t|Dependent variable: x

以上是 6 種組合條件下玉米產(chǎn)量均值之間兩兩比較的結(jié)果。對減肥區(qū)而言，由 t = –26.3053，P < 0.0001 說明農(nóng)夫樂施用量為 30 kg/hm2與農(nóng)夫樂施用量為 15 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義；由 t = –47.9471，P <0.0001 說明農(nóng)夫樂施用量為 30 kg/hm2與農(nóng)夫樂施用量為7.5 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義；由 t =–21.6418，P < 0.0001 說明農(nóng)夫樂施用量為 15 kg/hm2與農(nóng)夫樂施用量為 7.5 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義。對正常施肥區(qū)而言，由 t = –21.496，P < 0.0001 說明農(nóng)夫樂施用量為 30 kg/hm2與農(nóng)夫樂施用量為 15 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義；由 t = –42.5549，P <0.0001 說明農(nóng)夫樂施用量為 30 kg/hm2與農(nóng)夫樂施用量為7.5 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義；由 t =–21.0588，P < 0.0001 說明農(nóng)夫樂施用量為 15 kg/hm2與農(nóng)夫樂施用量為 7.5 kg/hm2之間玉米產(chǎn)量的差別有統(tǒng)計(jì)學(xué)意義。

由以上結(jié)果可以看出，減肥區(qū)和正常施肥區(qū)比較玉米產(chǎn)量之間的差別沒有統(tǒng)計(jì)學(xué)意義；但是無論是減肥區(qū)還是正常施肥區(qū)，農(nóng)夫樂施用量的不同造成的玉米產(chǎn)量之間的差別都具有統(tǒng)計(jì)學(xué)意義，玉米產(chǎn)量的高低順序?yàn)檗r(nóng)夫樂施用量30 kg/hm2> 15 kg/hm2> 7.5 kg/hm2。

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

[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:215-222. (in Chinese)胡良平. Windows SAS 6.12 & 8.0實(shí)用統(tǒng)計(jì)分析教程. 北京: 軍事醫(yī)學(xué)科學(xué)出版社, 2001: 215-222.

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

胡良平，Email：lphu812@sina.com