审查数据模型#

from copy import copy

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import seaborn as sns

from numpy.random import default_rng

%config InlineBackend.figure_format = 'retina'
rng = default_rng(1234)
az.style.use("arviz-darkgrid")

这个关于贝叶斯生存分析的示例笔记本涉及到了审查数据的问题。审查是一种缺失数据问题，其中大于某个阈值的观测值被截断到该阈值，或小于某个阈值的观测值被截断到该阈值，或两者兼有。这些分别称为右审查、左审查和区间审查。在这个示例笔记本中，我们考虑区间审查。

审查数据出现在许多建模问题中。两个常见的例子是：

生存分析：在研究某种医疗治疗对生存时间的影响时，不可能将研究延长到所有受试者都死亡为止。在研究结束时，许多患者收集到的唯一数据是他们在治疗后存活了一段时间 \(T\)：实际上，他们的真实生存时间大于 \(T\)。
传感器饱和：传感器可能有一个有限的范围，上下限仅仅是传感器能够报告的最高和最低值。例如，许多水银温度计只能报告一个非常窄的温度范围。

这个示例笔记本介绍了在PyMC3中处理截断数据的两

一个插补的删失模型，它将删失数据表示为参数，并为所有删失值生成合理的值。由于这种插补，该模型能够生成合理的插补值集，这些值原本会被删失。每个删失元素引入一个随机变量。
一个未插补的删失模型，其中删失数据被积分出去，仅通过对数似然来考虑。该方法更适当地处理大量删失数据，并且收敛更快。

为了建立一个基线，我们与未审查数据的未审查模型进行比较。

# Produce normally distributed samples
size = 500
true_mu = 13.0
true_sigma = 5.0
samples = rng.normal(true_mu, true_sigma, size)

# Set censoring limits
low = 3.0
high = 16.0


def censor(x, low, high):
    x = copy(x)
    x[x <= low] = low
    x[x >= high] = high
    return x


# Censor samples
censored = censor(samples, low, high)

# Visualize uncensored and censored data
_, ax = plt.subplots(figsize=(10, 3))
edges = np.linspace(-5, 35, 30)
ax.hist(samples, bins=edges, density=True, histtype="stepfilled", alpha=0.2, label="Uncensored")
ax.hist(censored, bins=edges, density=True, histtype="stepfilled", alpha=0.2, label="Censored")
[ax.axvline(x=x, c="k", ls="--") for x in [low, high]]
ax.legend();

../../../_images/5d33add778ae781898173c0065c76cd136bc93722573ac843c9f89db0a03f7c2.png

未审查模型#

def uncensored_model(data):
    with pm.Model() as model:
        mu = pm.Normal("mu", mu=((high - low) / 2) + low, sigma=(high - low))
        sigma = pm.HalfNormal("sigma", sigma=(high - low) / 2.0)
        observed = pm.Normal("observed", mu=mu, sigma=sigma, observed=data)
    return model

我们应该预测，在未审查的数据上运行未审查的模型，我们将得到合理的均值和方差的估计。

uncensored_model_1 = uncensored_model(samples)
with uncensored_model_1:
    idata = pm.sample()

az.plot_posterior(idata, ref_val=[true_mu, true_sigma], round_to=3);

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu, sigma]

100.00% [8000/8000 00:03<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.

../../../_images/8691c1e31edbedbb5a0e5e4aa7f5364043f254dcb1aae85490ac4509244432f6.png

而这正是我们所发现的。

然而，问题在于在审查数据的情况下，我们无法访问真实值。如果我们对审查数据使用相同的未审查模型，我们预计参数估计会有偏差。如果我们计算均值和标准差的点估计，那么我们可以看到我们可能会低估这个特定数据集和审查边界的均值和标准差。

print(f"mean={np.mean(censored):.2f}; std={np.std(censored):.2f}")

mean=12.32; std=3.76

uncensored_model_2 = uncensored_model(censored)
with uncensored_model_2:
    idata = pm.sample()

az.plot_posterior(idata, ref_val=[true_mu, true_sigma], round_to=3);

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu, sigma]

100.00% [8000/8000 00:03<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.

../../../_images/0d343ef387a38347fd570208b0137c0401559c61706d5045ac8577a63c2fc4d6.png

上图证实了这一点。