从具有 2 个预测变量的分层模型中绘制多个散点图

如何解决从具有 2 个预测变量的分层模型中绘制多个散点图

我是 python 和 pymc3 的新手。如果这个问题相当愚蠢，请原谅我。我有一个名为 toy_data.csv 的数据集：

batch_no	batch_id	点	pred1	pred2
12-150	1	1	70.26	10.766
12-150	1	2	65.72	10.512
12-150	1	6	55.44	7.78
12-150	1	4	63.554	12.552
12-150	1	5	61.66	10.95
12-150	1	3	65.44	10.808
12-190	2	3	55.68	6.36
12-190	2	4	67.36	7.38
12-190	2	1	78	9.32
12-190	2	2	70.12	8.554
12-190	3	5	60.68	9.568
12-190	3	3	68.694	9.634
12-190	3	2	56.118	11.474
12-190	3	4	58.54	9.31
12-190	3	6	65.08	10.604
13-999	4	6	56.73	9.754
13-999	4	3	55.492	9.4
13-999	4	2	51.556	9.802
13-999	4	1	53.748	8.85
13-999	4	5	59.054	9.204
13-999	4	4	49.77	9.468
13-999	4	7	58.266	9.954
13-999	4	8	57.78	9.38
14-140	5	2	69.68	12.086
14-140	5	1	68.5	10.438
14-140	5	6	71.7	11.624
14-140	5	3	63.68	11.058
14-140	5	4	62.02	10.498
14-140	5	5	61.94	10.95
14-140	5	9	57.22	10.164
14-140	5	7	54.44	9.798
14-140	5	8	60.82	10.784
14-290	6	1	56.2	9.566
14-290	6	1	50.06	9.23
14-290	6	2	50.76	10.646
14-290	6	2	46.98	8.38
14-700	7	2	92.8	11.532
14-700	7	1	81	9.522
14-700	7	3	75.62	10.004
14-700	8	2	71.44	10.076
14-889	8	1	55.2	8.1
14-889	8	3	71.18	9.412
14-889	8	4	53.24	7.138
14-900	9	3	50.08	7.744
14-900	9	1	47.138	8.294
14-900	9	4	68.56	11.362
14-900	9	1	69.34	11.584
14-900	9	2	63.14	10.77

我想从分层模型跟踪中获取相关批次和预测器来绘制它们并拟合 abline 线。最终输出是复制下面的图表，但不是灰度范围，我想显示每个批次的 abline。

例如：

Abline = y_best_fit_i=(mu_a+a_i)+[(mu_b+b_i)*pred_i] 使用来自相应 x 变量（pred1 或 pred2）的均值和标准差对所有参数进行回缩。仅考虑 pred1，后缩放公式为 mu_a = (mu_a*pred1.std()) + pred1.mean()​​​；这将应用于 mu_b 和所有 a_i 和 b_i 值。但是，按照 pred1 均值和 std dev 的后缩放仅适用于来自 pred1 的跟踪，而不适用于 pred2（即它将适用于 mu_a、mu_b 和 a[0,0] 和 b[0,0]，但不适用于a[0,1,0] 或 b[0,0]）。我想只适用于 mu_a,mu_b,a[i,0] 和 b[i,0] 在我应该是 batch_no 的地方（见下表），如果可能的话，让它自动化。下表是跟踪的摘要。

如果我能找到一种方法从跟踪本身获取数字以绘制最终图形，那就太好了。非常感谢任何帮助。

我用来建模的代码如下：

import pandas as pd
import numpy as np
import pymc3 as pm
from sklearn import preprocessing

seed = 68492
np.random.seed(seed)

data = pd.read_csv("toy_data.csv")    

# assigning variables
batch_no = data.iloc[:,0]
batch_id = (data.iloc[:,1])-1
point = data.iloc[:,2]
pred = data.iloc[:,3:5]
pred1 = data.iloc[:,3]
pred2 = data.iloc[:,4]


# creating indicators
n_pred = pred.shape[1]
n_batch = len(np.unique(batch_id))

# standarising variable
std_pred = preprocessing.scale(pred)


# Hierarchical model 
with pm.Model() as hierarchical_model1:

    # Hyperpriors
    mu_a = pm.normal('mu_a',mu=0,sigma=100)
    sigma_a = pm.Halfnormal('sigma_a',sigma=100)
    mu_b = pm.normal('mu_b',sigma=100)
    sigma_b = pm.Halfnormal('sigma_b',sigma=100)   
    

    # Intercept for each batch,distributed around group mean mu_a
    a = pm.normal('a',mu=mu_a,sd=sigma_a,shape=(n_batch,n_pred,1))
    # Slope for each batcht,distributed around group mean mu_b
    b = pm.normal('b',mu=mu_b,sd=sigma_b,1))
    
    # Model error
    eps = pm.HalfCauchy('eps',beta=5)
    
        
    # Expected value
    omega = (pm.math.dot(std_pred,b))[batch_id,1]
    tpoint_est = a[batch_id,1] + omega

         
    # Data likelihood
    y_like = pm.normal('y_like',mu=tpoint_est,sigma=eps,observed=point)
    

with hierarchical_model1:
    hier_trace = pm.sample(draws = 4000,tune = 1000,chains = 2)
    hier_out = pm.summary(hier_trace)

解决方法

这是我所拥有的，它应该可以解决批处理/预测循环问题并作为 你真正想做的事情的骨架。我无法遵循缩放和重新缩放。我也无法理解长度为 1 的 3dr 维度 但我一直把它放在那里。

请注意，我使用的是建立在 xarray 之上的 ArviZ（PyMC3 的依赖项） 使用基于标签的索引和自动广播。而且我实际上正在使用添加到 ArviZ（开发版）和 PyMC3 3.11.1 的最新功能。您可以使用 pip install git+git://github.com/arviz-devs/arviz.git 安装 ArviZ 开发版。

我首先对模型进行了一些更改，以定义具有命名维度的变量的形状：

coords = {
    "batch": np.arange(n_batch),"pred_num": ["pred1","pred2"],"one": [1]
}
with pm.Model(coords=coords) as hierarchical_model1:
    ...

    # Intercept for each batch,distributed around group mean mu_a
    a = pm.Normal('a',mu=mu_a,sd=sigma_a,dims=("batch","pred_num","one"))
    # Slope for each batcht,distributed around group mean mu_b
    b = pm.Normal('b',mu=mu_b,sd=sigma_b,"one"))
    
    ...

其余代码完全相同，除了在我添加 kwarg return_inferencedata=True 的地方进行采样时。

采样后，我从 pred 数组创建了一个 xarray 数据集，以便利用我上面提到的 xarray 提供的自动广播。因为我已经设置了 return_inferencedata=True 我已经有了所有的变量 在 xarray 对象的模型中。要创建一个 xarray 变量，我们必须提供 值、维度名称和坐标变量：

import xarray as xr
pred_ds = xr.DataArray(
    pred,dims=("x","pred_num"),coords={"pred_num": coords["pred_num"]}
).to_dataset(name="pred")
post = hier_trace.posterior
y_ds = (post["mu_a"]+post["a"])+(post["mu_b"]+post["b"])*pred_ds
y_ds # in jupyter you'll see a nice html interactive summary
# here is the plain text output
# <xarray.Dataset>
# Dimensions:   (batch: 9,chain: 4,draw: 1000,one: 1,pred_num: 2,x: 48)
# Coordinates:
#   * pred_num  (pred_num) <U5 'pred1' 'pred2'
#   * chain     (chain) int64 0 1 2 3
#   * draw      (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
#   * batch     (batch) int64 0 1 2 3 4 5 6 7 8
#   * one       (one) int64 1
# Dimensions without coordinates: x
# Data variables:
#     pred      (chain,draw,batch,pred_num,one,x) float64 1.938e+03 ... 67.98

然后我们可以使用 ArviZ 遍历所需的变量并生成一个子图 对于每个批次和 pred 的组合。虽然看起来很奇怪，因为 根本没有重新缩放。

import matplotlib.pyplot as plt
fig,axes = plt.subplots(n_pred,n_batch,figsize=(11,5),sharex=True,sharey="row")
x = np.arange(len(batch_id))

第一步是创建我们将用于绘图的 suplots，以及一个虚拟变量来表示 x。

from arviz.labels import BaseLabeller
labeller = BaseLabeller()
iterator = az.sel_utils.xarray_var_iter(
    y_ds.squeeze().transpose("chain","draw","x","batch"),combined=True,skip_dims={"x",}
)

然后我们创建这个迭代器对象，标签器纯粹是美学的（并且有详细的文档here）。第一步是转置 y_ds。这也是审美， 如果批次维度出现在之前，那么在绘图时我们不会在同一行看到所有的 pred1。 combined=True 和 skip_dims 表示我们不想迭代链或 x 维度中的值（默认情况下，ArviZ 从不迭代绘制。现在我们可以使用这个循环获取数组的迭代器 chain,x 维度，因此我们可以绘制 HDI，均值...并获得类似于问题中示例的图：

for ax,(_,sel,isel,values) in zip(axes.ravel(),iterator):
    az.plot_hdi(x,values,ax=ax)
    ax.plot(x,values.mean(axis=(0,1)))
    ax.set_title(labeller.make_label_vert("y",isel))
fig.tight_layout()

请注意，yaxis 仅在行之间共享，pred1 和 pred2 的比例是一个数量级不同的，即使它们与此绘图布局看起来相似。

旁注：我遇到了很多分歧，大多数rhat值都在推荐的1.01以上。我很确定 NUTS 没有收敛，你可能不应该相信这些结果。