隨機(jī)效應(yīng)模型用于污水?dāng)?shù)據(jù)分析

1、隨機(jī)效應(yīng)模型用于污水?dāng)?shù)據(jù)分析 摘 要：縱向數(shù)據(jù)分析在生物，醫(yī)學(xué)，經(jīng)濟(jì)，金融等領(lǐng)域具有廣泛的應(yīng)用，是近年來(lái)統(tǒng)計(jì)學(xué)的熱點(diǎn)課題之一，隨機(jī)效應(yīng)模型是縱向數(shù)據(jù)分析中最常用的模型之一，在過(guò)去的十年中對(duì)縱向數(shù)據(jù)隨機(jī)效應(yīng)模型的統(tǒng)計(jì)分析已經(jīng)引起統(tǒng)計(jì)學(xué)家的關(guān)注。在本文中應(yīng)用固定效應(yīng)模型、隨機(jī)效應(yīng)模型分析一組污水氨氮含量的數(shù)據(jù)，并分析組件交互作用。 關(guān)鍵詞：縱向數(shù)據(jù)；固定效應(yīng)；隨機(jī)效應(yīng)；aic；bic；氨氮 縱向數(shù)據(jù)longitudinal data（或panel data）是用來(lái)描述一個(gè)總體中給定樣本在一段時(shí)間的情況，并對(duì)樣本中每一個(gè)樣本單位都進(jìn)行多重觀察。這種多重觀察既包括對(duì)樣本單位在某一時(shí)刻上多個(gè)特性進(jìn)行觀察

2、，也包括對(duì)該樣本單位的這些特性在一段時(shí)間的連續(xù)觀察，連續(xù)觀察將得到數(shù)據(jù)集稱(chēng)為縱向數(shù)據(jù)。而時(shí)間序列數(shù)據(jù)是一維數(shù)據(jù)，是變量按時(shí)間得到的數(shù)據(jù)，截面數(shù)據(jù)是變量在截面空間上的數(shù)據(jù)。數(shù)據(jù)是同時(shí)在時(shí)間和截面空間上取得的二維數(shù)據(jù)。從橫截面（cross section）上看，是由若干個(gè)體在某一時(shí)刻構(gòu)成的截面觀測(cè)值，從縱剖面上看是一個(gè)時(shí)間序列。時(shí)間序列數(shù)據(jù)與縱向數(shù)據(jù)的顯著差別是時(shí)間序列數(shù)據(jù)研究的是很少的個(gè)體（通常只有一個(gè)個(gè)體）和重復(fù)測(cè)量許多次，且數(shù)據(jù)具有序列相關(guān)性。而縱向數(shù)據(jù)研究的是很多個(gè)體，重復(fù)測(cè)量次數(shù)并不多。但它們存在共性，重復(fù)測(cè)量是相關(guān)的，且相鄰數(shù)據(jù)相關(guān)性很大，而離得較遠(yuǎn)的數(shù)據(jù)也就是隨著時(shí)間間隔增大，相關(guān)性

3、減小?？v向數(shù)據(jù)通常用雙下標(biāo)或多下標(biāo)變量表示。某城市家庭及住宿業(yè)的污水中的氨氮含量數(shù)據(jù)，每個(gè)個(gè)體測(cè)量了三天每天測(cè)量三次。這里 表示第個(gè)i個(gè)體第t個(gè)時(shí)刻的氨氮含量。 一、固定效應(yīng)模型 在現(xiàn)實(shí)世界中，存在著大量的這樣的情況：兩個(gè)變量x，y有一些依賴(lài)關(guān)系。比如身高和體重，身高x大時(shí)的體重y也傾向于大的，但是x也不能?chē)?yán)格的決定y值。變量之間的這種關(guān)系稱(chēng)為“相關(guān)關(guān)系”，回歸模型是研究相關(guān)關(guān)系的一個(gè)有力工具。我們?cè)O(shè)想y的值由兩部分組成：一部分是由x能夠決定的部分，它是x的函數(shù)，這個(gè)函數(shù)關(guān)系是線性的或者是近似線性的 （1） 這里的參數(shù)0，1分別為固定的截距和斜率，yij代表第i個(gè)個(gè)體第j個(gè)時(shí)刻的觀測(cè)，xij代

4、表協(xié)變量，ij表示隨機(jī)誤差，并且認(rèn)為它服從正態(tài)分布。我們稱(chēng)這個(gè)線性回歸模型為固定效應(yīng)模型。在這里 ， 測(cè)量值的方差為2，協(xié)方差為0。 我們看到同一個(gè)體的不同時(shí)刻的測(cè)量值之間是不相關(guān)的，但是這與實(shí)際情況不符，當(dāng)我們對(duì)總體的均值感興趣，并且不關(guān)心個(gè)體之間差別的時(shí)候可以應(yīng)用這個(gè)模型。但是有些情況我們不但關(guān)心總體的均值，而且關(guān)心個(gè)體之間的差別，因此我們將模型改良引入下面的隨機(jī)截距模型。 二、隨機(jī)截距模型 進(jìn)一步考查簡(jiǎn)單模型（1），引進(jìn)隨機(jī)效應(yīng)項(xiàng)。則模型為 且 因?yàn)?這里xij=1，1是固定效應(yīng)截距項(xiàng)，1+bi為第i個(gè)個(gè)體的均值。因?yàn)殡S機(jī)效應(yīng)bi的均值為零，那么bi代表第i個(gè)個(gè)體的截距與總體截距的偏差。

5、yij的邊緣協(xié)方差由它與邊緣均值的偏差決定。對(duì)有隨機(jī)截距項(xiàng)的模型來(lái)說(shuō)邊緣方差為 同一個(gè)個(gè)體不同時(shí)刻的協(xié)方差為 因?yàn)閎i與ij相互獨(dú)立的，所以我們看到邊緣協(xié)方差陣有對(duì)稱(chēng)的形式 由此我們得到同一個(gè)個(gè)體不同時(shí)刻響應(yīng)變量的相關(guān)系數(shù)相關(guān)系數(shù)反映了隨機(jī)截距模型的一個(gè)重要方面：隨機(jī)效應(yīng)項(xiàng)bi可以看成重復(fù)測(cè)量值間的相關(guān)量。我們看到協(xié)方差陣為對(duì)稱(chēng)陣，這就意味著測(cè)量值之間的相關(guān)性不隨時(shí)間間隔的增加而變化，而實(shí)際情況是測(cè)量值之間的相關(guān)性是隨時(shí)間間隔的增加而減小的。所以隨機(jī)截距模型可以進(jìn)一步改進(jìn)，比如增加時(shí)間效應(yīng)。 用向量、矩陣來(lái)表示這個(gè)模型則有 這里 為（p *1）的固定效應(yīng)向量，bi是（q*1）的隨機(jī)效應(yīng)向量，

6、xi是（）的協(xié)變量陣，zi是一個(gè)（）的協(xié)變量陣（）。模型中的隨機(jī)效應(yīng)項(xiàng)bi與測(cè)量誤差ei都看作是服從正態(tài)分布的。則有 響應(yīng)yij的邊緣均值， 條件均值 方差 若 是一個(gè)對(duì)角矩陣，則 顯然具有非零的非對(duì)角元，這樣我們就看到了縱向數(shù)據(jù)研究中同一個(gè)體重復(fù)測(cè)量間的相關(guān)性。 三、污水?dāng)?shù)據(jù)處理 （一）9次重復(fù)測(cè)量 1.隨機(jī)截距模型。將所有數(shù)據(jù)拉直后進(jìn)行分析，每個(gè)個(gè)體測(cè)量了9次，對(duì)數(shù)據(jù)開(kāi)平方使之接近正態(tài)分布如表2，應(yīng)用隨機(jī)截距模型來(lái)擬和數(shù)據(jù)，個(gè)體隨機(jī)效應(yīng)為變量“來(lái)源”。 則模型為 i=1，2 j=120 k =19 yijk表示來(lái)自第。分別表示組效應(yīng)，個(gè)體效應(yīng)，時(shí)間效應(yīng)；為組與時(shí)刻的交互；ij表示組與個(gè)體的

7、交互，ijk表示隨機(jī)誤差。在這里且相互獨(dú)立。 應(yīng)用sas分析得到如下結(jié)果 效應(yīng) f值 prf group 3.09 0.0809 t 1.63 0.122 group*t 6.37 cov parm cov p1 cov p2 個(gè)體 0.8854 -0.00168 殘差 -0.00168 0.05564 從分析結(jié)果中看到組間差異近似顯著但是時(shí)刻差異不顯著，組與時(shí)刻交互非常顯著。隨機(jī)效應(yīng)方差的極大似然估計(jì)為： 2.5994 標(biāo)準(zhǔn)差為0.941=sqrt0.8854；測(cè)量方差： 2.0015標(biāo)準(zhǔn)差 0.2359=sqrt0.05564 。 二者之間的協(xié)方差為：-0.00618， 近似于零， 另外，

8、 2.5994/0.941=2.761.96， 2.0015/0.2359=8.481.96， 所以拒絕這些方差為0的假設(shè)。 aic （smaller is better） 663.0 bic （smaller is better） 665.0） 2.固定效應(yīng)模型。如果不考慮隨機(jī)截距項(xiàng)， 僅擬和固定效應(yīng)的模型。 yijk表示來(lái)自第。分別表示組效應(yīng)，個(gè)體效應(yīng)，時(shí)間效應(yīng)；為組與時(shí)刻的交互；ijk表示隨機(jī)誤差。 有如下結(jié)果： 效應(yīng) f prf group 17.06 t 0.71 0.6845 group*t 2.93 0.0044 aic （smaller is better） 750.1 bic

9、 （smaller is better） 762.0） 這時(shí)組間差異顯著，時(shí)刻差異仍不顯著，組與時(shí)刻交互仍顯著。但是比較兩者的aic&bic，傾向于選擇帶有隨機(jī)截距想的模型. 生活污水和住宿業(yè)污水氨氮的含量一天中隨不同的時(shí)刻變化，但是兩者的高峰期不一樣。另外，這種隨時(shí)刻的變化似乎和天沒(méi)有關(guān)系，。所以將一天中的時(shí)刻1， 時(shí)刻2， 時(shí)刻3合并起來(lái)進(jìn)行分析。 （二）按天3次重復(fù)測(cè)量 1.隨機(jī)效應(yīng)模型 i=1，2 j=120 k =1，2，3 yijk表示來(lái)自第。分別表示組效應(yīng)，個(gè)體效應(yīng)，時(shí)間效應(yīng)；為組與時(shí)刻的交互；ij表示組與個(gè)體的交互，ijk表示隨機(jī)誤差。在這里 且相互獨(dú)立。 sas分析結(jié)果以及分

10、組圖像如下 效應(yīng) f prf group 3.42 0.0725 day 1.52 0.2314 group*day 0.49 0.6179 cov parm cov p1 cov p2 個(gè)體 8.1103 -0.1066 殘差 -0.1066 0.3199 這里組間的差異不明顯顯著而時(shí)間的差異，組與時(shí)間的交互都不顯著。從結(jié)果上看污水中氨氮含量與天沒(méi)有關(guān)系。那么時(shí)刻對(duì)污水氨氮含量的影響是怎樣的呢？下面再將時(shí)刻1，2，3按照天合并，即為一天中三個(gè)時(shí)刻。分析時(shí)刻之間的差別。 三、按時(shí)刻3次重復(fù)測(cè)量 模型形式與上一目相同，應(yīng)用sas分析結(jié)果以及分組圖像如下。 效應(yīng) f prf group 3.20

11、0.0820 time 2.02 0.1477 group*time 14.78 cov parm cov p1 cov p2 個(gè)體 8.1640 -0.2407 殘差 -0.2407 0.7222 從結(jié)果看出， 時(shí)刻差異不顯著， 組間差異也不是一個(gè)非常顯著的變量， 但是它們的交互項(xiàng)非常顯著（p1.96。所以拒絕方差為零的假設(shè)。 aic （smaller is better） 275.5 bic （smaller is better） 277.5 天宇國(guó)際 天宇西域 雪林家園 天寶小區(qū) 雨荷家園 伊諾小區(qū) 地震局 科研貝貝小區(qū) 佳霖飯店 茉莉花飯店 元寶山賓館 開(kāi)心賓館 可可同賓館 青山賓館

12、口口鄉(xiāng)賓館 亞明賓館 如家賓館 遠(yuǎn)東賓館 林麟賓館 源和賓館 如果將個(gè)體效應(yīng)看成固定的則有如下結(jié)果 效應(yīng) f prf group 7.49 0.0084 time 0.67 0.5182 group*day 4.87 0.0113 aic （smaller is better） 298.1 aicc （smaller is better） 298.2 相比之下增加隨機(jī)效應(yīng)項(xiàng)的結(jié)果更好， 另外aic &bic也比有隨機(jī)效應(yīng)項(xiàng)的模型大。所以在某城市污水氨氮含量數(shù)據(jù)分析的過(guò)程中，發(fā)現(xiàn)隨機(jī)效應(yīng)模型分析縱向數(shù)據(jù)優(yōu)于固定效應(yīng)模型。 參考文獻(xiàn)： 1陳希孺，陳桂景，吳啟光，趙林城，線性模型的參數(shù)估計(jì)理論 m，科學(xué)出版社，北京，1985. 2王松桂，史建紅，尹素菊，吳密霞，線性模型引論，科學(xué)出版社，北京，2004. 3roderick j. a little & donald b. rubin 著，孫山譯，耿直審，缺失數(shù)據(jù)統(tǒng)計(jì)分析美，1987. 4茆詩(shī)松主編，高等數(shù)理統(tǒng)計(jì)，施普林格出版社. 5孫燕，柴根象，縱向數(shù)據(jù)的隨機(jī)效應(yīng)模型，同濟(jì)大學(xué)