為什麼 SAS PROC GLIMMIX 給我的隨機斜率與二項式 glmm 的 glmer (lme4) 非常不同
我是一個更熟悉 R 的用戶,並且一直試圖在 5 年內為四個棲息地變量估計大約 35 個人的隨機斜率(選擇係數)。響應變量是位置是“使用”(1) 還是“可用”(0) 棲息地(以下“使用”)。
我使用的是 Windows 64 位計算機。
在 R 版本 3.1.0 中,我使用下面的數據和表達式。PS、TH、RS 和 HW 是固定效應(標準化、測量到棲息地類型的距離)。lme4 V 1.1-7。
str(dat) 'data.frame': 359756 obs. of 7 variables: $ use : num 1 1 1 1 1 1 1 1 1 1 ... $ Year : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ... $ ID : num 306 306 306 306 306 306 306 306 162 306 ... $ PS: num -0.32 -0.317 -0.317 -0.318 -0.317 ... $ TH: num -0.211 -0.211 -0.211 -0.213 -0.22 ... $ RS: num -0.337 -0.337 -0.337 -0.337 -0.337 ... $ HW: num -0.0258 -0.19 -0.19 -0.19 -0.4561 ... glmer(use ~ PS + TH + RS + HW + (1 + PS + TH + RS + HW |ID/Year), family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))
glmer 為我提供了對我有意義的固定效應的參數估計,當我對數據進行定性調查時,隨機斜率(我將其解釋為每種棲息地類型的選擇係數)也有意義。該模型的對數似然為 -3050.8。
然而,大多數動物生態學研究不使用 R,因為對於動物位置數據,空間自相關會使標準錯誤容易出現 I 類錯誤。雖然 R 使用基於模型的標準誤差,但首選經驗(也是 Huber-white 或三明治)標準誤差。
雖然 R 目前不提供此選項(據我所知 - 請糾正我,如果我錯了),SAS 提供 - 雖然我無法訪問 SAS,但一位同事同意讓我借用他的計算機以確定標準錯誤使用經驗方法時發生顯著變化。
首先,我們希望確保在使用基於模型的標準誤差時,SAS 會產生與 R 相似的估計值——以確保在兩個程序中以相同的方式指定模型。我不在乎它們是否完全相同 - 只是相似。我試過(SAS V 9.2):
proc glimmix data=dat method=laplace; class year id; model use = PS TH RS HW / dist=bin solution ddfm=betwithin; random intercept PS TH RS HW / subject = year(id) solution type=UN; run;title;
我還嘗試了其他各種形式,例如添加線條
random intercept / subject = year(id) solution type=UN; random intercept PS TH RS HW / subject = id solution type=UN;
我試過沒有指定
solution type = UN,
或註釋掉
ddfm=betwithin;
無論我們如何指定模型(並且我們嘗試了很多方法),我都無法讓 SAS 中的隨機斜率與 R 的輸出很相似——即使固定效果足夠相似。當我的意思是不同時,我的意思是甚至符號都不相同。SAS 中的 -2 對數似然是 71344.94。
我無法上傳完整的數據集;所以我做了一個玩具數據集,只有三個人的記錄。SAS 在幾分鐘內給我輸出;在 R 中需要一個多小時。奇怪的。有了這個玩具數據集,我現在得到了對固定效應的不同估計。
我的問題:誰能解釋為什麼 R 和 SAS 之間的隨機斜率估計可能如此不同?我可以在 R 或 SAS 中做些什麼來修改我的代碼,以便調用產生類似的結果?我寧願更改 SAS 中的代碼,因為我“相信”我的 R 估計更多。
我真的很關心這些差異,並想深入了解這個問題!
我的玩具數據集的輸出僅使用了 R 和 SAS 的完整數據集中 35 個個體中的三個,被包含為 jpeg。
編輯和更新:
正如@JakeWestfall 幫助發現的那樣,SAS 中的斜率不包括固定效應。當我添加固定效果時,結果如下 - 將 R 斜率與 SAS 斜率進行比較,以獲得程序之間的一個固定效果“PS”:(選擇係數 = 隨機斜率)。請注意 SAS 中增加的變化。
根據 Zhang et al 2011 的說法,我似乎不應該期望包之間的隨機斜率相似。在他們的論文On Fitting Generalized Linear Mixed-effects Models for Binary Responses using different Statistical Packages中,他們描述了:
抽象的:
廣義線性混合效應模型 (GLMM) 是將橫截面數據模型擴展到縱向設置的流行範式。當應用於建模二進制響應時,不同的軟件包,甚至一個包中的不同程序可能會給出完全不同的結果。在本報告中,我們描述了這些不同程序背後的統計方法,並討論了它們在應用於擬合相關二元響應時的優勢和劣勢。然後,我們通過將在一些流行的軟件包中實現的這些程序應用於模擬和真實研究數據來說明這些考慮因素。我們的模擬結果表明,所考慮的大多數程序都缺乏可靠性,這對在實踐中應用此類流行的軟件包具有重要意義。
我希望@BenBolker 和團隊將我的問題視為對 R 的投票,以將經驗標準誤差和 Gauss-Hermite 正交能力用於具有多個隨機斜率項的模型與 glmer 相結合,因為我更喜歡 R 接口並且希望能夠應用該程序中的一些進一步分析。令人高興的是,即使 R 和 SAS 沒有可比的隨機斜率值,但總體趨勢是相同的。感謝大家的投入,我非常感謝您為此付出的時間和考慮!