將隨機效應添加到 GAM 的兩種方法會產生非常不同的結果。為什麼會這樣,應該使用哪一個?
mgcv 文檔的一個特定部分給出了將隨機效應合併到廣義加法模型中的多種方法。兩種方法是 1) 使用
bs="re"
in在類標籤中添加平滑項gam
;2)結合現有的功能,使用功能gamm
,其中包括類似的設施。然而,在模擬數據上,兩者給出了完全不同的模型擬合。為什麼會這樣,應該使用哪一個?lme``gam
x <- rnorm(1000) ID <- rep(1:200,each=5) y <- x for(i in 1:200) y[which(ID==i)] <- y[which(ID==i)] + rnorm(1) y <- y + rnorm(1000) ID <- as.factor(ID) m1 <- gam(y ~ x + s(ID,bs="re")) m2 <- gamm(y ~ x, random=list(ID=~1) ) mean( (fitted(m1)-fitted(m2$gam))^2 )
我懷疑差異在於您獲得的擬合值。如果你看一下我所說的模型擬合,係數估計,方差項,模型是相同的。比較
summary(m2$lme)
和。summary(m1)
_gam.vcomp(m1)
> summary(m1) Family: gaussian Link function: identity Formula: y ~ x + s(ID, bs = "re") Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.05234 0.07932 0.66 0.51 x 1.01375 0.03535 28.68 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(ID) 167.1 199 5.243 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.674 Deviance explained = 72.9% GCV = 1.2133 Scale est. = 1.0082 n = 1000 > summary(m2$lme) Linear mixed-effects model fit by maximum likelihood Data: strip.offset(mf) AIC BIC logLik 3218.329 3237.96 -1605.165 Random effects: Formula: ~1 | ID (Intercept) Residual StdDev: 1.025306 1.003452 Fixed effects: y.0 ~ X - 1 Value Std.Error DF t-value p-value X(Intercept) 0.0523358 0.07922717 799 0.660578 0.5091 Xx 1.0137531 0.03535887 799 28.670404 0.0000 Correlation: X(Int) Xx 0.014 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.80375873 -0.67702485 0.04245145 0.64026891 2.59257295 Number of Observations: 1000 Number of Groups: 200
我們看到估計為是 1.01375 和兩種型號均約為 0.05。還要注意受試者/組間方差的估計值(作為標準偏差),在 的輸出中給出為 1.025
summary(m2$lme)
。相同的信息可以使用 來從gam
模型中計算出來gam.vcomp()
,它給出m1
> gam.vcomp(m1) s(ID) 1.027795
這已經足夠讓我們不用擔心了。
因此,這些
fitted
方法必須返回不同的擬合值;如果我們從 生成擬合值m2$lme
,那麼我們得到的值與由 生成的值相同fitted(m1)
:> mean((fitted(m1)-fitted(m2$lme))^2) [1] 2.966927e-07
這是出於所有意圖和目的0。
fitted.lme``ID
被記錄以從總體水平(平均超過)和針對特定主題的組件返回貢獻。這就是fitted.gam
要做的,m1
因為它將隨機效果表示為樣條“固定”效果。在gamm
模型的情況下,fitted.gam
返回模型“固定”效應部分的擬合值,這將解釋差異。(我寫“固定”是因為使用樣條曲線,“固定”和“隨機”效果會變得有點模糊。)在他的書中,Simon Wood 在使用
gamm()
. 他談到使用resid()
和fitted()
在模型的$lme
組件上gamm
分別排除和包括隨機效應。他說這適用於此處使用的特定示例中的模型診斷。您需要哪一個取決於您的具體使用/研究問題的背景。
如果您所需要的只是像這樣的簡單隨機效應,並且您熟悉 GAM 和mgcv,那麼使用隨機效應樣條基礎可能會更簡單,
gam()
而不必處理作為 GAMM 的混合的奇怪輸出模型通過gamm()
. 如上所示,這兩個模型實際上是等效的,您報告的差異僅取決於擬合值是否包含或排除受試者(或 ID)特定效果。