平衡不完全區組設計 數據分析的SAS實踐

• • 數據的定性分析：平衡不完全區組設計的參數
  • ANOVA：使用催化劑是否能夠顯著提高產率？
  • ANOVA之后的分析：哪種催化劑效果更好？

某一個藥學研究團隊想要研究四種不同的催化劑對甲硝唑產率的影響。該團隊訂購了四批相同的原材料，以原材料的批次作為區組、催化劑作為處理變量，設計并完成了一個平衡不完全區組設計，數據記錄如下：

下面將用這個例子介紹平衡不完全區組試驗（Balanced Incompete Blocking Design，BIBD）結果的數據分析方法。

數據的定性分析：平衡不完全區組設計的參數

上例的BIBD中，處理變量共有4個取值（$a=4$）、區組一共有4個（$b=4$），因此一共有16中不同的處理，但在試驗中，只進行了12組處理（$N=12$），因此叫不完全試驗。處理變量的每個取值被分配到3個不同的區組中（$r=3$），每個區組會與3個處理變量的取值構成一組處理（$k=3$），平衡的含義是
$$N=ar=bk$$

并且每一對處理變量的取值，比如第一種和第二種催化劑，它們會同時出現在2個不同的區組中（$\lambda =2$），
$$\lambda (a?1)=r(k?1)$$

下面把數據輸入到SAS中：

data example; /*example是數據表名稱*/ input trt block resp @@; /*trt block resp表示數據表中的變量，@@表示每三個值作為一組分別賦給三個變量*/ cards;/*這里也可以把cards換成datalines*/ 1 1 73 1 2 74 1 4 71 2 2 75 2 3 67 2 4 72 3 1 73 3 2 75 3 3 68 4 1 75 4 3 72 4 4 75 ;proc print data = example; run;/*輸出數據表檢查是否有誤，SAS語句總是遇到run才輸出*/

輸入到SAS中的數據表是這樣的：

ANOVA：使用催化劑是否能夠顯著提高產率？

ANOVA是分析試驗設計的首選方法，它簡單方便能夠快速得出結論，但缺點也比較明顯，需要同方差、正態性假設，并且只能判斷處理變量對試驗結果是否具有顯著影響。下面是用SAS做ANOVA的code：

proc glm data=example;/*對數據表example使用GLM procedure*/ class block trt;/*將催化劑與區組作為類型變量*/ model resp = trt block; output out=myout r=res p=pred;/*輸出到數據表myout中，并記錄下殘差和擬合值*/ run;

下表是GLM procedure輸出的第一張表，紅框內是這個GLM模型整體的p值，很顯然這個模型在0.01的水平上是顯著的，說明催化劑和區組對產率是有顯著的影響的。

下表是模型的一些整體統計，第一個R方是0.959877，非常接近1，說明這個GLM模型的解釋力非常強。較大的R方是試驗研究的特點，因為我們可以通過一些試驗設計技術（比如隨機化、區組化等）控制我們不感興趣又可能影響試驗結果的因素。

下表展示的是平方和在Model的兩個變量間的分解，Type I SS是序貫的方差分析，變量進入模型的順序影響分解結果；Type III SS是不考進入次序的平方和分解。當試驗是完全的時候，二者是相等的；但對于BIBD，二者并不相等，因為處理變量與區組是沒有次序上的先后關系，所以我們應該參考Type III SS的結果。紅框內的結果是處理變量的p值，它說明在0.05的顯著性水平下，催化劑對產率具有顯著影響。

GLM procedure還輸出了一張圖，這張圖叫interaction plot，根據這張圖，我們可以直觀判斷哪一種催化劑效果更好。從下面的圖可以看出，在不同區組中，得到更多產率的依次是第4、3、2、1中催化劑。

因為ANOVA要在一定的假設下才可以使用，所以我們需要判斷ANOVA的假設是否成立，如果不成立的話需要用一些其他的手段來彌補，比如殘差圖中出現曲線模式，可以考慮用Box-Cox變換處理催化率；正態性不成立，可以考慮用非參的Kruskal-Willis檢驗代替ANOVA進行分析。

proc sgplot data=myout;/*殘差關于擬合值的散點圖*/ scatter y=res x=pred; refline 0; run;

這個殘差圖有一點點二次的pattern，但總體上沒太大問題。

proc univariate data=myout normal;/*一元統計分析，正態性分析*/ var res;/*指定正態性分析的對象*/ qqplot res/normal(mu=est sigma=est color=red L=1);/*做QQ圖，并擬合直線*/ run;

從QQ圖可以看出，正態性基本還是成立的，下圖是更正式的正態性的檢驗，四種檢驗的p值都非常大，說明不能顯著拒絕原假設（殘差服從正態分布），因此我們可以接受ANOVA的假設。

ANOVA之后的分析：哪種催化劑效果更好？

先看下面這一段SAS代碼：

proc glm data=example; class block trt; model resp = trt block; lsmeans trt / tdiff pdiff adjust=bon stderr; lsmeans trt / pdiff adjust=tukey stderr; contrast 'a' trt 1 -1 0 0; estimate 'b' trt 1 -1 0 0; output out=myout r=res p=pred; run;

這一段代碼與上面做ANOVA的代碼相比，多了一些語句。第一類語句是lsmean，lsmean語句lsmean trt的作用是基于催化劑的效應的最小二乘估計量做推斷；第一句中的tdiff和pdiff表示對不同催化劑的效應進行兩兩比較，返回t值和p值，adjust=bon表示按Bonferroni方法對一共六組兩兩比較進行familywise調整，stderr表示返回標準誤；第二句也是做兩兩比較，但用Tukey方法進行。這兩個語句的輸出如下：
Bonferroni方法

第一張表是四種催化劑的效應的最小二乘估計量及其t檢驗，第二張表表示兩兩比較的結果，可以看出催化劑4和1、2、3的效應之間存在顯著差異，而1、2、3之間的差異不顯著。

這張圖是對上表結果的可視化，橫豎各四條灰色的線，一共有16個交點，因為只有六組兩兩比較，所以我們只需要看在對角虛線的左上部分的六個交點，最上面的三個交點表示催化劑4與另外三個催化劑的比較，三條藍線與對角線垂直，線段范圍表示兩兩比較的t檢驗導出的置信區間，藍線段與對角線不相交，說明催化劑4的效果顯著好于另外3種催化劑。三條紅線也可以類似分析。

Tukey方法
Tukey方法的結果與Bonferroni方法沒有太大區別。

在上面的代碼中，第二類語句是contrast和estimate。contrast ‘a’ trt 1 -1 0 0;表示以1 -1 0 0為系數構造contrast，將其命名為a，估計contrast的效應，并做推斷：

estimate ‘b’ trt 1 -1 0 0;表示以1 -1 0 0為系數構造效應的線性組合，并對其進行推斷：

這兩個語句的共同點是都在比較催化劑1和催化劑2的效應，并且取得的相同的p值，發現兩種催化劑對產率的作用沒有顯著差別。

總結

以上是生活随笔為你收集整理的平衡不完全区组设计 数据分析的SAS实践的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。