16. 随机化回应调查#

16.1. 概述#

社会污名可能会阻止人们承认潜在的令人尴尬的行为或观点。

当人们不愿意参与关于个人敏感问题的抽样调查时，他们可能会拒绝参与。即使参与了，他们也可能会对敏感问题提供不真实的答案。

这些问题会导致选择偏差，给调查的解释和设计带来挑战。

为了说明社会科学家如何思考估计此类令人尴尬的行为和观点的普遍程度，本讲座描述了S. L. Warner [Warner, 1965]的一种经典方法。

Warner使用基础概率论构建了一种方法，在保护调查受访者个人隐私的同时，仍能估算出在一个群体中具有社会污名特征或从事社会污名活动的人数比例。

Warner的想法是在受访者的答案与调查制作者最终收到的信号之间添加噪声。

了解噪声的结构可以让受访者确信调查制作者无法观察到他的答案。

噪声注入程序的统计特性为受访者提供了合理的可否认性。

相关理念构成了现代差分隐私系统的基础。

(参见 https://baike.baidu.com/item/差分隐私/22415732)

16.2. Warner的策略#

像往常一样，让我们导入将要使用的Python模块。

import numpy as np
import pandas as pd

假设人群中的每个人要么属于 A 组，要么属于 B 组。

我们想要估计属于 A 组的人口比例 \(\pi\)，同时保护个人受访者的隐私。

Warner [Warner, 1965] 提出并分析了以下程序：

  • 从人群中有放回地抽取 \(n\) 个随机样本，并对每个人进行访谈。

  • 准备一个随机转盘，该转盘指向字母 A 的概率为 \(p\)，指向字母 B 的概率为 \((1-p)\)。

  • 每个受试者转动随机转盘，看到一个访谈者看不到的结果（A 或 B）。

  • 受试者说明自己是否属于转盘所指向的组。

  • 如果转盘指向受试者所属的组，受试者回答”是”；否则回答”否”。

  • 受试者如实回答问题。

Warner构建了一个用于估计总体中集合A所占比例的最大似然估计量。

设

  • \(\pi\) : 总体中A的真实概率

  • \(p\) : 指针指向A的概率

  • \(X_{i}=\begin{cases}1,\text{ 如果第}i\text{个受试者回答是}\\0,\text{ 如果第}i\text{个受试者回答否}\end{cases}\)

将样本集编号，使得前\(n_1\)个报告”是”，而后\(n-n_1\)个报告”否”。

样本集的似然函数为

(16.1)#\[ L=\left[\pi p + (1-\pi)(1-p)\right]^{n_{1}}\left[(1-\pi) p +\pi (1-p)\right]^{n-n_{1}} \]

似然函数的对数为：

(16.2)#\[ \log(L)= n_1 \log \left[\pi p + (1-\pi)(1-p)\right] + (n-n_{1}) \log \left[(1-\pi) p +\pi (1-p)\right] \]

关于 \(\pi\) 最大化对数似然函数的一阶必要条件是：

\[ \frac{(n-n_1)(2p-1)}{(1-\pi) p +\pi (1-p)}=\frac{n_1 (2p-1)}{\pi p + (1-\pi)(1-p)} \]

或

(16.3)#\[\pi p + (1-\pi)(1-p)=\frac{n_1}{n} \]

如果 \(p \neq \frac{1}{2}\)，则最大似然估计量（MLE）\(\pi\) 为：

(16.4)#\[ \hat{\pi}=\frac{p-1}{2p-1}+\frac{n_1}{(2p-1)n} \]

我们计算最大似然估计量 \(\hat \pi\) 的均值和方差为：

(16.5)#\[\begin{split} \begin{aligned} \mathbb{E}(\hat{\pi})&= \frac{1}{2 p-1}\left[p-1+\frac{1}{n} \sum_{i=1}^{n} \mathbb{E} X_i \right] \\ &=\frac{1}{2 p-1} \left[ p -1 + \pi p + (1-\pi)(1-p)\right] \\ &=\pi \end{aligned} \end{split}\]

以及

(16.6)#\[\begin{split} \begin{aligned} Var(\hat{\pi})&=\frac{n Var(X_i)}{(2p - 1 )^2 n^2} \\ &= \frac{\left[\pi p + (1-\pi)(1-p)\right]\left[(1-\pi) p +\pi (1-p)\right]}{(2p - 1 )^2 n^2}\\ &=\frac{\frac{1}{4}+(2 p^2 - 2 p +\frac{1}{2})(- 2 \pi^2 + 2 \pi -\frac{1}{2})}{(2p - 1 )^2 n^2}\\ &=\frac{1}{n}\left[\frac{1}{16(p-\frac{1}{2})^2}-(\pi-\frac{1}{2})^2 \right] \end{aligned} \end{split}\]

方程 (16.5) 表明 \(\hat{\pi}\) 是 \(\pi\) 的无偏估计量，而方程 (16.6) 告诉我们估计量的方差。

为了计算置信区间，首先将 (16.6) 重写为：

(16.7)#\[ Var(\hat{\pi})=\frac{\frac{1}{4}-(\pi-\frac{1}{2})^2}{n}+\frac{\frac{1}{16(p-\frac{1}{2})^2}-\frac{1}{4}}{n} \]

这个方程表明 \(\hat{\pi}\) 的方差可以表示为抽样方差和随机设备方差的和。

从上述表达式中我们可以发现：

  • 当 \(p\) 为 \(\frac{1}{2}\) 时，表达式 (16.1) 退化为一个常数。

  • 当 \(p\) 为 \(1\) 或 \(0\) 时，随机化估计退化为不含随机抽样的估计量。

我们将只讨论 \(p \in (\frac{1}{2},1)\) 的情况

（\(p \in (0,\frac{1}{2})\) 的情况是对称的）。

从表达式 (16.5) 和 (16.7) 我们可以推导出：

  • 随着 \(p\) 的增加，\(\hat{\pi}\) 的均方误差(MSE)减小。

16.3. 比较两种调查设计#

让我们比较前面的随机化回答法与一个简化的非随机化回答法。

在我们的非随机化回答法中，我们假设：

  • A组成员以概率 \(T_a\) 说真话，而B组成员以概率 \(T_b\) 说真话

  • \(Y_i\) 为 \(1\) 或 \(0\)，取决于样本中第 \(i\) 个成员的报告是否属于A组。

那么我们可以估计 \(\pi\) 为：

(16.8)#\[ \hat{\pi}=\frac{\sum_{i=1}^{n}Y_i}{n} \]

我们计算该估计量的期望值、偏差和方差为：

(16.9)#\[\begin{split} \begin{aligned} \mathbb{E}(\hat{\pi})&=\pi T_a + \left[ (1-\pi)(1-T_b)\right]\\ \end{aligned} \end{split}\]
(16.10)#\[\begin{split} \begin{aligned} Bias(\hat{\pi})&=\mathbb{E}(\hat{\pi}-\pi)\\ &=\pi [T_a + T_b -2 ] + [1- T_b] \\ \end{aligned} \end{split}\]
(16.11)#\[ \begin{align}\begin{aligned} \begin{aligned}\\Var(\hat{\pi})&=\frac{ \left[ \pi T_a + (1-\pi)(1-T_b)\right] \left[1- \pi T_a -(1-\pi)(1-T_b)\right] }{n} \end{aligned} \end{aligned}\end{align} \]

定义一个

\[ \text{MSE 比率}=\frac{\text{随机化的均方误差}}{\text{常规的均方误差}} \]

我们可以计算不同参数值下不同调查设计的MSE比率。

以下Python代码计算了我们想要观察的对象，以便在不同的\(\pi_A\)和\(n\)值下进行比较：

class Comparison:
    def __init__(self, A, n):
        self.A = A
        self.n = n
        TaTb = np.array([[0.95,  1], [0.9,   1], [0.7,    1], 
                         [0.5,   1], [1,  0.95], [1,    0.9], 
                         [1,   0.7], [1,   0.5], [0.95, 0.95], 
                         [0.9, 0.9], [0.7, 0.7], [0.5,  0.5]])
        self.p_arr = np.array([0.6, 0.7, 0.8, 0.9])
        self.p_map = dict(zip(self.p_arr, [f"MSE 比率: p = {x}" for x in self.p_arr]))
        self.template = pd.DataFrame(columns=self.p_arr)
        self.template[['T_a','T_b']] = TaTb
        self.template['Bias'] = None
    
    def theoretical(self):
        A = self.A
        n = self.n
        df = self.template.copy()
        df['Bias'] = A * (df['T_a'] + df['T_b'] - 2) + (1 - df['T_b'])
        for p in self.p_arr:
            df[p] = (1 / (16 * (p - 1/2)**2) - (A - 1/2)**2) / n / \
                    (df['Bias']**2 + ((A * df['T_a'] + (1 - A) * (1 - df['T_b'])) * (1 - A * df['T_a'] - (1 - A) * (1 - df['T_b'])) / n))
            df[p] = df[p].round(2)
        df = df.set_index(["T_a", "T_b", "Bias"]).rename(columns=self.p_map)
        return df
        
    def MCsimulation(self, size=1000, seed=123456):
        A = self.A
        n = self.n
        df = self.template.copy()
        np.random.seed(seed)
        sample = np.random.rand(size, self.n) <= A
        random_device = np.random.rand(size, n)
        mse_rd = {}
        for p in self.p_arr:
            spinner = random_device <= p
            rd_answer = sample * spinner + (1 - sample) * (1 - spinner)
            n1 = rd_answer.sum(axis=1)
            pi_hat = (p - 1) / (2 * p - 1) + n1 / n / (2 * p - 1)
            mse_rd[p] = np.sum((pi_hat - A)**2)
        for inum, irow in df.iterrows():
            truth_a = np.random.rand(size, self.n) <= irow.T_a
            truth_b = np.random.rand(size, self.n) <= irow.T_b
            trad_answer = sample * truth_a + (1 - sample) * (1 - truth_b)
            pi_trad = trad_answer.sum(axis=1) / n
            df.loc[inum, 'Bias'] = pi_trad.mean() - A
            mse_trad = np.sum((pi_trad - A)**2)
            for p in self.p_arr:
                df.loc[inum, p] = (mse_rd[p] / mse_trad).round(2)
        df = df.set_index(["T_a", "T_b", "Bias"]).rename(columns=self.p_map)
        return df

让我们使用以下参数值来运行代码

  • \(\pi_A=0.6\)

  • \(n=1000\)

我们可以使用上述公式理论上生成MSE比率。

我们也可以对MSE比率进行蒙特卡洛模拟。

cp1 = Comparison(0.6, 1000)
df1_theoretical = cp1.theoretical()
df1_theoretical
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.03 5.45 1.36 0.60 0.33
0.90 1.00 -0.06 1.62 0.40 0.18 0.10
0.70 1.00 -0.18 0.19 0.05 0.02 0.01
0.50 1.00 -0.30 0.07 0.02 0.01 0.00
1.00 0.95 0.02 9.82 2.44 1.08 0.60
0.90 0.04 3.41 0.85 0.37 0.21
0.70 0.12 0.43 0.11 0.05 0.03
0.50 0.20 0.16 0.04 0.02 0.01
0.95 0.95 -0.01 18.25 4.54 2.00 1.11
0.90 0.90 -0.02 9.70 2.41 1.06 0.59
0.70 0.70 -0.06 1.62 0.40 0.18 0.10
0.50 0.50 -0.10 0.61 0.15 0.07 0.04 
df1_mc = cp1.MCsimulation()
df1_mc
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.030060 5.76 1.36 0.63 0.35
0.90 1.00 -0.060045 1.73 0.41 0.19 0.1
0.70 1.00 -0.179530 0.21 0.05 0.02 0.01
0.50 1.00 -0.300077 0.07 0.02 0.01 0.0
1.00 0.95 0.019770 10.59 2.5 1.15 0.64
0.90 0.040050 3.63 0.86 0.39 0.22
0.70 0.120052 0.46 0.11 0.05 0.03
0.50 0.199746 0.17 0.04 0.02 0.01
0.95 0.95 -0.010137 18.65 4.41 2.02 1.12
0.90 0.90 -0.020103 10.48 2.48 1.14 0.63
0.70 0.70 -0.060488 1.71 0.4 0.19 0.1
0.50 0.50 -0.099341 0.66 0.16 0.07 0.04

理论计算很好地预测了蒙特卡洛结果。

我们看到在许多情况下，特别是当偏差不小时，随机抽样方法的均方误差比非随机抽样方法要小。

随着 \(p\) 的增加，这些差异变得更大。

通过调整参数 \(\pi_A\) 和 \(n\)，我们可以研究不同情况下的结果。

例如，对于Warner [Warner, 1965]描述的另一种情况：

  • \(\pi_A=0.5\)

  • \(n=1000\)

我们可以使用以下代码

cp2 = Comparison(0.5, 1000)
df2_theoretical = cp2.theoretical()
df2_theoretical
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.025 7.15 1.79 0.79 0.45
0.90 1.00 -0.050 2.27 0.57 0.25 0.14
0.70 1.00 -0.150 0.27 0.07 0.03 0.02
0.50 1.00 -0.250 0.10 0.02 0.01 0.01
1.00 0.95 0.025 7.15 1.79 0.79 0.45
0.90 0.050 2.27 0.57 0.25 0.14
0.70 0.150 0.27 0.07 0.03 0.02
0.50 0.250 0.10 0.02 0.01 0.01
0.95 0.95 0.000 25.00 6.25 2.78 1.56
0.90 0.90 0.000 25.00 6.25 2.78 1.56
0.70 0.70 0.000 25.00 6.25 2.78 1.56
0.50 0.50 0.000 25.00 6.25 2.78 1.56 
df2_mc = cp2.MCsimulation()
df2_mc
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.025230 7.0 1.69 0.75 0.44
0.90 1.00 -0.050279 2.23 0.54 0.24 0.14
0.70 1.00 -0.149866 0.27 0.07 0.03 0.02
0.50 1.00 -0.250211 0.1 0.02 0.01 0.01
1.00 0.95 0.024410 7.38 1.78 0.79 0.46
0.90 0.049839 2.26 0.54 0.24 0.14
0.70 0.149769 0.27 0.07 0.03 0.02
0.50 0.249851 0.1 0.02 0.01 0.01
0.95 0.95 -0.000260 24.29 5.86 2.59 1.52
0.90 0.90 -0.000109 25.73 6.2 2.74 1.61
0.70 0.70 -0.000439 25.75 6.21 2.74 1.61
0.50 0.50 0.000768 24.91 6.01 2.65 1.56

我们还可以重新审视Warner [Warner, 1965]结论部分的一个计算，其中

  • \(\pi_A=0.6\)

  • \(n=2000\)

我们使用以下代码

cp3 = Comparison(0.6, 2000)
df3_theoretical = cp3.theoretical()
df3_theoretical
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.03 3.05 0.76 0.33 0.19
0.90 1.00 -0.06 0.84 0.21 0.09 0.05
0.70 1.00 -0.18 0.10 0.02 0.01 0.01
0.50 1.00 -0.30 0.03 0.01 0.00 0.00
1.00 0.95 0.02 6.03 1.50 0.66 0.37
0.90 0.04 1.82 0.45 0.20 0.11
0.70 0.12 0.22 0.05 0.02 0.01
0.50 0.20 0.08 0.02 0.01 0.00
0.95 0.95 -0.01 14.12 3.51 1.55 0.86
0.90 0.90 -0.02 5.98 1.49 0.66 0.36
0.70 0.70 -0.06 0.84 0.21 0.09 0.05
0.50 0.50 -0.10 0.31 0.08 0.03 0.02 
df3_mc = cp3.MCsimulation()
df3_mc
MSE 比率: p = 0.6 MSE 比率: p = 0.7 MSE 比率: p = 0.8 MSE 比率: p = 0.9
T_a T_b Bias
0.95 1.00 -0.030316 3.27 0.8 0.34 0.19
0.90 1.00 -0.060352 0.91 0.22 0.09 0.05
0.70 1.00 -0.180087 0.11 0.03 0.01 0.01
0.50 1.00 -0.299849 0.04 0.01 0.0 0.0
1.00 0.95 0.019734 6.7 1.64 0.69 0.39
0.90 0.039766 2.01 0.49 0.21 0.12
0.70 0.119789 0.24 0.06 0.02 0.01
0.50 0.200138 0.09 0.02 0.01 0.0
0.95 0.95 -0.010475 14.78 3.61 1.53 0.85
0.90 0.90 -0.020373 6.32 1.54 0.65 0.36
0.70 0.70 -0.059945 0.92 0.23 0.1 0.05
0.50 0.50 -0.100103 0.34 0.08 0.03 0.02

显然，随着 \(n\) 的增加，随机化回应法在更多情况下表现更好。

16.4. 结束语#

这个QuantEcon讲座描述了一些其他的随机化回应调查方法。

该讲座介绍了Lars Ljungqvist [Ljungqvist, 1993]对这些替代方案进行的功利主义分析。