好久沒有更新文章了，因為同學們咨詢的問題有點多，另一個原因就是自己實在太懶。。。。
今天繼續給大家寫廣義估計方程式GEE。
In statistics, a generalized estimating equation (GEE) is used to estimate the parameters of a generalized linear model with a possible unknown correlation between outcomes.
上面的英文解釋告訴我們，當我們想用廣義線性模型的時候，突然發現我們的結局變量恐怕是有某種相關性的（比如重複測量，比如嵌套），我們就可以考慮廣義估計方程了。
執行個體說明
還是來看例子：現在我手上有來自兩個病區患有呼吸系統疾病共111個病人，分為治療組和安慰劑組，進行縱向随訪4次，每一次随訪我們會記錄病人情況（1 = good, 0 = poor）
于是得到我們的資料，大概長這樣：
其中，我們的響應變量有4個，因為我們随訪了4次嘛，拿到上面的資料，我現在問你治療到底有沒有作用？
你該如何分析上面的資料？
首先我們明确我們的研究問題是治療效果，就是說我們想看治療組A和安慰劑組P是不是在響應變量上有差異。
我們遇到的問題在于響應變量有4個，怎麼辦呢？
反正你肯定知道隻選其中一個響應進行分析肯定是不對的，把4次響應求均值肯定也是不對的，4次響應都當作獨立的觀測肯定也是不對的。
我們必須想辦法利用全部的響應變量，而且還需要考慮同一個人的4次響應之間的相關性。
嘿嘿，怎麼做呢？
混合效應模型？
可以，不過今天我們的主角是廣義估計方程式。
廣義估計方程和混合模型的根本差別就在于它把并不是很關注組内的相關，而是把重點放在了對效應均值的模組化上，為了更好地了解GEE我們先看兩個概念：
The very crux of GEE is instead of attempting to model the within-subject covariance structure, to treat it as a nuisance and simply model the mean response.
方差函數與均值函數
假設我們有獨立響應變量
與自變量矩陣
，且
與
有如下關系：
其中
被稱作均值函數，你把這個μ(⋅)想辦法移到等号左邊它就成了GLM中的連接配接函數。同時，我們假設
稱
為方差函數。
那麼，我們的主要目标就是估計β向量。大家也可以看出來相比廣義線性模型，廣義估計方程式，就是多了方差函數。
邊際模型與條件模型
這個有點難哦，因為沒怎麼學過矩陣運算，是以我隻能嘗試性地給大家解釋一波，不一定哈。
之前有給大家寫很多混合模型的文章，回憶一下混合模型是如何控制組間相關性的，是通過估計随機效應實作的：
混合模型的模型表達是這樣：
有固定效應β，有随機效應u，那麼這麼一個模型也可以看作是一個條件模型。
怎麼講？
就是如果我們令殘差e的變異為σ 2 I，在随機效應u存在的情況下，y的分布為：
這個應該沒問題吧，因為我們已經考慮了u，考慮u之後我們認為方差齊的情況下，因變量就是獨立的啦，因為混合效應的假設之一就是随機效應服從正态分布嘛。
再看邊際模型
所謂邊際模型在我的了解中就是把條件去掉，哈哈，簡單粗暴的了解（不一定對哈）。
在條件去掉之後我還是希望我的因變量y是互相獨立的，這時我們就需要給我們的條件，也就是随機效應u加一個變異Var(u) = G，那麼此時y的邊際分布就為：
上面式子啥意思？
上面的式子告訴我們y和V = ZGZ0 + σ 2 I是有很大關系的，y的變化是依賴于固定效應以V為标準差進行變化的。
我們把上面的式子簡化一下：
就是y的分布如上式，y是依賴于固定效應後以V為标準差進行變化的，是以我們要對y進行模組化我們隻需要設定好V，然後直接估計β就行？GEE就是這麼一個思想。
寫到這再給大家聊聊GEE和混合模型的差別：
GEE估計的整體樣本的均值，而混合模型也會對個體變異進行模組化，一個是population average，一個是subject specific，這兒給大家推薦一篇文獻，大家自己去琢磨琢磨（文獻裡面有寫什麼時候用GEE什麼時候用混合模型哦）
Szmaragd, C., Clarke, P., & Steele, F. (2013). Subject specific and population average models for binary longitudinal data: a tutorial. Longitudinal and life course studies, 4(2), 147-165.
邊際模型的建構
對于上面提到的響應y有相關的資料，我們可以通過邊際模型進行處理，我們需要設定的參數有兩個：
以混合模型的思想了解，第一個就是固定效應如何，第二個就是變異如何。
做GEE的時候，我們認為不同個案之間的測量是獨立的，同一個個案的不同測量是可以有相關的，和混合模型一樣，隻有systematic part（固定效應）才是我們關注的點，而Random part（随機效應）指數為了控制住統一觀測不同測量之間的相關，這一部分是服務于固定效應估計的，并不是我們關注的點。（看下面的圖檔）
在GEE中對随機效應的控制Random part是通過設定作業相關矩陣實作的。
為啥叫作業相關矩陣呢，其實我們本來想要的是不同時點響應y的真實相關，但是這個不好得到，是以我們為了工作的繼續進行，自己進行自認為合理的相關矩陣設定，是以叫作業相關矩陣，那麼這個相關矩陣肯定就可以有很多種了，我們用例子資料一個一個來看：
一階自相關AR(1)
在我們的例子中y是測量了4次的，我現在假定測量之間的間隔時間越遠其相關性α越低，那麼我們的時點t和時點k的相關系數可以表示如下：
對于我們的例子，因為隻有4個時點我們完全可以寫出完整的相關結構矩陣：
這個矩陣依然需要我們從資料中進行估計，比如我估計出來α為0.61，那麼在的情況下我們的相關結構矩陣為下圖：
一階自相關的意思就是隔得越遠，相關越小。
等相關或者叫可交換的相關（exchangeable）
大家肯定會說，你剛剛隻是假定測量之間的間隔時間越遠其相關性α越低，那麼我是不是可以假定所有時間點的相關性都一樣呢？（能問出這個問題的同學一看就是動腦子了）
當然可以啦，比如我就假定所有時點的相關都是相等的，為α：
那麼，此時我們假定α都一樣時從資料中估計出來的α為0.46。
非确定相關（unstructured）
相應地，你還可以假定每個不同時間點α都不一樣，那麼此時的相關結構就可以如下表示：
獨立（independent）
還有一種理想情況就是重複測量之間是獨立的：
作業相關矩陣怎麼選？
因為我們拿到資料的時候，我隻知道這是一個縱向資料，我要用GEE，但是我肯定不是很确定作業相關矩陣到底是怎樣的。
是以我們需要根據理論預先設定好，或者就是根據資料試，一個一個試
比如你測人的道德水準，一個月測一次，測了4次，你完全可以假定等相關嘛，因為短短4個月我們認為道德水準是不會變化的。
如果要根據資料試，我們就有5種選擇（見下圖）：1.所有的觀測都是獨立的(普通回歸做就行)，和4種不同的工作矩陣結構。
對于每種選擇，我們都相應地進行模型的系數估計，就可以得到下表的結果：
大家可以很清晰的看到：當我們将作業相關矩陣設定為獨立（independent）和等相關（exchangeable）時，模型得到的系數估計就和混合模型的系數一樣了。但是GEE的系數标準誤要大一點。
這兒提醒大家的是，作業相關矩陣基本上根據常識設就行，不用糾結，影響不大的。
如何在R種做GEE
這部分依然是用文章開頭的例子給大家說明，來自兩個病區（center=1，2）患有呼吸系統疾病共111個病人，分為治療組(A)和安慰劑組（P），進行縱向随訪4次，每一次随訪我們會記錄病人情況（1 = good, 0 = poor）
我們收來的資料大概長這樣：
其中，我們的響應變量有4個，因為我們随訪了4次嘛，但是我們的協變量center，treat，sex,age都是不随時間變化的。
如果我們不考慮個體之内的相關性，我們可以把上面的資料轉化為長形資料直接做個邏輯斯蒂回歸：
m.glm <- glm(outcome ~ baseline + center + sex + treat + age + I(age^2),
data = respiratory, family = binomial) 當然，上面的方法肯定是不對的嘛，因為沒有考慮組内相關嘛。
我們可以考慮用GEE，需要用到的函數為geeglm，用法和glm其實差不多，我們需要通過family參數設定連接配接函數和方差函數，還需要通過corstr參數手動設定作業相關矩陣，參數的說明和取值見下圖：
我們不是要跑5個模型嘛，我這人給大家跑兩個示例，一個是作業相關矩陣是exchangeable，另一個是independence，大家可以自行去跑别的模型：
m.ex <- geeglm(outcome ~ baseline + center + sex + treat + age + I(age^2),
data = respiratory, id = interaction(center, id),
family = binomial, corstr = "exchangeable") m.ind <- geeglm(outcome ~ baseline + center + sex + treat + age + I(age^2),
data = respiratory, id = interaction(center, id),
family = binomial, corstr = "independence") 你就按照上面的方法把所有的模型跑完，然後就可以得到模型系數如下：
可以看到，基本上混合模型和不同作業相關矩陣GEE跑出來的模型系數都是差不多的
是以你完全不用糾結作業矩陣到底選哪個，如果你願意，你可以用anova函數對不同作業相關矩陣的模型進行比較。一般情況下我們在自變量變動的時候才進行模型比較，改變自變量後比較兩個模型的拟合情況，可以用anova函數，比如我想把age這個變量去掉：
m.ex.0 <- update(m.ex, . ~ . - age - I(age^2)) 然後把這個模型和全模型m.ex進行對比，我就可以用如下代碼：
anova(m.ex, m.ex.0) 從結果可以看出來，age還是對拟合模型有作用的，删除不得。
小結
今天給大家寫了廣義估計方差的粗淺知識和做法示例，感謝大家耐心看完，自己的文章都寫的很細，代碼都在原文中，希望大家都可以自己做一做，請關注後私信回複“資料連結”擷取所有資料和本人收集的學習資料。如果對您有用請先收藏，再點贊轉發。
也歡迎大家的意見和建議，大家想了解什麼統計方法都可以在文章下留言，說不定我看見了就會給你寫教程哦。
如果你是一個大學大學生或研究所學生，如果你正在因為你的統計作業、資料分析、論文、報告、考試等發愁，如果你在使用SPSS,R，Python，Mplus, Excel中遇到任何問題，都可以聯系我。因為我可以給您提供好的，詳細和耐心的資料分析服務。
如果你對Z檢驗，t檢驗，方差分析，多元方差分析，回歸，卡方檢驗，相關，多水準模型，結構方程模型，中介調節，量表信效度等等統計技巧有任何問題，請私信我，擷取詳細和耐心的指導。
If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #reports, #composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.
Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??
Then Contact Me. I will solve your Problem...
本文參考文獻：Halekoh, U., Højsgaard, S., & Yan, J. (2006). The R package geepack for generalized estimating equations. Journal of statistical software, 15(2), 1-11.
掃描二維碼擷取
更多精彩