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