可下载资源

本文分析的数据、代码、报告分享至会员群

在此，固定效应 (z\_i\beta) 通过与区域 (i) 对应的截距和 (p) 个协变量的向量来表示，即 (z\_i=(1, z_{i1}, \cdots, z_{ip}))，系数向量为 (\beta = (\beta\_0, \cdots, \beta\_p)’)。

该模型包含一个空间随机效应 (u\_i)，它用于解释结果之间的空间依赖性，表明彼此接近的区域可能具有相似的值，还有一个无结构可交换分量 (v\_i) 用于对不相关噪声进行建模。空间随机效应 (u_i) 可使用内在条件自回归模型（CAR）进行建模，该模型会根据特定的邻域结构对数据进行平滑处理。具体而言：

（一）贝叶斯推断与INLA

贝叶斯分层模型可通过多种方法进行拟合，如集成嵌套拉普拉斯近似（INLA）（Rue, Martino, and Chopin 2009）和马尔可夫链蒙特卡罗（MCMC）方法（Gelman et al. 2000）。INLA是一种在潜在高斯模型中进行近似贝叶斯推断的计算方法，它涵盖了广泛的模型，如广义线性混合模型、空间和时空模型等。INLA通过结合解析近似和数值积分来获得参数的近似后验分布，与MCMC方法相比速度非常快。

​

然后，我们使用 `[nb2INLA()](https://r-spatial.github.io/spdep/reference/nb2INLA.html "nb2INLA()")` 函数将 nb 列表转换为一个名为 map.adj 的文件，该文件包含了R-INLA所需的邻域矩阵表示形式。map.adj 文件保存在当前工作目录中，可通过 `[getwd()](https://rdrr.io/r/base/getwd.html "getwd()")` 函数获取。然后，我们使用R-INLA的 `[inla.read.graph()](https://rdrr.io/pkg/INLA/man/read.graph.html "inla.read.graph()")` 函数读取 map.adj 文件，并将其存储在对象 g 中，稍后我们将使用该对象通过R-INLA来指定空间模型。

（四）模型公式与 inla() 调用

我们通过包含响应变量、~ 符号以及固定效应和随机效应来指定模型公式。默认情况下有一个截距，所以我们不需要在公式中包含它。在公式中，随机效应通过 `[f()](https://rdrr.io/pkg/INLA/man/f.html "f()")` 函数来指定。该函数的第一个参数是一个索引向量，用于指定适用于每个观测值的随机效应元素，第二个参数是模型名称。对于空间随机效应 (u_i)，我们使用 model = "besag" 并将邻域矩阵设为 g。选项 scale.model = TRUE 用于使具有不同CAR先验的模型的精度参数具有可比性（Freni-Sterrantino, Ventrucci, and Rue 2018）。对于无结构效应 (v_i)，我们选择 model = "iid"。随机效应的索引向量分别为 re_u 和 re_v，它们是为 (u_i) 和 (v_i) 创建的，这些向量的值为 (1, \cdots, n)，其中 (n) 是区域数量。

map$re_u <- 1:nrow(map)
map$re_v <- 1:nrow(map)
formula <- vble ~ CRIM + RM +  f(re_u, model = "besag", graph = g, scale.model = TRUE) +  f(re_v, model = "iid")

需要注意的是，在R-INLA中，BYM模型也可通过 model = "bym" 来指定，这将同时包含空间和无结构组件。或者，我们也可以使用BYM2模型（Simpson et al. 2017），它是BYM模型的一种新参数化形式，使用了缩放后的空间组件 (u^) 和无结构组件 (v^)： ​

从上述结果中我们可以观察到，截距项  的估计值为 ，其  可信区间为 (, )。而对于协变量犯罪率（CRIM），其系数估计值  为 ，对应的  可信区间为 (, )，这表明犯罪率与房价之间存在显著的负相关关系。房间数（RM）的系数  为 ， 可信区间为 (, )，这意味着房间数与房价之间存在显著的正相关关系。由此可见，犯罪率和房间数这两个因素在解释房价的空间分布模式方面都起着重要作用。

​

 2. 响应变量后验分布概要

我们还可以通过输入 res$summary.fitted.values 来获取每个区域响应变量  的后验分布概要信息。其中，“mean” 列表示后验均值，“0.025quant” 列和 “0.975quant” 列分别表示  可信区间的下限和上限，它们代表了所获得估计值的不确定性程度。其输出结果如下：

summary(res$summary.fitted.values)

​ ​

随后，我们使用 mapview 包来创建这些变量的地图。在创建地图过程中，我们为这三张地图指定了一个通用的图例，并使用一个弹出式表格，其中包含区域名称、房价对数、协变量以及后验均值和  可信区间等信息。同时，我们利用 leafsync 包中的 sync() 函数来绘制同步地图，具体步骤如下：

​

# 通用图例
at <- seq(min(c(map$PM, map$LL, map$UL)),
          max(c(map$PM, map$LL, map$UL)),
          length.out = 8)

然后，创建弹出式表格：

# 弹出式表格
popuptable <- leafpop::popupTable(dplyr::mutate_if(map,
                                                    is.numeric, round, digits = 2),
                                  zcol = c("TOWN", "vble", "CRIM", "RM", "PM", "LL", "UL"),
                                  row.numbers = FALSE, feature.id = FALSE)

​

4. 获取原始房价规模的估计值

为了得到原始房价规模的估计值，我们需要对房价对数的估计值进行转换。首先，使用 inla.tmarginal() 函数来获取价格的边际分布，通过 exp(log(price)) 的方式进行转换。然后，再利用 inla.zmarginal() 函数来获取这些边际分布的概要信息。最后，创建变量 PMoriginal、LLoriginal 和 ULoriginal，它们分别对应原始房价后验分布的后验均值以及  可信区间的下限和上限，具体代码如下：

# 对第一个区域的边际分布进行转换（示例）
# inla.tmarginal(function(x) exp(x),
#                res$marginals.fitted.values[[1]])

# 对所有边际分布进行转换
marginals <- lapply(res$marginals.fitted.values,
                    FUN = function(marg){inla.tmarginal(function(x) exp(x), marg)})

# 获取边际分布的概要信息
marginals_summaries <- lapply(marginals,
                              FUN = function(marg){inla.zmarginal(marg)})

# 后验均值和95%可信区间
map$PMoriginal <- sapply(marginals_summaries, '[[', "mean") 
map$LLoriginal <- sapply(marginals_summaries, '[[', "quant0.025")
map$ULoriginal <- sapply(marginals_summaries, '[[', "quant0.975")

同样地，我们可以基于这些变量创建地图来展示原始房价估计值及其  可信区间，以便更好地理解波士顿房价的空间分布模式以及估计值的不确定性。具体创建地图的步骤与上述类似，只是这里使用的是 PMoriginal、LLoriginal 和 ULoriginal 变量，生成的地图可参考下图

​

通过上述一系列的分析和可视化操作，我们能够较为全面地了解基于贝叶斯分层模型对波士顿房价数据进行建模分析的结果，以及各因素对房价空间分布的影响和估计值的不确定性情况。

​  ​