R ginla GAM 集成嵌套拉普拉斯逼近牛顿增强

说明

使用 Wood (2019) 中说明的 INLA 变体，将集成嵌套拉普拉斯逼近 (INLA, Rue et al. 2009) 应用于可通过 gam 或 bam 估计的模型。为每个系数、选定系数或系数向量的线性变换生成边际后验密度。

用法

ginla(G,A=NULL,nk=16,nb=100,J=1,interactive=FALSE,int=0,approx=0)

参数

G	预拟合 gam 对象，由 gam(...,fit=FALSE) 或 bam(...,discrete=TRUE,fit=FALSE) 生成。
A	感兴趣的系数的变换矩阵，或者感兴趣的参数的索引数组。如果NULL，则为所有系数生成分布。
nk	用于评估其对数边后验密度的每个系数的值的数量。然后对这些点进行样条插值。
nb	评估系数后验密度以作为网格函数返回的点数。
J	在对数行列式近似步骤中要采取多少个行列式更新步骤。不建议增加此值。
interactive	如果这是 >0 或 TRUE，则每个近似后验都以红色绘制，覆盖在简单的高斯近似后验上。如果2则等待用户在每个绘图之间按回车键。对于判断使用 INLA 方法是否有所收获很有用。
int	0 跳过积分并仅使用后验模态平滑参数。 >0 用于使用 Rue 等人提出的 CCD 方法进行积分。 (2009)。
approx	0 表示完全近似； 1 更新Hessian，但使用近似模式； 2 作为 1 并假设常数 Hessian。查看具体信息。

细节

让\beta,\theta和y分别表示模型系数、超参数/平滑参数和响应数据。原则上，INLA 采用拉普拉斯近似\pi(\beta_i|\theta,y)和\pi(\theta|y)然后得到边后验分布\pi(\beta_i|y)通过将近似值积分\pi(\beta_i|\theta,y)\pi(\theta|y)相对于\theta(也可以生成超参数的边际)。在实践中，拉普拉斯近似为\pi(\beta_i|\theta,y)计算每个的成本太高\beta_i并且本身必须是近似的。为此，必须计算两个量：后验模式\beta^*|\beta_i以及联合对数密度的 Hessian 矩阵的行列式\log \pi(\beta,\theta,y)w.r.t.\beta在模式下。鲁等人。 (2009)最初通过对后验的简单高斯近似所隐含的条件模式来近似后验条件模式\pi(\beta|y)。然后，他们将 Hessian 矩阵的对数行列式近似为以下函数：\beta_i使用一阶泰勒展开式，对于他们使用的稀疏模型表示来说，计算起来很便宜，但在使用密集的低秩基础展开式时则不然gam。他们还提供了一种更昂贵的替代近似，基于计算仅与以下元素相关的对数行列式：\beta具有足够高的相关性\beta_i根据简单的高斯后验近似：效率似乎再次取决于稀疏性。 Wood (2018) 建议精确计算所需的后验模式，并将对数行列式近似基于无条件模型下 Hessian 的 BFGS 更新。无论有或没有稀疏性，后者都是有效的，而前者是“免费”的改进。这两个步骤都很有效，因为获得 Cholesky 因子的成本很低H[-i,-i]从那H- 看choldrop。这是该例程所采用的方法。

approx 参数允许两个进一步的近似来加速计算。对于approx==1，不使用精确的后验条件模式，而是使用简单高斯后验近似隐含的条件模式。对于approx==2，模态使用相同的近似值，并且假定 Hessian 矩阵为常数。后者非常快，因为不需要对数联合密度梯度评估。

请注意，对于许多模型，INLA 估计值非常接近通常的后验高斯近似，interactive 参数对于研究此问题非常有用。

bam 模型仅支持 disrete=TRUE 选项。 discrete=FALSE 方法效率太低。不支持 AR1 模型(相关参数将被忽略)。

值

包含元素 beta 和 density 的列表，这两个元素都是矩阵。每一行与一个感兴趣的系数(或线性系数组合)相关。两个矩阵都有nb 列。如果int!=0，则另一个元素reml给出CCD积分中使用的积分权重，首先给出中心点权重。

警告

此例程仍处于实验阶段，因此细节可能会发生变化。此外，目前并非所有步骤都具有最佳效率。

该例程是为相对专业的用户编写的。

ginla 并非旨在处理排名不足的模型。

例子

  require(mgcv); require(MASS)

  ## example using a scale location model for the motorcycle data. A simple plotting
  ## routine is produced first...

  plot.inla <- function(x,inla,k=1,levels=c(.025,.1,.5,.9,.975),
               lcol = c(2,4,4,4,2),lwd = c(1,1,2,1,1),lty=c(1,1,1,1,1),
	       xlab="x",ylab="y",cex.lab=1.5) {
    ## a simple effect plotter, when distributions of function values of
    ## 1D smooths have been computed
    require(splines)
    p <- length(x) 
    betaq <- matrix(0,length(levels),p) ## storage for beta quantiles 
    for (i in 1:p) { ## work through x and betas
      j <- i + k - 1
      p <- cumsum(inla$density[j,])*(inla$beta[j,2]-inla$beta[j,1])
      ## getting quantiles of function values...
      betaq[,i] <- approx(p,y=inla$beta[j,],levels)$y
    }
    xg <- seq(min(x),max(x),length=200)
    ylim <- range(betaq)
    ylim <- 1.1*(ylim-mean(ylim))+mean(ylim)
    for (j in 1:length(levels)) { ## plot the quantiles
      din <- interpSpline(x,betaq[j,])
      if (j==1) {
        plot(xg,predict(din,xg)$y,ylim=ylim,type="l",col=lcol[j],
             xlab=xlab,ylab=ylab,lwd=lwd[j],cex.lab=1.5,lty=lty[j])
      } else lines(xg,predict(din,xg)$y,col=lcol[j],lwd=lwd[j],lty=lty[j])
    }
  } ## plot.inla

  ## set up the model with a `gam' call...

  G <- gam(list(accel~s(times,k=20,bs="ad"),~s(times)),
            data=mcycle,family=gaulss(),fit=FALSE)
  b <- gam(G=G,method="REML") ## regular GAM fit for comparison

  ## Now use ginla to get posteriors of estimated effect values
  ## at evenly spaced times. Create A matrix for this...
  
  rat <- range(mcycle$times)
  pd0 <- data.frame(times=seq(rat[1],rat[2],length=20))
  X0 <- predict(b,newdata=pd0,type="lpmatrix")
  X0[,21:30] <- 0
  pd1 <- data.frame(times=seq(rat[1],rat[2],length=10))
  X1 <- predict(b,newdata=pd1,type="lpmatrix")
  X1[,1:20] <- 0
  A <- rbind(X0,X1) ## A maps coefs to required function values

  ## call ginla. Set int to 1 for integrated version.
  ## Set interactive = 1 or 2 to plot marginal posterior distributions
  ## (red) and simple Gaussian approximation (black).
 
  inla <- ginla(G,A,int=0)

  par(mfrow=c(1,2),mar=c(5,5,1,1))
  fv <- predict(b,se=TRUE) ## usual Gaussian approximation, for comparison

  ## plot inla mean smooth effect...
  plot.inla(pd0$times,inla,k=1,xlab="time",ylab=expression(f[1](time))) 

  ## overlay simple Gaussian equivalent (in grey) ...
  points(mcycle$times,mcycle$accel,col="grey")
  lines(mcycle$times,fv$fit[,1],col="grey",lwd=2)
  lines(mcycle$times,fv$fit[,1]+2*fv$se.fit[,1],lty=2,col="grey",lwd=2)
  lines(mcycle$times,fv$fit[,1]-2*fv$se.fit[,1],lty=2,col="grey",lwd=2)

  ## same for log sd smooth...
  plot.inla(pd1$times,inla,k=21,xlab="time",ylab=expression(f[2](time)))
  lines(mcycle$times,fv$fit[,2],col="grey",lwd=2)
  lines(mcycle$times,fv$fit[,2]+2*fv$se.fit[,2],col="grey",lty=2,lwd=2)
  lines(mcycle$times,fv$fit[,2]-2*fv$se.fit[,2],col="grey",lty=2,lwd=2)

  ## ... notice some real differences for the log sd smooth, especially
  ## at the lower and upper ends of the time interval.

作者

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

参考

Rue, H, Martino, S. & Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B. 71: 319-392.

Wood (2019) Simplified Integrated Laplace Approximation. In press Biometrika.