4: The R-INLA package-----2019 Geospatial Health Data-Modeling and Visualization with R-INLA and Shi

4.2 .The inla() function

inla()函数是用来调整模型的。inla()的主要参数如下。

• formula：指定线性预测器的公式对象。

• data：包含数据的数据框。如果我们希望预测某些观测值的响应变量，我们需要将这些观测值的响应变量指定为NA。

• family：字符串或字符串的向量，表示似然族，如高斯、泊松或二项式。默认情况下，族是高斯族。输入names(inla.models()$likelihood)可以看到可能的替代列表，输入inla.doc(“familyname”)可以看到单个族的细节。

• control.compute：包含几个计算变量规格的列表，如dic是一个布尔变量，表示是否应该计算模型的DIC。

• control.predictor：包含几个预测变量的列表，如link是模型的链接函数，compute是一个布尔变量，表示是否应该计算线性预测的边际密度。

4.3 先验参数规格化

输入names(inla.models()$prior)可以看到R-INLA中可用的先验名称，输入inla.models()$prior可以看到包含每个先验选项的列表。关于特定先验的文档可以用inla.doc("priorname")查看。

默认情况下，模型的截距被分配了一个高斯先验，其平均值和精度都等于0。其余的固定影响被分配了高斯先验，其平均值等于0，精度等于0.001。这些值可以通过inla.set.control.fixed.default()[c("mean.intercept", "prec.intercept", "mean", "prec")]查看。这些先验值可以在inla()的control.fixed参数中改变，指定一个包含高斯分布的平均值和精度的列表。

具体来说，这个列表包含mean.intercept和prec.intercept，代表截距的先验平均值和精度。而mean和prec则代表了除截距以外的所有先验均值和精度。

prior.fixed < - list(mean.intercept = <>,  
	prec.intercept = <>, 
	mean = <>, prec = <>) 
res <- inla(formula,
	data = d,
	control.fixed = prior.fixed 
)

超参数θ的先验值被分配在f()的参数hyper中。

formula <- y ~ 1 + f(<>, model = <>, hyper = prior.f)

在inla()的parameter control.family 中分配似然参数的先验值。

res <- inla(formula, 
	data = d, 
	control.fixed = prior.fixed, 
	control.family = list(..., hyper = prior.l) 
)

hyper接受一个命名的列表，其名称等同于每个超参数，其值等同于一个标明先验参数的列表。具体来说，该列表包含以下值。

• initial：超参数的初始值（好的初始值可以使推理过程更快）。
• prior：先验分布的名称（例如，“iid”、“bym2”）。
• param：先验分布的参数值的向量。
• fixed：表示超参数是否为固定值的布尔变量。

prior.prec <- list(initial = <>, prior = <>, 
	param = <>, fixed = <>) 
prior <- list(prec = prior.prec)

先验需要在超参数的内部尺度中设置。例如，iid模型定义了一个精度为τ的独立高斯分布的随机变量向量。我们可以通过输入inla.doc(“iid”)来检查这个模型的文档，看到精度τ在对数尺度上被表示为log(τ)。因此，先验需要在对数精度log(τ)上进行定义。

4.4 Example

这里我们展示了一个示例，该示例演示了如何指定和适合一个模型，以及如何使用实际数据集和R-INLA检查结果。具体来说，我们对12家医院手术后的死亡率进行了建模。这项分析的目的是利用外科手术死亡率来评估每家医院的表现，并确定每家医院的表现是否异常出色或异常糟糕。

4.4.1 Data

我们使用的数据包括手术的数量和死亡的数字在12家医院进行心脏手术的婴儿。Surg是一个有三列的数据框架，即医院表示医院，n表示每家医院在一年内进行的手术数量，r表示每家医院在手术后30天内的死亡人数。

Surg

4.4.2 Model

我们指定一个模型来获得每个医院的死亡率。我们假设每家医院的死亡人数为二项式可能性，Y i，死亡率为pi

Y i ∼ Binomial(n i , p i ), i = 1, . . . , 12.

我们还假设各个医院的死亡率在某种程度上是相似的，并为真实死亡率pi指定了一个随机效应模型

logit(p i ) = α + u i , u i ∼ N(0, σ 2 ).

缺省情况下，为α指定一个非信息先验，即种群对数死亡率α ∼ N(0, 1/τ), τ = 0.
在R-INLA中，随机效应u i的精度的默认先验值为1/σ 2 ~ Gamma(1,5 × 10−5)。我们可以通过在标准差σ之前设置一个惩罚复杂度(PC)来改变这个先验。例如，我们可以指定σ大于1的概率小等于0.01:P(σ > 1) = 0.01。在R-INLA中，此先验被指定为

prior.prec <- list(prec = list(prior = "pc.prec", 
	param = c(1, 0.01)))

将模型转化为R代码，用下面的公式:

formula <- r ~ f(hospital, model = "iid", 
	hyper = prior.prec)

然后，我们调用inla()指定公式、数据、族和试验次数。我们添加了控制。predictor = list(compute = TRUE)计算参数的后边缘，control.compute = list(dic = TRUE)表示应该计算dic。

res <- inla(formula,
	data = Surg, family = "binomial", Ntrials = n,
	control.predictor = list(compute = TRUE),
	control.compute = list(dic = TRUE)
)

4.4.3 Results

当执行inla()时，我们得到一个类inla的对象，该对象包含了拟合模型的信息，包括固定效应、随机效应、超参数、线性预测因子和拟合值的摘要和后缘密度。返回对象res的摘要可以通过summary(res)查看。

summary(res)

res$dic$dic

结果：
[1] 74.93
可以通过输入res$summary.fixed来获得固定影响的摘要。这将返回一个数据框架，包括后验的平均值、标准差、2.5、50和97.5百分位数以及模型。kld列代表对称的Kullback-Leibler差异（Kullback and Leibler, 1951），描述了每个后验的高斯和简化或完全拉普拉斯近似之间的差异。

res$summary.fixed

res$summary.random

res$summary.hyperpar

在执行inla()时，如果在control.predictor中设置compute = TRUE，返回的对象还包括以下对象。

• summary.linear.predictor：包含线性预测器的平均值、标准差和定量的数据框。

• summary.fitted.values：数据框，包含通过链接函数的逆向转换线性预测器而得到的fitted值的平均值、标准差和量值。

• marginals.linear.predor：包含线性预测因子的后验边际的列表。

• marginals.fitted.values：通过对线性预测器进行反向链接函数转换而得到的fitted值的后验边际的列表。

请注意，如果一个观察值是NA，那么使用的链接函数是相同的。如果我们希望summary.fit.values和marginals.fit.values包含转化后的规模的拟合值，我们需要在control.predor中设置适当的链接。另外，我们也可以使用inla.tmarginal()函数手动转换inla对象中的边际。

我们的例子中预测的死亡率可以通过res$summary.fitted.values获得。

res$summary.fitted.values

mean列显示了，医院2、8和11是死亡率后验平均值最高的医院。
列 0.025quant和0.975quant两栏包含了死亡率95%可信区间的下限和上限，并提供了不确定性的测量。

我们还可以通过输入res$marginals.fixed获得一个包含固定effects的后验边缘的列表，通过输入marginals.random和marginals.hyperpar分别获得包含随机effects和超参数的后验边缘的列表。边际是指包含两列矩阵的列表。第x列代表参数的值，第y列是密度。R-INLA集成了几个函数来操作后验边际。例如，inla.emarginal()和inla.qmarginal()分别计算后验边际的期望值和量值。inla.smarginal()可用于获得样条平滑，inla.tmarginal()可用于转换边际，inla.zmarginal()提供汇总统计。

在我们的例子中，固定effects的后验边际的第一个元素，res$marginals.fixed[[1]]，包含截距α的后验元素。我们可以应用inla.smarginal()来获得边际密度的样条平滑，然后用ggplot2软件包的ggplot()函数来绘制（图4.1）。

library(ggplot2) 
alpha <- res$marginals.fixed[[1]] ggplot(data.frame(inla.smarginal(alpha)), aes(x, y)) +
	geom_line() + theme_bw()

分位数和分布函数分别由inla.qmarginal()和inla.pmarginal()给出。我们可以得到α的分位数0.05，然后画出α低于这个分位数的概率如下:

quant <- inla.qmarginal(0.05, alpha) 
quant

[1] -2.782

inla.pmarginal(quant, alpha)

[1] 0.05

α的概率低于0.05分位数的图可以创建如下(图4.2):

ggplot(data.frame(inla.smarginal(alpha)), aes(x, y)) +
	geom_line() + 
	geom_area(data = subset(data.frame(inla.smarginal(alpha)),
		x < quant),fill = "black") + 
	theme_bw()

函数inla.dmarginal()计算特定值处的密度。例如，-2.5处的密度可以计算如下(图4.3):

inla.dmarginal(-2.5, alpha)

[1] 2.989

ggplot(data.frame(inla.smarginal(alpha)), aes(x, y)) +
	geom_line() +
	geom_vline(xintercept = -2.5, linetype = "dashed") +
	theme_bw()

如果我们想要得到边际的变换，我们可以用 inla.tmarginal()。例如，如果我们希望得到随机效应u i的方差，我们可以得到精度τ的边际，然后应用逆函数。

marg.variance <- inla.tmarginal(function(x) 1/x, res$marginals.hyperpar$"Precision for hospital")

随机效应u i的方差的后验图如图4.4所示。

ggplot(data.frame(inla.smarginal(marg.variance)), aes(x, y)) +
	geom_line() + 
	theme_bw()

现在，如果我们想要得到方差的平均后验，我们可以使用inla. margininal()。

m <- inla.emarginal(function(x) x, marg.variance) 
m

[1] 0.1465
我们也可以使用inla.zmarginal()来获取边际的汇总统计信息。

inla.zmarginal(marg.variance)

还有几个图和绘图函数未翻译。