Random-Effects-Model

將隨機效應添加到 GAM 的兩種方法會產生非常不同的結果。為什麼會這樣,應該使用哪一個?

  • February 22, 2016

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)特定效果。

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

comments powered by Disqus