R jagam 另一個 Gibbs Additive Modeller：JAGS 對 mgcv 的支持。

說明

auto-generate 模型規範代碼和關聯數據的設施，可在 JAGS(或 BUGS)中使用 GAM 進行模擬。這對於推斷具有最好在 JAGS 中編碼的複雜隨機效應結構的模型非常有用。這是對標準 GAM 進行推斷的一種非常低效的方法。這個想法是，jagam 為模型的平滑部分生成模板 JAGS 代碼和關聯數據。然後直接編輯該模板以包含其他隨機分量。在使用所得模型進行仿真之後，提供了用於使用模型平滑組件進行繪圖和預測的工具。

用法

jagam(formula,family=gaussian,data=list(),file,weights=NULL,na.action,
offset=NULL,knots=NULL,sp=NULL,drop.unused.levels=TRUE,
control=gam.control(),centred=TRUE,sp.prior = "gamma",diagonalize=FALSE)

sim2jam(sam,pregam,edf.type=2,burnin=0)

參數

formula	GAM 公式(請參閱 formula.gam 和 gam.models )。這與 GLM 的公式完全相同，隻是可以將平滑項 s 、 te 、 ti 和 t2 添加到右側以指定線性預測器依賴於預測器的平滑函數(或這些的線性泛函)。
family	這是一個族對象，指定要使用的分布和鏈接函數。有關更多詳細信息，請參閱glm 和family。目前僅支持高斯、泊鬆、二項式和伽瑪族，但用戶可以輕鬆修改 JAGS 代碼中的假設分布。
data	包含模型響應變量和公式所需的協變量的 DataFrame 或列表。默認情況下，變量取自 environment(formula) ：通常是調用 jagam 的環境。
file	應寫入 JAGS 模型規範代碼的文件的名稱。設置和查詢當前工作目錄請參見setwd。
weights	數據的先驗權重。
na.action	一個函數，指示當數據包含“NA”時應該發生什麽。默認值由 ‘options’ 的“na.action”設置設置，如果未設置，則為“na.fail”。 “factory-fresh” 默認為“na.omit”。
offset	可用於提供模型偏移以用於擬合。請注意，與 formula 中包含的偏移量不同，預測時此偏移量將始終被完全忽略：這符合 lm 和 glm 的行為。
control	用於替換 gam.control 返回的默認值的擬合控製參數列表。任何未提供的控製參數都保留其默認值。對 jagam 影響不大。
knots	這是一個可選列表，包含用戶指定的用於基礎構造的結值。對於大多數底座，用戶隻需提供要使用的結，該結必須與提供的 k 值匹配(請注意，結的數量並不總是隻是 k )。請參閱 tprs 了解 "tp"/"ts" 情況下發生的情況。不同的項可以使用不同數量的結，除非它們共享協變量。
sp	此處可以提供平滑參數向量。必須按照平滑項在模型公式中出現的順序提供平滑參數(不要忘記零空間懲罰)。負元素表示應該估計參數，因此固定參數和估計參數的混合是可能的。如果平滑共享平滑參數，則 length(sp) 必須對應於底層平滑參數的數量。
drop.unused.levels	默認情況下，未使用的級別會在擬合之前從因子中刪除。對於某些涉及因子變量的平滑，您可能需要將其關閉。僅當您知道自己在做什麽時才這樣做。
centred	是否應該像 GAMS 中那樣對平滑應用居中約束？僅當您確切知道自己在做什麽時才將其設置為FALSE。如果FALSE 每個平滑都有一個(通常是全局的)截距。
sp.prior	平滑參數先於"gamma" 或"log.uniform"？請檢查 JAGS 代碼中的默認參數是否適合您的模型。
diagonalize	是否應該重新參數化平滑以具有 i.i.d.高斯先驗(如果可能)？對於高斯數據，這允許使用高效的共軛采樣器，並且如果加載了 JAGS "glm" 模塊，它也可以與 GLM 很好地配合使用，但否則最好按塊更新平滑器，而不是這樣做。
sam	jags 樣本對象，至少包含字段b(係數)和rho(對數平滑參數)。還可能包含字段mu，其中包含受監控的預期響應。
pregam	標準mgcv GAM 設置數據，在jagam 返回列表中返回。
edf.type	由於 EDF 不是唯一定義的，並且可能會受到添加到 JAGS 模型文件中的隨機結構的影響，因此提供了 3 個選項。查看具體信息。
burnin	從模擬鏈中丟棄的燒錄量。限製為鏈長的 0.9。

細節

使用平滑係數的多元正態先驗，可以輕鬆地將平滑合並到 JAGS 模型中。平滑參數和平滑懲罰矩陣直接指定先驗多元正態精度矩陣。通常，平滑懲罰並不對應於滿秩精度矩陣，這意味著不適當的先驗不適合吉布斯采樣。為了糾正這個問題，Marra 和 Wood (2011) 中建議的零空間懲罰被添加到通常的懲罰中。

在加性建模環境中，通常將平滑居中，以避免與每個平滑項的截距(除了全局截距之外)相關的可識別性問題。在使用 JAGS 進行吉布斯采樣的情況下，技術上可以省略這種居中，因為無論如何我們都會強製先驗的適當性，而這種適當性意味著形式模型的可識別性。然而，在大多數情況下，這種形式上的可識別性是相當人為的，並不意味著具有統計意義的可識別性。相反，它隻會大量誇大置信區間，因為多個截距項無法從數據中識別出來，而隻能從先驗中識別出來。默認情況下，jagam 對所有平滑施加標準 GAM 可識別性約束。 centred 參數確實允許您關閉此函數，但不建議這樣做。如果您確實設置了centred=FALSE，那麽鏈收斂和混合檢查應該特別嚴格。

模型設置的最後一個技術問題是係數和平滑參數的初始條件的設置。采取的方法是采用 mgcv 在其他地方使用的默認初始平滑參數值，並對這些平滑參數執行單個 PIRLS 擬合步驟，以獲得平滑係數的起始值。在完全共軛更新的設置中，係數的初始值並不重要，並且在不提供它們的情況下也可以獲得良好的結果。但在通常的設置中，至少某些更新需要切片采樣，那麽有時在沒有初始值的情況下會獲得非常差的結果，因為采樣器根本無法找到後驗模式的區域。

sim2jam函數取部分gam目的 (pregam) 從jagam以及標準的模擬輸出rjags形式並創建一個簡化版本gam對象，適合繪製和預測模型的平滑分量。sim2gam計算每個平滑的有效自由度，但應該注意的是，在具有複雜隨機效應結構的模型的背景下，有多種可能性可以執行此操作。最簡單的方法(edf.type=0) 是計算如果平滑是未加權高斯加性模型的一部分，則該平滑將具有的自由度。如果模型已被修改，使得響應分布和/或鏈接不是指定的那些，則可以選擇使用此選項jagam。第二個選項是(edf.type=1) 使用 edf 計算得出gam如果它產生了這些估計 - 在 JAGS 模型修改都是關於修改隨機效應結構的情況下，這相當於簡單地將所有隨機效應設置為零以進行有效自由度計算。默認選項(edf.type=2) 是將 EDF 基於樣本協方差矩陣，Vp，模型係數。如果模擬輸出(sim)包括一個mu字段，那麽這將用於生成權重矩陣W在XWX = t(X)%*%W%*%X，其中 EDF 計算自rowSums(Vp*XWX)*scale.如果mu未提供，則根據模型矩陣進行估計X和模擬係數的平均值，但結果W可能不嚴格兼容Vp在這種情況下矩陣。在擬合模型的結構與生成的模板的回歸模型有很大不同的情況下jagam那麽默認選項可能沒有意義，實際上最好使用選項 0。

值

對於 jagam 包含三個項目的列表

pregam	標準mgcv GAM 設置數據。
jags.data	要提供給 JAGS 的參數列表，其中包含模型規範中引用的信息。
jags.ini	平滑係數和平滑參數的初始化數據。

對於 sim2jam 類的對象 "jam" ： mgcv gamObject 的部分版本，適合繪圖和預測。

警告

對於標準 GAM 來說，吉布采樣是一種非常慢的推理方法。僅當需要複雜的隨機效應結構高於直接 GAMM 方法所能達到的水平時，這才可能是值得的。

檢查參數的先驗參數是否適合您的目的。

例子

## the following illustrates a typical workflow. To run the 
## 'Not run' code you need rjags (and JAGS) to be installed.
require(mgcv)
  
set.seed(2) ## simulate some data... 
n <- 400
dat <- gamSim(1,n=n,dist="normal",scale=2)
## regular gam fit for comparison...
b0 <- gam(y~s(x0)+s(x1) + s(x2)+s(x3),data=dat,method="REML")

## Set directory and file name for file containing jags code.
## In real use you would *never* use tempdir() for this. It is
## only done here to keep CRAN happy, and avoid any chance of
## an accidental overwrite. Instead you would use
## setwd() to set an appropriate working directory in which
## to write the file, and just set the file name to what you
## want to call it (e.g. "test.jags" here). 

jags.file <- paste(tempdir(),"/test.jags",sep="") 

## Set up JAGS code and data. In this one might want to diagonalize
## to use conjugate samplers. Usually call 'setwd' first, to set
## directory in which model file ("test.jags") will be written.

jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file=jags.file,
            sp.prior="gamma",diagonalize=TRUE)

## In normal use the model in "test.jags" would now be edited to add 
## the non-standard stochastic elements that require use of JAGS....

## Not run: 
require(rjags)
load.module("glm") ## improved samplers for GLMs often worth loading
jm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
list.samplers(jm)
sam <- jags.samples(jm,c("b","rho","scale"),n.iter=10000,thin=10)
jam <- sim2jam(sam,jd$pregam)
plot(jam,pages=1)
jam
pd <- data.frame(x0=c(.5,.6),x1=c(.4,.2),x2=c(.8,.4),x3=c(.1,.1))
fv <- predict(jam,newdata=pd)
## and some minimal checking...
require(coda)
effectiveSize(as.mcmc.list(sam$b))

## End(Not run)

## a gamma example...
set.seed(1); n <- 400
dat <- gamSim(1,n=n,dist="normal",scale=2)
scale <- .5; Ey <- exp(dat$f/2)
dat$y <- rgamma(n,shape=1/scale,scale=Ey*scale)
jd <- jagam(y~s(x0)+te(x1,x2)+s(x3),data=dat,family=Gamma(link=log),
            file=jags.file,sp.prior="log.uniform")

## In normal use the model in "test.jags" would now be edited to add 
## the non-standard stochastic elements that require use of JAGS....

## Not run: 
require(rjags)
## following sets random seed, but note that under JAGS 3.4 many
## models are still not fully repeatable (JAGS 4 should fix this)
jd$jags.ini$.RNG.name <- "base::Mersenne-Twister" ## setting RNG
jd$jags.ini$.RNG.seed <- 6 ## how to set RNG seed
jm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
list.samplers(jm)
sam <- jags.samples(jm,c("b","rho","scale","mu"),n.iter=10000,thin=10)
jam <- sim2jam(sam,jd$pregam)
plot(jam,pages=1)
jam
pd <- data.frame(x0=c(.5,.6),x1=c(.4,.2),x2=c(.8,.4),x3=c(.1,.1))
fv <- predict(jam,newdata=pd)

## End(Not run)

作者

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

參考

Wood, S.N. (2016) Just Another Gibbs Additive Modeller: Interfacing JAGS and mgcv. Journal of Statistical Software 75(7):1-15 doi:10.18637/jss.v075.i07)

Marra, G. and S.N. Wood (2011) Practical variable selection for generalized additive models. Computational Statistics & Data Analysis 55(7): 2372-2387

Here is a key early reference to smoothing using BUGS (although the approach and smooths used are a bit different to jagam)

Crainiceanu, C. M. D Ruppert, & M.P. Wand (2005) Bayesian Analysis for Penalized Spline Regression Using WinBUGS Journal of Statistical Software 14.

也可以看看

gam , gamm , bam