轉(zhuǎn)自個人微信公粽號【易學(xué)統(tǒng)計】的統(tǒng)計學(xué)習(xí)筆記：R語言：廣義估計方程(GEE)
01 解決何種問題
在臨床研究中骡技，經(jīng)常會比較兩種治療方式對患者結(jié)局的影響逆济，并且多次測量結(jié)局。例如為了研究兩種降血糖藥對血糖的控制結(jié)果是否存在差異观游，研究者會在兩組人群服藥后不同的時間點記錄血糖值，然后評價降血糖結(jié)果驮俗。
為評價抗癲癇藥物的作用懂缕，觀察并記錄兩組不同用藥的人群在8周內(nèi)，每2周發(fā)病的次數(shù)王凑，分析該藥物是否有抑制癲癇發(fā)作的作用搪柑。
另外吮蛹，有的研究只收集一個時間點的數(shù)據(jù)，但一個研究者會提供多個部位的數(shù)據(jù)點拌屏，例如研究者想評價冠心病患者在冠脈搭橋術(shù)后應(yīng)用阿司匹林是否能有效降低患者血管再堵塞的風(fēng)險潮针，而一個患者可能在同一次手術(shù)中對多條冠脈動脈進(jìn)行搭橋。
這類數(shù)據(jù)存在相關(guān)性倚喂，不滿足線性模型和廣義線性模型所要求的的獨立性每篷，且計數(shù)類數(shù)據(jù)，不滿足正態(tài)性端圈，怎么辦焦读？
02 方法說明
廣義估計方程GEE是在傳統(tǒng)的模型基礎(chǔ)上對相關(guān)性數(shù)據(jù)做了校正，可以擬合Logistic回歸舱权、泊松回歸矗晃、一般線性回歸等回歸模型。專門用于處理這種重復(fù)測量資料且非獨立宴倍，不同于重復(fù)測量方差分析(要求結(jié)局變量為連續(xù)性數(shù)據(jù)且滿足球形假設(shè)檢驗)张症，既可以處理連續(xù)性結(jié)局變量也可以處理分類型結(jié)局變量，其利用連接函數(shù)將多種分布的因變量擬合為相應(yīng)的統(tǒng)計模型鸵贬。
該模型的一個重要概念是作業(yè)相關(guān)矩陣俗他，表示因變量的各次重復(fù)測量值兩兩之間的相關(guān)性的大小，常用Ri(a)表示阔逼，是個t*t維的對角陣兆衅，t是總測量次數(shù)，其第s行t列的元素表示嗜浮，Yis和Yit的相關(guān)羡亩，近似的表示個體之間平均的相關(guān)。常見的形式有這幾種危融，獨立畏铆，等相關(guān)，相鄰相關(guān)专挪，自相關(guān)和不確定相關(guān)及志。
其中獨立是指同一患者貢獻(xiàn)的各數(shù)據(jù)點數(shù)據(jù)彼此獨立，無相關(guān)性寨腔。等相關(guān)是各數(shù)據(jù)點的相關(guān)性均相等速侈，相鄰相關(guān)是僅臨近的M+1個數(shù)據(jù)有相關(guān)，自相關(guān)是用于不同時間點的數(shù)據(jù)迫卢，相鄰時間點相關(guān)性最大倚搬，時間間隔越大相關(guān)性越小，不確定相關(guān)是不限定相關(guān)結(jié)構(gòu)乾蛤，由數(shù)據(jù)本身決定每界。
GEE模型估計涉及到的數(shù)學(xué)原理比較復(fù)雜捅僵，讀者可參考相關(guān)文獻(xiàn)。下面講述GEE在R中的操作眨层。
首先我們以第三個案例作為例子庙楚，講解在R中的操作：
03 加載數(shù)據(jù)
表中變量的含義：
04 R代碼
library(gee)
dt <- read.csv('C:\\GEE.csv',stringsAsFactors = F)
fit <- gee(Outcome~Treatment+Sex,id=id,data=dt,corstr = 'exchangeable',family = 'binomial')
summary(fit)
##返回結(jié)果
##Model:
## Link: Logit
## Variance to Mean Relation: Binomial
## Correlation Structure: Exchangeable
##Call:
##gee(formula = Outcome ~ Treatment + Sex, id = id, data = dt,
## family = "binomial", corstr = "exchangeable")
##Summary of Residuals:
## Min 1Q Median 3Q Max
##-0.7573351 -0.4696669 0.2426649 0.4005292 0.5759983
##
##Coefficients:
## Estimate Naive S.E. Naive z Robust S.E. Robust z
##(Intercept) 1.1381242 0.5546159 2.0520945 0.3807494 2.9891690
##Treatment -1.0767769 0.5992072 -1.7970025 0.5388111 -1.9984310
##Sex -0.3677144 0.5885404 -0.6247904 0.5868946 -0.6265424
##
##Estimated Scale Parameter: 1.089952
##Number of Iterations: 2
##
##Working Correlation
## [,1] [,2] [,3]
##[1,] 1.0000000 -0.1617976 -0.1617976
##[2,] -0.1617976 1.0000000 -0.1617976
##[3,] -0.1617976 -0.1617976 1.0000000
05結(jié)果解讀
1.R中加載gee包，用該包中的gee()函數(shù)建模趴樱。
2.該數(shù)據(jù)集有39個觀測值馒闷，每名患者貢獻(xiàn)的數(shù)據(jù)量不一致，2號患者對兩根血管進(jìn)行搭橋手術(shù)叁征，1號患者只對一根血管行搭橋手術(shù)纳账。
3.gee()構(gòu)建模型的寫法跟之前線性回歸的寫法是一樣的，需要注意的是id用來設(shè)置數(shù)據(jù)集的標(biāo)識變量捺疼，即觀測個體疏虫，這里表示每個患者。corstr參數(shù)是指數(shù)據(jù)的相關(guān)結(jié)構(gòu)啤呼，這里設(shè)置的是等相關(guān)卧秘，還可以設(shè)置其他類型。
4.返回的結(jié)果中#model是采用logit轉(zhuǎn)換媳友，對應(yīng)的二分類斯议。#coefficients表示系數(shù)估計，其結(jié)果解釋和LR回歸一致醇锚。Treatment的系數(shù)是-1.077，exp(B)為0.314坯临，即OR值焊唬，說明術(shù)后使用阿司匹林的患者發(fā)生血管再堵塞的風(fēng)險是使用安慰劑組的0.314倍，具有統(tǒng)計學(xué)意義看靠，阿司匹林具有保護(hù)作用赶促。而納入性別變量的校正OR值是0.692。#working correlation是作業(yè)相關(guān)矩陣挟炬，對角線上都是1鸥滨，因為選的是等相關(guān)，除了對角線上谤祖，其他相關(guān)系數(shù)都是-0.162婿滓。
再以一個結(jié)局為連續(xù)變量的案例作為例子，講解在R中的操作粥喜。
06R代碼
library(gee)
head(warpbreaks)
## breaks wool tension
##1 26 A L
##2 30 A L
##3 54 A L
##4 25 A L
##5 70 A L
##6 52 A L
table(warpbreaks$wool,warpbreaks$tension)
## L M H
## A 9 9 9
## B 9 9 9
##等相關(guān)
summary(gee(breaks~tension,id=wool,data=warpbreaks,corstr = 'exchangeable'))
#自相關(guān)
summary(gee(breaks~tension,id=wool,data=warpbreaks,corstr = 'AR-M',Mv=1))
##返回結(jié)果1
##Model:
## Link: Identity
## Variance to Mean Relation: Gaussian
## Correlation Structure: Exchangeable
## coefficients:
## Estimate Naive S.E. Naive z Robust S.E. Robust z
##(Intercept) 36.38889 3.069434 11.855246 5.774705 6.301428
##tensionM -10.00000 3.910008 -2.557539 7.463905 -1.339781
##tensionH -14.72222 3.910008 -3.765266 3.731952 -3.94491
##返回結(jié)果2
##Model:
## Link: Identity
## Variance to Mean Relation: Gaussian
## Correlation Structure: AR-M , M = 1
#Coefficients:
## Estimate Naive S.E. Naive z Robust S.E. Robust z
##(Intercept) 36.41361 2.889306 12.602890 5.788146 6.291065
##tensionM -10.09043 4.077388 -2.474730 7.513413 -1.342989
##tensionH -14.72471 4.086060 -3.603646 3.758972 -3.917218
07 結(jié)果解讀
1.table()函數(shù)觀察數(shù)據(jù)集的基本結(jié)構(gòu)凸主，該數(shù)據(jù)集是R包中自帶數(shù)據(jù)集。觀察兩種羊毛A和B分別在三種不同強(qiáng)度下额湘，即L,M,H卿吐，對應(yīng)低中高三種韌度旁舰，在紡織中崩斷的次數(shù)(breaks)。每種羊毛有27個觀測值嗡官，共計54個觀測箭窜。
2.gee()構(gòu)建模型的寫法跟之前線性回歸的寫法是一樣的，需要注意的是id用來設(shè)置數(shù)據(jù)集的標(biāo)識變量衍腥，即觀測個體磺樱，這里表示每種羊毛。corstr參數(shù)同上紧阔。
3.返回的結(jié)果1中#model是Identity表示原始數(shù)據(jù)沒有進(jìn)行任何轉(zhuǎn)化坊罢，擬合線性回歸模型。#coefficients表示系數(shù)估計擅耽，其結(jié)果解釋和Lm回歸一致活孩，tensionM的系數(shù)是-10，相對于強(qiáng)度低L的羊毛乖仇，強(qiáng)度M要少斷10次憾儒，tensionH的系數(shù)是-14.7，相對于強(qiáng)度低L的羊毛乃沙，強(qiáng)度H要少斷14.7次起趾。
4.返回的結(jié)果2是改變數(shù)據(jù)相關(guān)結(jié)構(gòu)類型為自相關(guān)MR，此時一定要設(shè)置Mv參數(shù)警儒，設(shè)置為1训裆，對模型進(jìn)行估計，可看到系數(shù)結(jié)果和等相關(guān)基本一致蜀铲。
08總結(jié)
GEE具有非常好的特性边琉，既能處理連續(xù)型結(jié)局變量，又能處理分類型結(jié)局變量记劝。它的優(yōu)勢在于：
1.建模穩(wěn)定变姨，即使設(shè)置的數(shù)據(jù)相關(guān)結(jié)構(gòu)與實際有偏差，但在樣本量足夠大的時候厌丑，其的參數(shù)估計仍然具有無偏性定欧。其自變量的系數(shù)估計準(zhǔn)確性論證高于多水平模型。
2.充分利用資料信息怒竿。對于多次重復(fù)測量的數(shù)據(jù)砍鸠，充分利用每次測量結(jié)果，減少信息損失愧口。
以上就是本次分享的內(nèi)容了睦番。后面還有更多高分統(tǒng)計方法分享，請持續(xù)關(guān)注哦~
如果您覺得有用，請點贊托嚣，轉(zhuǎn)發(fā)哦~
更多統(tǒng)計小知識巩检，請關(guān)看 公粽號 易學(xué)統(tǒng)計
更多閱讀
R語言|兩因素重復(fù)測量方差分析
R語言|基于Cox模型pec包深度驗證
R語言|中位生存時間列線圖繪制
R語言|Cox模型校準(zhǔn)度曲線繪制
R語言|中位生存時間列線圖繪制
基于Lasso回歸篩選變量構(gòu)建Cox模型并繪制Nomogram
R語言Logistic回歸模型驗證及Nomogram繪制
如何進(jìn)行高維變量篩選和特征選擇(一)？Lasso回歸
