空间数据的 Besag-York-Mollie 模型#

import os

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt

from scipy import sparse
from scipy.linalg import solve
from scipy.sparse.linalg import spsolve

注意

此 notebook 使用的库不是 PyMC 的依赖项，因此需要专门安装才能运行此 notebook。打开下面的下拉菜单以获取更多指导。

额外的依赖项安装说明

为了运行此 notebook（本地或在 binder 上），您不仅需要安装可用的 PyMC 以及所有可选依赖项，还需要安装一些额外的依赖项。 有关安装 PyMC 本身的建议，请参阅 安装

您可以使用您喜欢的包管理器安装这些依赖项，我们提供 pip 和 conda 命令作为示例如下。

$ pip install nutpie networkx

请注意，如果您想（或需要）从 notebook 内部而不是命令行安装软件包，您可以通过运行 pip 命令的变体来安装软件包

import sys

!{sys.executable} -m pip install nutpie networkx

您不应运行 !pip install，因为它可能会将软件包安装在不同的环境中，即使已安装也无法从 Jupyter notebook 中使用。

另一种选择是使用 conda

$ conda install nutpie networkx

当使用 conda 安装科学 python 软件包时，我们建议使用 conda forge

# these libraries are not dependencies of pymc

import networkx as nx
import nutpie

RANDOM_SEED = 8926
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

为什么要使用 Besag-York-Mollie 模型？#

此 notebook 解释了为什么以及如何在 PyMC 中部署 Besag-York-Mollie (BYM) 模型。BYM 模型是解决许多空间统计问题的有吸引力的方法。它很灵活 - 一旦添加了 BYM 组件，其余工作流程就像任何其他贝叶斯广义线性模型一样进行。您可以添加预测变量来估计因果效应。您可以交换链接函数和结果分布来处理不同的数据类型。您可以混合和匹配任何效果最佳的采样器。

BYM 在大型数据集上也具有良好的可扩展性。空间模型的常见问题是，随着数据集大小的增加，其计算成本迅速增长。PyMC 的 CAR 模型 就是这种情况。 使用 BYM 模型，计算成本的增长几乎是线性的。

BYM 模型适用于面状数据，例如相邻的州、县或人口普查区。对于涉及空间点或连续距离度量的问题，请考虑使用 高斯过程 代替。

由 ICAR 驱动#

驱动 BYM 模型的主要引擎是内在条件自回归先验 (ICAR)。ICAR 是一种特殊的多元正态分布，它假设相邻区域共同变化。

在谈论空间建模时，采用一些图论的词汇会有所帮助。图由节点和边组成。节点代表空间中的区域，边代表邻近关系。在这种类型的问题中，我们在共享边界的两个区域之间绘制一条边。

假设我们有一个如图所示的图，该图由邻接矩阵 W 构建。

# adjacency matrix

W = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 0], [0, 1, 0, 0]])

# plotting routine for graphs

G = nx.Graph(W)
colors = np.array([0.8, 1, 1.5, 2])
nx.draw_networkx(G, node_color=colors, vmin=0, vmax=2, node_size=500, cmap="plasma")
../_images/f9978a554480ea0aa2fb3ab3d68855b1bc7fd0b0956807394252b4c7d84a3b93.png

邻接矩阵编码了哪些节点连接到哪些其他节点。如果节点 i 和 j 连接，则在第 i 行、第 j 列的位置将为 1。否则，将为零。例如，节点 1 和 3 连接。矩阵中第 3 行、第 1 列的位置为 1。但是节点 3 未连接到节点 2，因此第 3 行第 2 列的位置为 0。我们当然必须记住 python 索引。第一行是 0，最后一行是 3。邻接矩阵也是对称的 - 如果加拿大与美国相邻，则美国与加拿大相邻。

ICAR 的密度函数接受邻接矩阵 W 和方差 \(\sigma\)。我们通常假设 \(\sigma = 1\) 并在其他方面处理方差，因此我现在将忽略第一个分数。

\[f(\phi| W,\sigma) = -\frac{1}{2\sigma^{2}} \sum_{i\sim j} (\phi_{i} - \phi_{j})^2\]

每个 \(\phi_{i}\) 都根据其每个邻居的平方距离进行惩罚。符号 \(i \sim j\) 表示对 \(\phi_{i}\) 的所有邻居求和。

因此，例如，想象一下颜色的强度表示每个节点处变量的值。节点 1 连接到所有人。节点 1 和 0 具有相当相似的颜色，因此惩罚会很小。但是节点 2 的颜色与节点 1 的颜色非常不同，因此它会受到很大的惩罚。如果我们将注意力转移到节点 3，它只有一个邻居，并且只会根据与节点 1 的距离受到一次大的惩罚。

通过这种方式，ICAR 编码了空间统计的核心假设 - 邻近区域应比遥远区域更相似。最可能的结果是每个节点都具有相同值的图。在这种情况下，邻居之间的平方距离始终为零。图中相邻区域之间发生突然变化越多，对数密度越低。

ICAR 还有一些其他特殊功能：它受到约束，因此所有 \(\phi\) 的总和为零。这也意味着 \(\phi\) 的均值为零。将 ICAR 值视为类似于 z 分数可能会有所帮助。它们表示以 0 为中心的相对偏差。ICAR 通常也仅用作较大模型的子组件。模型的其他部分通常调整尺度（使用方差参数）或位置（使用截距参数）。关于 ICAR 背后的数学及其与 CAR 关系的易于理解的讨论可以在 此处 或学术论文版本 [Morris et al., 2019] 中找到。

随机效应的灵活性#

统计建模的典型目标之一是将数据的方差划分为三个类别：由感兴趣的原因解释的方差、结构化方差和非结构化方差。在我们的例子中，ICAR 模型旨在捕获（空间）结构化方差。添加预测变量可以处理第一个类别。BYM 模型使用随机效应 \(\theta\) 来处理第三个类别。随机效应是长度为 n 的随机变量向量，其中 n 是区域的数量。它旨在捕获所有未由空间或因果效应解释的剩余方差。

构建包含结构化方差和非结构化方差的模型可能很棘手。朴素的方法经常遇到可识别性问题。原则上，每个组件都可以独立解释方差。因此，拟合算法可能无法在参数空间中的小邻域中稳定下来。

BYM 模型经过精心设计以解决可识别性问题。它使用混合分布，其中参数 \(\rho\) 控制结构化方差与非结构化方差的平衡。BYM 模型看起来像这样

\[\beta + \sigma(\sqrt{1 - \rho} \theta + \sqrt{\rho / \text{s}} \phi)\]

当 \(\rho\) 接近 1 时，大部分方差是空间结构化的。当 \(\rho\) 接近 0 时，大部分方差是非结构化的。

\(\sigma\) 是 \(\theta\) 和 \(\phi\) 共享的尺度参数。 \(\theta\) 和 \(\phi\) 都以零为中心，方差为 1。因此，它们都像 z 分数一样起作用。\(\sigma\) 可以拉伸或缩小效应的混合，使其适合实际数据。\(\beta\) 是一个共享截距，它重新居中混合以拟合数据。最后，\(\text{s}\) 是比例因子。它是从邻接矩阵计算的常数。它重新缩放 \(\phi\)，使其具有与 \(\theta\) 相同的预期方差。关于为什么这有效 的更详细讨论将在下面出现。

拟合此模型解决了将方差划分为结构化和非结构化组件的挑战。剩下的唯一挑战是确定预测变量，这是一个因案例而异的挑战。 Riebler et al. [2016] 提出了 BYM 模型的这种特定方法，并提供了更多关于为什么 BYM 模型的这种参数化既可解释又可识别的解释，而以前的 BYM 模型参数化通常不可解释和不可识别。

在纽约市行人事故数据集上演示 BYM 模型#

我们将在记录纽约市行人交通事故数量的数据集上演示 BYM 模型。数据组织成大约 2000 个普查区，提供了我们的空间结构。我们的目标是证明我们可以将方差划分为已解释的、空间结构化的和非结构化的组件。

设置数据#

空间数据以 边列表 的形式出现。除了邻接矩阵之外，边列表是在计算机上表示面状数据的另一种常用技术。边列表是一对列表，用于存储有关图中边的信息。假设 i 和 j 是两个节点的名称。如果节点 i 和节点 j 连接，则一个列表将包含 i，另一个列表将在同一行包含 j。例如，在下面的数据框中，节点 1 连接到节点 1452 以及节点 1721。

try:
    df_edges = pd.read_csv(os.path.join("..", "data", "nyc_edgelist.csv"))
except FileNotFoundError:
    df_edges = pd.read_csv(pm.get_data("nyc_edgelist.csv"))

df_edges
N N_edges node1 node2
0 1921 5461 1 1452
1 1921 5461 1 1721
2 1921 5461 2 3
3 1921 5461 2 4
4 1921 5461 2 5
... ... ... ... ...
5456 1921 5461 1918 1919
5457 1921 5461 1918 1921
5458 1921 5461 1919 1920
5459 1921 5461 1919 1921
5460 1921 5461 1920 1921

5461 行 × 4 列

但是，为了实际运行我们的模型，我们需要将边列表转换为邻接矩阵。下面的代码执行了该任务以及一些其他清理任务。

# convert edgelist to adjacency matrix

# extract and reformat the edgelist

nodes = np.stack((df_edges.node1.values, df_edges.node2.values))
nodes = nodes.T

# subtract one for python indexing

nodes = nodes - 1

# convert the number of nodes to a integer

N = int(df_edges.N.values[0])

# build a matrix of zeros to store adjacency
# it has size NxN where N is the number of
# areas in the dataset

adj = np.zeros((N, N))

# loop through the edgelist and assign 1
# to the location in the adjacency matrix
# to represent the edge
# this will only fill in the upper triangle
# of the matrix

for node in nodes:
    adj[tuple(node)] = 1

# add the transpose to make the adjacency
# matrix symmetrical

W_nyc = adj.T + adj

我们将计算比例因子。这将需要一个相当特殊的函数。对该函数的正确解释将使我们远离 NYC 案例研究，因此我将 比例因子的讨论留到以后。

def scaling_factor_sp(A):
    """Compute the scaling factor from an adjacency matrix.
    This function uses sparse matrix computations and is most
    efficient on sparse adjacency matrices. Used in the BYM2 model.
    The scaling factor is a measure of the variance in the number of
    edges across nodes in a connected graph.
    Only works for fully connected graphs. The argument for scaling
    factors is developed by Andrea Riebler, Sigrunn H. Sørbye,
    Daniel Simpson, Havard Rue in "An intuitive Bayesian spatial
    model for disease mapping that accounts for scaling"
    https://arxiv.org/abs/1601.01180"""

    # Computes the precision matrix in sparse format
    # from an adjacency matrix.

    num_neighbors = A.sum(axis=1)
    A = sparse.csc_matrix(A)
    D = sparse.diags(num_neighbors, format="csc")
    Q = D - A

    # add a small jitter along the diagonal

    Q_perturbed = Q + sparse.diags(np.ones(Q.shape[0])) * max(Q.diagonal()) * np.sqrt(
        np.finfo(np.float64).eps
    )

    # Compute a version of the pseudo-inverse

    n = Q_perturbed.shape[0]
    b = sparse.identity(n, format="csc")
    Sigma = spsolve(Q_perturbed, b)
    A = np.ones(n)
    W = Sigma @ A.T
    Q_inv = Sigma - np.outer(W * solve(A @ W, np.ones(1)), W.T)

    # Compute the geometric mean of the diagonal on a
    # precision matrix.

    return np.exp(np.sum(np.log(np.diag(Q_inv))) / n)


scaling_factor = scaling_factor_sp(W_nyc)
scaling_factor

0.7136574058611103

第一个 .csv 文件只有空间结构位。其余数据是单独提供的 - 在这里我们将提取事故数量 y 和普查区的人口规模 E。我们将使用人口规模作为偏移量 - 我们应该预期人口稠密的地区由于微不足道的原因会发生更多事故。更有趣的是与某个地区相关的超额风险。

最后，我们还将探索一个预测变量，即社会碎片化指数。该指数由空置住房单元数量、独居人口、租房者和前一年内搬家的人口等指标构成。这些社区往往融合度较低，社会支持系统较弱。社会流行病学界对生态变量如何渗透到公共卫生的各个方面感兴趣。因此，我们将看看社会碎片化是否可以解释交通事故的模式。该指标已标准化为均值为零，标准差为 1。

try:
    df_nyc = pd.read_csv(os.path.join("..", "data", "nyc_traffic.csv"))
except FileNotFoundError:
    df_nyc = pd.read_csv(pm.get_data("nyc_traffic.csv"))

y = df_nyc.events_2001.values
E = df_nyc.pop_2001.values
fragment_index = df_nyc.fragment_index.values

# Most census tracts have huge populations
# but a handful have 0. We round
# those up to 10 to avoid triggering an error
# with the log of 0.

E[E < 10] = 10
log_E = np.log(E)
area_idx = df_nyc["census_tract"].values
coords = {"area_idx": area_idx}

我们可以通过可视化邻接矩阵来了解空间结构。下图仅捕获了普查区的相对位置。它没有考虑绝对位置，因此看起来不像纽约市。此表示突出了城市如何由几个均匀连接的区域、一些具有大量连接的中心枢纽以及一些狭窄的走廊组成。

# build the positions of the nodes. We'll only
# generate the positions once so that we can
# compare visualizations from the data to
# the model predictions.
# I found that running spectral layout first
# and passing it to spring layout makes easy to read
# visualizations for large datasets.

G_nyc = nx.Graph(W_nyc)
pos = nx.spectral_layout(G_nyc)
pos = nx.spring_layout(G_nyc, pos=pos, seed=RANDOM_SEED)

# the spread of the data is pretty high. Most areas have 0 accidents.
# one area has 300. Color-gradient based visualization doesn't work
# well under those conditions. So for the purpose of the color
# we'll cap the accidents at 30 using vmax
#
# however, we'll also make the node size sensitive to the real
# number of accidents. So big yellow nodes have way more accidents
# than small yellow nodes.

plt.figure(figsize=(10, 8))
nx.draw_networkx(
    G_nyc,
    pos=pos,
    node_color=y,
    cmap="plasma",
    vmax=30,
    width=0.5,
    alpha=0.6,
    with_labels=False,
    node_size=20 + 3 * y,
)
../_images/9e90708d64b17b75736c36089b5e622a76dd48f5be82a2f4de7f457d964f18dc.png

该地图还显示，有许多热点是大多数事故发生的地方，它们在空间上聚集，即右下角、底部中心和中心左侧区域。这是一个很好的迹象，表明空间自相关模型是一个合适的选择。

我们还可以可视化社会碎片化的空间布局。您会注意到，有一个社会碎片化社区与交通事故地图重叠。下面的统计分析将帮助我们了解这种重叠的真实强度。

plt.figure(figsize=(10, 8))
nx.draw_networkx(
    G_nyc,
    pos=pos,
    node_color=fragment_index,
    cmap="plasma",
    width=0.5,
    alpha=0.6,
    with_labels=False,
    node_size=40 + 5 * fragment_index,
)
../_images/e358b2dae33b327a31a314b6d384dc026c26c2982474034ecb7082b7c27ab667.png

使用 PyMC 指定 BYM 模型#

BYM 的所有参数已在 第 1 节 中介绍。现在只需分配一些先验。 \(\theta\) 的先验很挑剔 - 我们需要分配均值为 0 和标准差为 1，以便我们可以将其解释为与 \(\phi\) 可比。否则，先验分布将提供纳入领域专业知识的机会。在这个问题中，我将选择一些弱信息先验。

最后，我们将使用泊松结果分布。交通事故的数量是计数结果，最大可能值非常大。为了确保我们的预测保持为正，我们将在将线性模型传递给泊松分布之前对其进行指数化。

with pm.Model(coords=coords) as BYM_model:
    # intercept
    beta0 = pm.Normal("beta0", 0, 1)

    # fragmentation effect
    beta1 = pm.Normal("beta1", 0, 1)

    # independent random effect
    theta = pm.Normal("theta", 0, 1, dims="area_idx")

    # spatially structured random effect
    phi = pm.ICAR("phi", W=W_nyc, dims="area_idx")

    # joint variance of random effects
    sigma = pm.HalfNormal("sigma", 1)

    # the mixing rate is rho
    rho = pm.Beta("rho", 0.5, 0.5)

    # the bym component - it mixes a spatial and a random effect
    mixture = pm.Deterministic(
        "mixture", pt.sqrt(1 - rho) * theta + pt.sqrt(rho / scaling_factor) * phi, dims="area_idx"
    )

    # exponential link function to ensure
    # predictions are positive
    mu = pm.Deterministic(
        "mu", pt.exp(log_E + beta0 + beta1 * fragment_index + sigma * mixture), dims="area_idx"
    )

    y_i = pm.Poisson("y_i", mu, observed=y)

模型采样#

# if you haven't installed nutpie, it's okay to to just delete
# 'nuts_sampler="nutpie"'. The default sampler took roughly 12 minutes on
# my machine.

with BYM_model:
    idata = pm.sample(1000, nuts_sampler="nutpie", random_seed=rng)

采样器进度

链总数：4

活动链：0

已完成链：4

采样 31 秒

预计完成时间：now

进度 抽取 发散 步长 梯度/抽取
2000 0 0.13 255
2000 0 0.13 127
2000 0 0.13 63
2000 0 0.13 63

我们可以通过多种方式评估采样器。首先，看起来我们所有的链都已收敛。所有参数的 rhat 值都非常接近于 1。

rhat = az.summary(idata, kind="diagnostics").r_hat.values
sum(rhat > 1.03)

0

其次，所有主要参数的迹图看起来都是平稳且混合良好的。它们还显示 rho 的均值在 0.50 左右，表明数据中可能存在空间效应。

az.plot_trace(idata, var_names=["beta0", "beta1", "sigma", "rho"])
plt.tight_layout();
../_images/eeb50ab49cbea694089ed58ee0567f63d75614c4ab8a5cf0bd3028c99c12bcc0.png

我们的迹图还表明，社会碎片化对交通事故有很小的影响，后验质量的大部分在 0.06 到 0.12 之间。

后验预测检查#

所有这些工作的回报是，我们现在可以可视化将方差分解为解释性部分、空间部分和非结构化部分意味着什么。使此生动形象的一种方法是单独检查模型的每个组件。我们将看到如果空间效应是方差的唯一来源，模型认为纽约市应该是什么样子，然后我们将转向解释性效应，最后是随机效应。

在第一种情况下，我们将仅可视化来自模型空间组件的预测。换句话说，我们假设 \(\rho = 1\)，并且我们忽略 \(\theta\) 和社会碎片化。

然后，我们将我们的预测覆盖到相同的 我们之前构建的邻接图 上。

# draw posterio

with pm.do(BYM_model, {"rho": 1.0, "beta1": 0}):
    y_predict = pm.sample_posterior_predictive(
        idata, var_names=["mu", "mixture"], predictions=True, extend_inferencedata=False
    )

y_spatial_pred = y_predict.predictions.mu.mean(dim=["chain", "draw"]).values

plt.figure(figsize=(10, 8))
nx.draw_networkx(
    G_nyc,
    pos=pos,
    node_color=y_spatial_pred,
    cmap="plasma",
    vmax=30,
    width=0.5,
    alpha=0.6,
    with_labels=False,
    node_size=20 + 3 * y_spatial_pred,
)

Sampling: []





../_images/693fcaaf5427da5ad28684383f941c1c51d9ee1a9b3392b7fe5e653f3585f0a9.png





结果图片称为空间平滑。附近区域往往非常相似，从而导致风险的独特邻域。在深紫色区域，方差很小，预测的事故数量也很少，接近于零。


空间平滑对于预测特别有用。想象一下，有一个低事故区被高事故邻域包围。假设您想预测未来哪里会发生高事故数量，以便您可以针对这些区域进行干预。仅关注过去事故数量高的区域环可能是一个错误。中间的焦点低事故区过去可能只是运气好。将来，该区域可能会更像其邻居，而不是其过去。空间平滑依赖于与部分合并背后的相同原理 - 我们可以通过合并来自附近区域的信息来纠正异常，从而学到更多。


我们可以注意到，有三个风险邻域，以大型黄色集群表示，它们被很好地捕获。这表明，对交通事故的许多解释与未识别但空间结构化的原因有关。相比之下，社会碎片化指数仅解释了地图底部中心的一个风险邻域（以及其他地方的一些小口袋）。






with pm.do(
    BYM_model,
    {
        "sigma": 0.0,
    },
):
    y_predict = pm.sample_posterior_predictive(
        idata, var_names=["mu", "mixture"], predictions=True, extend_inferencedata=False
    )

y_frag_pred = y_predict.predictions.mu.mean(dim=["chain", "draw"]).values

plt.figure(figsize=(10, 8))
nx.draw_networkx(
    G_nyc,
    pos=pos,
    node_color=y_frag_pred,
    cmap="plasma",
    vmax=30,
    width=0.5,
    alpha=0.6,
    with_labels=False,
    node_size=20 + 3 * y_frag_pred,
)









Sampling: []









../_images/96923dd7ebd9eb8c250591104f96f7e87b1c81befbc8f0329d85b95de0e482dc.png





最后，我们可以通过假设 \(\rho = 0\) 来查看非结构化方差。如果我们的模型成功地划分了方差，则非结构化方差中不应剩下太多空间集群。相反，方差应该分散在整个地图上。






with pm.do(BYM_model, {"rho": 0.0, "beta1": 0}):
    y_predict = pm.sample_posterior_predictive(
        idata, var_names=["mu", "mixture"], predictions=True, extend_inferencedata=False
    )

y_unspatial_pred = y_predict.predictions.mu.mean(dim=["chain", "draw"]).values

plt.figure(figsize=(10, 8))
nx.draw_networkx(
    G_nyc,
    pos=pos,
    node_color=y_unspatial_pred,
    cmap="plasma",
    vmax=30,
    width=0.5,
    alpha=0.6,
    with_labels=False,
    node_size=20 + 3 * y_unspatial_pred,
)









Sampling: []









../_images/96b46861f838cf4e0ce37c9548ed6c36c5e5dac093ebbc8a35d4e97a42a542e3.png











比例因子实际上是做什么的？#


关于 BYM 模型的讨论通常省略了对比例因子的过多细节讨论。这有充分的理由。如果您的主要兴趣是流行病学，您真的不需要了解它。用户可以让它只是一个黑匣子。比例因子的计算还涉及线性代数中一些相当晦涩难懂的思想。我不会在这里介绍计算，但我会尝试为它在 BYM 模型中扮演的角色提供一些直觉。


看看这两个图。






W1 = np.array([[0, 1, 1, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]])

G = nx.Graph(W1)
nx.draw(G)









c:\Users\dsaun\Anaconda3\envs\spatial_pymc_env\Lib\site-packages\IPython\core\events.py:89: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
  func(*args, **kwargs)
c:\Users\dsaun\Anaconda3\envs\spatial_pymc_env\Lib\site-packages\IPython\core\pylabtools.py:152: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
  fig.canvas.print_figure(bytes_io, **kw)




../_images/e584fb36b50577a5f27a53569643363b73ac48a9b0075f51fa5db86100be6497.png









W2 = np.array([[0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 0, 1], [1, 1, 1, 0]])

G = nx.Graph(W2)
nx.draw(G)








../_images/23d437b93ec7f67bff2673405b2918337e4aa12fd429b1caae135a482c2c1fb8.png





如果节点之间存在很强的空间协方差，我们应该预期第一个图允许比第二个图更大的方差。在第二个图中，每个节点都对每个其他节点施加影响。因此，结果应该相对均匀。


比例因子是衡量特定邻接矩阵隐含多少方差的指标。如果我们计算上面两个矩阵的比例因子，它会证实我们的直觉。第一个图允许比第二个图更大的方差。






scaling_factor_sp(W1), scaling_factor_sp(W2)









(0.31249999534338707, 0.18749999767169354)









第二个例子可以真正强调这一点。这些是两个优先连接图 - 少数节点有很多边，而大多数节点只有很少的边。唯一的区别是最小边数。在第一个图中，每个节点至少获得两条边。在第二个图中，每个节点至少有一条边。






G = nx.barabasi_albert_graph(50, 2)
nx.draw(G)
W_sparse = nx.adjacency_matrix(G, dtype="int")
W = W_sparse.toarray()
print("scaling factor: " + str(scaling_factor_sp(W)))









scaling factor: 0.41031348581585425




../_images/f6a41f34f809c85c93fc4f4f453e3433da652e4b3e7d4e2df40065e4002e654a.png









G = nx.barabasi_albert_graph(50, 1)
nx.draw(G)
W_sparse = nx.adjacency_matrix(G, dtype="int")
W = W_sparse.toarray()
print("scaling factor: " + str(scaling_factor_sp(W)))









scaling factor: 1.7744552306768082




../_images/c6da00dc95ad939d2b5838802c4dc2fceffd74eebf1c71130b9f2fc96600a65b.png





第二个图具有更高的比例因子，因为它不太均匀连接。有更多机会形成具有独特特征的小口袋。在第一个图中，连接的规律性缓和了这种机会。同样，比例因子证实了我们的直觉。


这在很大程度上澄清了比例因子衡量的标准。但是为什么我们需要使用它呢？让我们回顾一下 BYM 组件的数学描述


\[\beta + \sigma(\sqrt{1 - \rho} \theta + \sqrt{\rho / \text{s}} \phi)\]


BYM 模型的目的是我们将两种不同类型的随机效应混合在一起，然后 \(\sigma\) 提供混合物的总体方差。这意味着我们需要非常小心每个随机效应的个体方差 - 它们都需要大约等于 1。很容易确保 \(\theta \approx 1\) 的方差。我们可以将其指定为先验的一部分。获得 \(\phi \approx 1\) 的方差更难，因为方差来自数据（空间结构），而不是来自先验。


比例因子是确保 \(\phi\) 的方差大致等于 1 的技巧。当空间结构隐含的方差非常小时，例如，小于 1，将 \(\rho\) 除以比例因子将得到一些大于 1 的数字。换句话说，我们扩展 \(\phi\) 的方差，直到它等于 1。现在所有其他参数都将正常运行。\(\rho\) 表示两个相似事物之间的混合，\(\sigma\) 表示来自随机效应的联合方差。


理解比例因子目的的最后一种方法是想象一下如果我们不包含它会发生什么。假设该图暗示了非常大的方差，就像上面的第一个优先连接图一样。在这种情况下，混合参数 \(\rho\) 可能会更多地拉入 \(\phi\)，因为数据具有很大的方差，并且模型正在搜索它可以找到的方差来解释它。但这使得结果的解释具有挑战性。 \(\rho\) 趋向于 \(\phi\) 是因为实际上存在很强的空间结构吗？还是因为它比 \(\theta\) 具有更高的方差？除非我们重新缩放 \(\phi\)，否则我们无法分辨。






作者#







参考文献#







[1]

Andrea Riebler, Sigrunn H Sørbye, Daniel Simpson, 和 Håvard Rue. 一种直观的贝叶斯空间模型，用于疾病绘图，可解释缩放。《统计医学研究方法》，25(4):1145–1165, 2016。





[2]

Mitzi Morris, Katherine Wheeler-Martin, Dan Simpson, Stephen J. Mooney, Andrew Gelman, 和 Charles DiMaggio. 贝叶斯分层空间模型：在 Stan 中实现 Besag York Mollié 模型。《空间和时空流行病学》，31:100301, 2019。URL：https://www.sciencedirect.com/science/article/pii/S1877584518301175, doi:https://doi.org/10.1016/j.sste.2019.100301。












水印#






%load_ext watermark
%watermark -n -u -v -iv -w -p xarray









Last updated: Thu Aug 01 2024

Python implementation: CPython
Python version       : 3.11.0
IPython version      : 8.8.0

xarray: 2024.3.0

matplotlib: 3.6.2
pytensor  : 2.14.2
networkx  : 3.0
nutpie    : 0.13.0
numpy     : 1.24.3
pymc      : 5.16.0
pandas    : 2.2.2
scipy     : 1.10.0
arviz     : 0.17.1

Watermark: 2.3.1











许可声明#


此示例画廊中的所有 notebook 均根据 MIT 许可证 提供，该许可证允许修改和再分发以用于任何用途，前提是保留版权和许可声明。






引用 PyMC 示例#


要引用此 notebook，请使用 Zenodo 为 pymc-examples 存储库提供的 DOI。




重要提示


许多 notebook 都改编自其他来源：博客、书籍……在这种情况下，您也应该引用原始来源。


另请记住引用您的代码使用的相关库。




这是一个 bibtex 中的引用模板


@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}





一旦呈现，它可能看起来像