R

指定混合效應模型中不同組的相關結構 (lme4/nlme)

  • February 17, 2017

我試圖在 R 中的線性混合效應模型中解釋空間自相關,並及時重複測量。在4 年的時間裡,每150 種不同的產品中BodyMass收集一次。時間自相關應該可以忽略不計,因為體重測量是從死去的動物身上進行的。Year``Sites

由於 中的模型correlation=不接受表單的相關結構,因此我嘗試從較舊的包中取而代之。lmer()``lme4``lme()``nlme

我的基線模型如下所示:

library(nlme)
M1 <- lme(BodyMass~Year, data=df, random=~1|Site)

現在我想考慮空間自相關,考慮分組的空間位置(由於 4 年的周期,每個坐標組合重複 4 次)。但是,此語法不起作用:

M2 <- lme(BodyMass~Year, data=df, random=~1|Site, correlation=corAR1(form=~x+y|Year))

它給出了以下錯誤:

Error in lme.formula(BodyMass ~ 1, data = df, random = ~1 | Site,  :
incompatible formulas for groups in 'random' and 'correlation'

根據 Ben Bolker在此處的回答,該錯誤的原因是

lme()堅持隨機效應和相關性的分組變量是相同的。

然而,在我的模型中,使用相同的分組變量是沒有意義的,因為年份正在對相關結構進行分組(通過重複坐標),而年份不適合作為隨機因素,因為它們僅包含 4 個級別(如這裡討論)。

正如這裡所建議的那樣,使用gls()代替對我來說不是真正的選擇,因為我需要考慮隨機因素並且喜歡使用 R 平方,這只能在普通最小二乘 (OLS) 方法中使用。Site

(1) 你知道在這種情況下如何解釋空間自相關,保持我的隨機因子以及 OLS 框架能夠計算 R 平方值嗎?

非常感謝您的幫助!

編輯: AdamO回答lme4nlme鏈接)之間的差異時寫道:

“三巨頭”相關結構:獨立、可交換和 AR-1 相關結構很容易被這兩個包處理。

(2) 如何將這些“三大”相關結構包含在 中lme4

我認為如果實際上可以在其中處理空間自相關,則可以nlme通過使用來規避分組變量問題。lme4

**編輯2:**我發現如下更改相關結構的分組變量可以使模型正常工作而不會出現錯誤:

M3 <- lme(BodyMass~Year, data=df, random=~1|Site, correlation=corAR1(form=~x+y|Site/Year))

(3) 相關結構中分組變量的這個表達式是否有意義?

我不完全確定它是否確實如此,因為在一個隨機因素中,|Site/Year這意味著該年份嵌套在站點內,但我認為年份和站點寧願交叉,因為每個站點每年都有代表。

另一種非常靈活的方法是貝葉斯。您可以使用 JAGS 在 R 中實現它(除了一個包之外,您還必須通過一些步驟來下載)。如果你這樣做,你可以指定任何你想要的相關結構。

要以這種方式構建它,您可以 1) 將空間相關的結果視為多元正態模型的一部分(現在 y 有 2 個維度,即結果和空間)。或2) 為具有自己相關結構的模型添加另一個空間隨機分量。

例如,您可以修改 (2) 的代碼並在 R 中構建一個隨機截距和斜率模型。每個人都由 i 索引(例如 anx 1 向量),您還可以提供他們的空間索引的另一個 nx1 向量(稱為 site_indicator) . 您還需要提供站點總數。

model_RIAS_MVN<-"
model{
#Likelihood
for(i in 1:N_tot) {# all obs
# outcome is normally distributed
bodymass ~ dnorm(mu, sigmainverse)

# outcome model
mu <- b[1] + RI[subj_num[i], 1] + b[2]*Age[i] + RI[subj_num[i], 2]*Age[i] + RandomSpace[site_indicator[i]]
}

# Prior for random intercepts and slopes 
# this allows them to be correlated 
for (j in 1:N_people) {
RI[j, 1:2] ~ dmnorm(meanU[1:2], G.inv[1:2, 1:2])
}

# CHANGE HERE FOR NEW SITE CORRELATION #
# change number_sites in meanspace and G.invspace to actual number bc it # could throw error
for (j in 1:number_sites) {
RandomSpace[j] ~ dmnorm(meanspace[1:number_sites], G.invspace[1:number_sites, 1:number_sites])
}

for(i in 1:2){
meanU[i] <- 0 # zero mean for random components
}

G.inv[1:2, 1:2] ~ dwish(G0[1:2, 1:2], Gdf)
G[1:2, 1:2] <- inverse(G.inv[1:2, 1:2])

# whatever structure you want for correlation
G.invspace[1:number_sites, 1:number_sites] ~ dwish(G0inv[1:number_sites, 1:number_sites], Gdf_space)
Gspace[1:number_sites, 1:number_sites] <- inverse(G.invspace[1:number_sites, 1:number_sites])

sigmainverse ~ dgamma(1,1)

# informative priors for fixed effects
b[1] ~ dnorm(20, 0.25) 
b[2] ~ dnorm(1, 4) 

# uncomment for uninformative priors
# b[1] ~ dnorm(0, 0.01) 
# b[2] ~ dnorm(0, 0.01) 
}

這只是工作代碼,需要對其進行調整,但希望您了解靈活性以及如何以這種方式指定相關結構。

引用自:https://stats.stackexchange.com/questions/262573

comments powered by Disqus