R gam.mh 具有 gam 擬合的簡單後驗模擬

說明

GAM 係數可以直接從高斯近似到係數的後驗進行模擬，或者使用簡單的 Metropolis Hastings 采樣器。另請參閱ginla。

用法

gam.mh(b,ns=10000,burn=1000,t.df=40,rw.scale=.25,thin=1)

參數

b	來自 gam 的擬合模型對象。不支持bam 擬合。
ns	要生成的樣本數。
burn	要丟棄的任何初始燒錄期的長度(除了代碼之外)。
t.df	靜態多元 t 建議的自由度。對於較重的尾部提案，較低。
rw.scale	生成隨機遊走提案時縮放後驗協方差矩陣的因子。負或非有限跳過隨機遊走步驟。
thin	僅保留每個 thin 樣本。

細節

後驗模擬對於推斷模型係數的非線性函數特別有用。模擬後驗隨機抽取，計算每次抽取的函數，然後從該函數的後驗中進行抽取。在許多情況下，模型係數後驗的高斯近似是準確的，並且從中生成的樣本可以被視為來自係數後驗的樣本。請參閱下麵的示例代碼。這種方法在計算上非常高效。

在其他情況下，高斯近似可能會變得很差。一個典型的例子是在具有 log 或 logit 鏈接的空間模型中，當存在大麵積的觀測值僅包含零時。在這種情況下，線性預測器很難識別，高斯近似可能變得毫無用處(下麵提供了一個示例)。在這種情況下，使用 Metropolis Hastings 采樣器從後驗模擬有時會很有用。一種簡單的方法將基於後驗高斯近似的固定提案與基於近似後驗協方差矩陣的縮小版本的隨機遊走提案交替使用。 gam.mh 實現了這一點。固定提案通常促進快速混合，而隨機遊走組件確保鏈不會陷入固定高斯提案密度遠低於後驗密度的區域。

該函數報告兩種類型步驟的接受率。如果隨機遊走接受概率高於四分之一，那麽 rw.step 可能應該增加。同樣，如果接受率太低，則應降低接受率。隨機遊走步驟可以完全關閉(見上文)，但如果這樣做，檢查鏈是否有卡住的部分很重要。

值

包含矩陣 bs 中保留的模擬係數和接受概率的兩個條目的列表。

例子

library(mgcv)
set.seed(3);n <- 400

############################################
## First example: simulated Tweedie model...
############################################

dat <- gamSim(1,n=n,dist="poisson",scale=.2)
dat$y <- rTweedie(exp(dat$f),p=1.3,phi=.5) ## Tweedie response
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=tw(),
          data=dat,method="REML")

## simulate directly from Gaussian approximate posterior...
br <- rmvn(1000,coef(b),vcov(b))

## Alternatively use MH sampling...
br <- gam.mh(b,thin=2,ns=2000,rw.scale=.15)$bs
## If 'coda' installed, can check effective sample size
## require(coda);effectiveSize(as.mcmc(br))

## Now compare simulation results and Gaussian approximation for
## smooth term confidence intervals...
x <- seq(0,1,length=100)
pd <- data.frame(x0=x,x1=x,x2=x,x3=x)
X <- predict(b,newdata=pd,type="lpmatrix")
par(mfrow=c(2,2))
for(i in 1:4) {
  plot(b,select=i,scale=0,scheme=1)
  ii <- b$smooth[[i]]$first.para:b$smooth[[i]]$last.para
  ff <- X[,ii]%*%t(br[,ii]) ## posterior curve sample
  fq <- apply(ff,1,quantile,probs=c(.025,.16,.84,.975))
  lines(x,fq[1,],col=2,lty=2);lines(x,fq[4,],col=2,lty=2)
  lines(x,fq[2,],col=2);lines(x,fq[3,],col=2)
}

###############################################################
## Second example, where Gaussian approximation is a failure...
###############################################################

y <- c(rep(0, 89), 1, 0, 1, 0, 0, 1, rep(0, 13), 1, 0, 0, 1, 
       rep(0, 10), 1, 0, 0, 1, 1, 0, 1, rep(0,4), 1, rep(0,3),  
       1, rep(0, 3), 1, rep(0, 10), 1, rep(0, 4), 1, 0, 1, 0, 0, 
       rep(1, 4), 0, rep(1, 5), rep(0, 4), 1, 1, rep(0, 46))
set.seed(3);x <- sort(c(0:10*5,rnorm(length(y)-11)*20+100))
b <- gam(y ~ s(x, k = 15),method = 'REML', family = binomial)
br <- gam.mh(b,thin=2,ns=2000,rw.scale=.4)$bs
X <- model.matrix(b)
par(mfrow=c(1,1))
plot(x, y, col = rgb(0,0,0,0.25), ylim = c(0,1))
ff <- X%*%t(br) ## posterior curve sample
linv <- b$family$linkinv
## Get intervals for the curve on the response scale...
fq <- linv(apply(ff,1,quantile,probs=c(.025,.16,.5,.84,.975)))
lines(x,fq[1,],col=2,lty=2);lines(x,fq[5,],col=2,lty=2)
lines(x,fq[2,],col=2);lines(x,fq[4,],col=2)
lines(x,fq[3,],col=4)
## Compare to the Gaussian posterior approximation
fv <- predict(b,se=TRUE)
lines(x,linv(fv$fit))
lines(x,linv(fv$fit-2*fv$se.fit),lty=3)
lines(x,linv(fv$fit+2*fv$se.fit),lty=3)
## ... Notice the useless 95% CI (black dotted) based on the
## Gaussian approximation!

作者

Simon N. Wood simon.wood@r-project.org

參考

Wood, S.N. (2015) Core Statistics, Cambridge