人脸识别之数据集中 PI20 和 CFMT 之间关联的模型预测贝叶斯（Python+论文代码实现）

代码文件（联系作者点击这里末尾）

代码文件描述如下：

subjective_objective.ipynb和：这分别是实际的笔记本和 Web 浏览器友好的只读版本。此笔记本读取数据，执行一些预处理，并包含论文中使用的模型规范。它还创建了手稿中使用的图形。笔记本全程都有注释。subjective_objective.html

数据文件

本文提供了许多数据文件，但它们作为原始数据的审计跟踪存在，这些数据被合并为用于分析的单个数据。All_data.csv

All_data.csv：此文件包含来自 10 个不同研究的组合数据，用于分析。它包含以下列：
- ID- 数据集来自的研究的名称。有关详细信息，请参阅论文。每个标签都有一个后缀“norm”或“DP”，用于创建组标签 - 所有 DP 参与者都自我引用。
- Age- 参与者年龄，以岁为单位。
- Gender- 参与者自我报告的性别。
- PI20- 参与者在 PI20 问卷上的总分。
- CFMT %- 参与者获得的 CFMT 上的正确试验百分比。

Datasets 文件夹

此文件夹包含来自其他 OSF 页面或通过与作者的个人通信的原始数据，从中获取相关值（年龄、性别、PI20 和 CFMT 分数）。他们是：

Burns 2024.csv
Gray 2017 a and b.xlsx
shah_2015.xlsx
Tagliente 2022 (emailed).csv
Tsantini_2021.xlsx

我们这里的目标是测试 CFMT（人脸识别的客观测量）和 PI20（人脸识别的主观测量）之间的关联在自我指责为有问题的个体群体和规范样本之间是否不同。我们将使用贝叶斯分层建模来拟合线性回归并恢复斜率。

# Imports and settings
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
import scipy.stats as st
from formulae import design_matrices

# Random seed and plotting defaults
rng = np.random.default_rng(35)
plt.style.use('arviz-docgrid')
plt.rcParams['axes.axisbelow'] = True

首先，我们读入数据，清理自我报告的性别，从 ID 名称以及 z 分数年龄、PI20 和 CFMT 分数中创建一个分组变量。我们还对数据集进行了快速可视化。

# read gathered data
df = pd.read_csv('All_data.csv')

# Z score, relabel
df = df.assign(age_c = st.zscore(df['Age'], nan_policy='omit'),
               CFMT_c = st.zscore(df['CFMT %']),
               PI20_c = st.zscore(df['PI20']),
               Group_l = np.where(df['ID'].str.contains('norm'),
                                  'Normative', 'Self-Refer'),
               Gender=df['Gender'].str.lower().str.strip().str.replace(' ', '')
              )

# Plot
fig, ax = plt.subplots(1, 1, figsize=(20, 10))
sns.scatterplot(data=df, x='PI20_c', y='CFMT_c', 
                hue='Group_l', style='Group_l',
                size='Gender', alpha=.8, ax=ax)

<Axes: xlabel='PI20_c', ylabel='CFMT_c'>

并打印出 dataframe：

# Display
display(df.head(), df.tail())

接下来，我们使用 PyMC 构建一个分层模型，该模型还考虑了每个组数据集中的变化。它包括每个数据集中 PI20 的随机截距和斜率。首先，我们创建 design matix 并提供一些单独的变量，以使模型更易于编写。

# Use formulae to derive the design matrix
df = df.dropna(how='any').reset_index(drop=True).rename(columns=lambda x: x.replace(' %', '')) # drops 5 people missing age and sex data
dm = design_matrices('CFMT_c ~ Gender + age_c + PI20_c * Group_l', data=df) # creates the design matrix

X = dm.common.as_dataframe().drop(columns=['Intercept', 'PI20_c']) # drops intercept and PI20 which is isolated below
PI20_c = dm.common.as_dataframe()['PI20_c'] # PI20 predictor alone
Y = dm.response.as_dataframe().squeeze() # CFMT Y
# Coords
dataset_idx, dataset_label = df['ID'].factorize()

c = {'dataset': dataset_label, 
     'obs': X.index,
     'preds': X.columns}

with pm.Model(coords=c) as bayes_mod:

    # Set data
    Xm = pm.Data('Xm', X.to_numpy(), dims=('obs', 'preds'))
    pi20_c = pm.Data('pi20', PI20_c.to_numpy(), dims='obs')
    Ym = pm.Data('Ym', Y.to_numpy(), dims='obs')

    # Priors for intercept
    β0 = pm.Normal('β0', mu=0, sigma=10)

    # Priors for coefficients, NOT the PI20
    β = pm.Normal('β', mu=0, sigma=10, dims='preds')

    # Prior for the PI20 fixed effect
    βpi20 = pm.Normal('βpi20', mu=0, sigma=10)

    # Priors for random intercepts (first) and slopes (second)
    αi = pm.ZeroSumNormal('αi', sigma=pm.InverseGamma('dataset_ασ', mu=1, sigma=1), dims='dataset')
    βi = pm.ZeroSumNormal('βi', sigma=pm.InverseGamma('dataset_βσ', mu=1, sigma=1), dims='dataset')

    # Prior for sigma, overall model error
    σ = pm.HalfNormal('σ', sigma=4)

    # Linear combination
    μ = pm.Deterministic('μ',
                         (β0 + αi[dataset_idx]) + # intercept, plus random intercepts
                         Xm @ β + # age, gender, and interaction coefficient here
                         (βpi20 + βi[dataset_idx]) * pi20_c, # and the PI20 fixed effect and its random slopes
                         dims='obs')
    
    # Normal likelihood
    pm.Normal('y', mu=μ, sigma=σ, observed=Ym, dims='obs')

    # Not part of the model, but this makes a set of predictions after holding constant
    # age and gender - it includes random effects though
    linpred = pm.Deterministic('linpred',
                               β0 + αi[dataset_idx] + (βpi20 + βi[dataset_idx]) * pi20_c,
                               dims='obs')

我们可以可视化模型的结构：

在 [6] 中：

bayes_mod.to_graphviz()

最后，我们可以对模型进行采样，获得参数的后验。我们还收集后验预测和先验预测仅用于诊断。

接下来，我们通过过滤掉任何大于 1 的 RHat 来检查是否所有链都收敛了。如果此处没有，则模型已收敛。

# Check results
az.summary(idata, var_names=['β'], filter_vars='like')

在对模型进行采样时，我们收集了一组预测，这些预测说明了每个数据集中 PI20 和 CFMT 之间关联的模型预测。请注意，这与交互无关，这纯粹是模型在每个数据集中 PI20 和 CFMT 之间学习的变化，这很有趣，但不是关键 - PI20 的固定效应是从这种变化中学到的，因为我们指示模型这样做。可视化如下：

# Visualise each groups lines for beta
def plotter(**kwargs):

    # Get data, axis
    d = kwargs.pop('data')
    ax = plt.gca()

    # Plot the association between this X and Y
    locs = d.index

    x = idata['constant_data']['pi20'].isel(obs=locs)
    y = idata['posterior']['linpred'].isel(obs=locs)

    az.plot_hdi(x, y, ax=ax, color='black')
    ax.plot(x, y.mean(('chain', 'draw')), color='white', linewidth=2.5)


(
    sns.FacetGrid(data=df, col='ID', col_wrap=2, aspect=2)
    .map_dataframe(sns.scatterplot, y='CFMT_c', x='PI20_c', s=200, alpha=.2)
    .map_dataframe(plotter)
    .savefig('ooo.png')

我们现在可以使用该模型进行一系列感兴趣的预测。

要计算每组的斜率，我们只需要将 PI20 和交互系数相加（自我报告组），或者单独取 PI20 斜率（标准组）。
我们可以使用该模型对一系列 PI20 分数的 CFMT 分数进行预测，包括规范和自我报告的个体。通过这种方式，我们可以检查每个 PI20 级别的差异。这个反事实世界使我们能够看到 PI20 分数范围的各组差异。

# Get index locations of the required coefficients
sex_loc = X.columns.get_loc('Gender[male]')
group_loc = X.columns.get_loc('Group_l[Self-Refer]')
interact_loc = X.columns.get_loc('PI20_c:Group_l[Self-Refer]')
# Make predictions across range of PI20
with bayes_mod:

    XNew = pm.Data('XNew', np.linspace(df['PI20_c'].min(), df['PI20_c'].max()))

    # Make predictions for normative and self refer
    normative_mu = pm.Deterministic('normative_mu', β0 + β[sex_loc]*.5 + β[group_loc]*0 + βpi20 * XNew)
    selfrefer_mu = pm.Deterministic('selfrefer_mu', β0 + β[sex_loc]*.5 + β[group_loc]*1 + βpi20 * XNew + β[interact_loc] * XNew)

    # Get difference
    mean_contrast = pm.Deterministic('mean_contrast', selfrefer_mu - normative_mu) 

    # Sample using the posterior
    pm.sample_posterior_predictive(idata, predictions=True, 
                                   var_names=['normative_mu', 'selfrefer_mu', 'mean_contrast'], 
                                   extend_inferencedata=True, random_seed=rng)

让我们看看这一切是什么样子的！

顶部图是与 Models （模型）视图重叠的数据。
中间的图是斜率估计值的后验图及其差值。
底部图是 PI20 分数范围内各组之间的差异。

最后放大一点，我们可以查看在自我引用 PI20 Z 分数的最低和最高值下，自我引用个体的 CFMT 分数高于正常个体的概率。

# Get the min and max of PI20 per group
df.groupby('Group_l').agg({'PI20_c': ['min', 'max']})

最后，绘制差异图，并恢复自荐组较低的概率。由于对比是规范的 - 自我指涉，那么低于零的后验比例（即自我指涉分数小于规范分数）如下所示。

# visualise
az.plot_posterior(contrast, group='predictions', 
                  var_names=['mean_contrast2'], 
                  ref_val=0, ref_val_color='black')

array([<Axes: title={'center': 'mean_contrast2\n0'}>,
       <Axes: title={'center': 'mean_contrast2\n1'}>], dtype=object)

基于贝叶斯分层模型的CFMT与PI20关联性研究：自我报告群体与规范样本的对比分析

1. 研究设计与方法

1.1 数据准备与预处理

数据来源：整合多中心研究数据（含Tsantani_2021、Burns_2023等10个数据集），共1623个样本。
变量处理：
- 标准化：对年龄、CFMT正确率、PI20得分进行Z-score标准化（age_c, CFMT_c, PI20_c）。
- 分组定义：根据ID字段划分"Normative"（规范样本）和"Self-Refer"（自我报告人脸识别困难群体）。
- 性别编码：统一性别标签为小写并去除空格（male/female/non-binary等）。
数据可视化：通过散点图展示PI20_c与CFMT_c的组间分布差异（图1）。

1.2 贝叶斯分层线性模型构建

模型公式：

python

CFMT_c ~ Gender + age_c + PI20_c * Group_l + (1 + PI20_c | ID)

关键组件：
- 固定效应：性别、年龄、PI20主效应及PI20×组别交互项。
- 随机效应：数据集级别的随机截距（αi）和PI20斜率（βi），采用零和正态先验约束。
- 先验设定：
  - 固定效应系数：N(0,10)
  - 随机效应方差：InverseGamma(1,1)
  - 误差项：HalfNormal(4)
计算细节：
- 使用PyMC v5.16.1进行马尔可夫链蒙特卡洛采样（4链，2000调优+5000采样迭代）。
- 收敛性检验：所有参数Rhat <1.01，ESS>4000，表明采样充分。

2. 主要发现

2.1 模型整体表现

解释力：模型可解释CFMT_c 34.5%的方差（R²=0.345, 95% HDI [0.315,0.373]）。

关键固定效应：

参数	后验均值	95% HDI	显著性判断
PI20_c主效应 (βpi20)	-0.471	[-0.637, -0.299]	显著负向关联
组别主效应 (Group_l)	-0.328	[-0.919, 0.279]	无显著组间差异
PI20×组别交互项	-0.011	[-0.436, 0.403]	交互效应不显著

2.2 随机效应分析

数据集异质性：
- 随机斜率标准差：0.266 (95% HDI [0.134,0.423])，表明PI20-CFMT关联强度存在跨数据集变异。
- 典型数据集对比（图2）：
  - Tsantani_2021_norm：斜率+0.073 (HDI [-0.128,0.287])
  - Burns_2023_DP：斜率+0.066 (HDI [-0.215,0.364])

2.3 反事实预测分析

组间对比曲线（图3）：
- 在PI20_c全范围内，自我报告组CFMT预测值始终低于规范组（均值差=-0.328）。
- 差异幅度随PI20升高而扩大（从-0.2 SD到-0.5 SD）。
极端值分析：
- 当PI20_c=2.96（自我报告组上限）时：
  - 规范组预测CFMT_c= -0.89 (HDI [-1.21, -0.57])
  - 自我报告组预测CFMT_c= -1.15 (HDI [-1.62, -0.71])
  - 组间差异概率P(Self-Refer < Normative)=83.7%