在 pymc3 中创建三级逻辑回归模型

Posted 2023-02-19

【中文标题】在 pymc3 中创建三级逻辑回归模型

【英文标题】：Creating a three-level logistic regression model in pymc3

【发布时间】：2017-04-13 22:04:18

【问题描述】：

我正在尝试在 pymc3 中创建一个三级逻辑回归模型。有顶层、中层和个体层，其中中层系数是根据顶层系数估计的。但是，我很难为中级指定正确的数据结构。

这是我的代码：

with pm.Model() as model:
    # Hyperpriors
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    

    # Priors
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
    mid_level = [pm.Normal('mid_level_'.format(j),
                           mu=top_level[mid_to_top_idx[j]],
                           tau=mid_level_tau)
                 for j in range(k_mid)]

    intercept = pm.Normal('intercept', mu=0., sd=100.)

    # Model prediction
    yhat = pm.invlogit(mid_level[mid_to_bot_idx] + intercept)

    # Likelihood
    yact = pm.Bernoulli('yact', p=yhat, observed=y)

我收到错误 "only integer arrays with one element can be converted to an index"（第 16 行），我认为这与 mid_level 变量是一个列表，而不是一个正确的 pymc 容器这一事实有关。 （我也没有在pymc3源代码中看到Container类。）

任何帮助将不胜感激。

编辑：添加一些模拟数据

y = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0])
mid_to_bot_idx = np.array([0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 2])
mid_to_top_idx = np.array([0, 0, 1, 1])
k_top = 2
k_mid = 4

编辑#2：

似乎有几种不同的方法可以解决这个问题，尽管没有一种方法是完全令人满意的：

1) 可以将模型重构为：

with pm.Model() as model:
    # Hyperpriors
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    

    # Priors
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
    mid_level = pm.Normal('mid_level', mu=0., tau=mid_level_tau, shape=k_top)
    intercept = pm.Normal('intercept', mu=0., sd=100.)

    # Model prediction
    yhat = pm.invlogit(top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept)

    # Likelihood
    yact = pm.Bernoulli('yact', p=yhat, observed=y)

这似乎可行，尽管我不知道如何将其扩展到所有中级组的中级方差不是恒定的情况。

2) 可以使用 theano.tensor.stack 将中级系数包装成 Theano 张量：即，

import theano.tensor as tt
mid_level = tt.stack([pm.Normal('mid_level_'.format(j),
                           mu=top_level[mid_to_top_idx[j]],
                           tau=mid_level_tau)
                 for j in range(k_mid)])

但这似乎在我的实际数据集（30k 观察）上运行得非常缓慢，并且它使绘制不方便（每个 mid_level 系数都使用pm.traceplot 获得自己的跟踪）。

无论如何，我们将不胜感激开发人员的一些建议/意见。

【问题讨论】：

@gung 现在看起来还可以吗？

谢谢，太好了。关于 Python 编码的问题不在此处讨论，但可以在 Stack Overflow 上讨论。如果您等待，我们将尝试将您的问题迁移到那里。

我不同意这是题外话：这不是一个通用的 Python 编码问题。这个问题是关于使用成熟的 python 统计包实现统计模型——答案很可能是以不同的方式表示模型。

我相信这个问题属于 stats.stackexchange.com

你的模型中没有预测变量，应该是yhat = pm.invlogit(top_level[top_to_bot_idx] * x + mid_level[mid_to_bot_idx] * x + intercept)吗？

【参考方案1】：

原来解决方案很简单：似乎任何分布（如pm.Normal）都可以接受均值向量作为参数，因此替换此行

mid_level = [pm.Normal('mid_level_'.format(j),
                       mu=top_level[mid_to_top_idx[j]],
                       tau=mid_level_tau)
             for j in range(k_mid)]

有了这个

mid_level = pm.Normal('mid_level',
                       mu=top_level[mid_to_top_idx],
                       tau=mid_level_tau,
                       shape=k_mid)

有效。同样的方法也可用于为每个中级组指定单独的标准差。

编辑：修正错字

【参考方案2】：

变化不大（注意我把 yhat 改成了 theta）：

theta = pm.Deterministic( 'theta', pm.invlogit(sum(mid_level[i] for i in mid_to_bot_idx)+intercept) )
yact = pm.Bernoulli( 'yact', p=theta, observed=y )

【讨论】：

我认为这不是我想要的（尽管我可能会弄错）。这似乎将所有系数相加，以获得所有观测值的单个 theta 值。我希望每次观察都有不同的 thetas，具体取决于 top_level 和 mid_level。换句话说，theta 应该是与 y 长度相同的向量，而不是标量。例如，看这个模型：nbviewer.jupyter.org/github/aflaxman/pymc-examples/blob/master/…

@vbox 是的，我倾向于同意你的观点。您的原始代码中的 mid_level 数组由 mid_to_bot_idx 数组简单地重新排序（并且一些值重复）。我完全按照您的代码中显示的方式实现了。

invlogit 的参数通常类似于 (x * beta + intercept) 其中 x 是特征，beta 是回归系数，intercept 是偏差项。

【参考方案3】：

在您的问题中，您声明了yhat。 您可以避免它并将方程传递给Bernoulli 的logit_p 参数。

注意 - 您可以传递 p 或 logit_p。

就我而言，使用logit_p 可以加快我的采样过程。

代码-

# Likelihood
yact = pm.Bernoulli('yact', logit_p=top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept, observed=y)