21. Randomized Response Surveys

21.1. Overview

Social stigmas can inhibit people from confessing potentially embarrassing activities or opinions.

When people are reluctant to participate a sample survey about personally sensitive issues, they might decline to participate, and even if they do participate, they might choose to provide incorrect answers to sensitive questions.

These problems induce selection biases that present challenges to interpreting and designing surveys.

To illustrate how social scientists have thought about estimating the prevalence of such embarrassing activities and opinions, this lecture describes a classic approach of S. L. Warner [War65].

Warner used elementary probability to construct a way to protect the privacy of individual respondents to surveys while still estimating the fraction of a collection of individuals who have a socially stigmatized characteristic or who engage in a socially stigmatized activity.

Warner’s idea was to add noise between the respondent’s answer and the signal about that answer that the survey maker ultimately receives.

Knowing about the structure of the noise assures the respondent that the survey maker does not observe his answer.

Statistical properties of the noise injection procedure provide the respondent plausible deniability.

Related ideas underlie modern differential privacy systems.

(See https://en.wikipedia.org/wiki/Differential_privacy)

21.2. Warner’s Strategy

As usual, let’s bring in the Python modules we’ll be using.

import numpy as np
import pandas as pd

Suppose that every person in population either belongs to Group A or Group B.

We want to estimate the proportion \(\pi\) who belong to Group A while protecting individual respondents’ privacy.

Warner [War65] proposed and analyzed the following procedure.

• A random sample of \(n\) people is drawn with replacement from the population and each person is interviewed.

• Draw \(n\) random samples from the population with replacement and interview each person.

• Prepare a random spinner that with \(p\) probability points to the Letter A and with \((1-p)\) probability points to the Letter B.

• Each subject spins a random spinner and sees an outcome (A or B) that the interviewer does not observe.

• The subject states whether he belongs to the group to which the spinner points.

• If the spinner points to the group that the spinner belongs, the subject reports “yes”; otherwise he reports “no”.

• The subject answers the question truthfully.

Warner constructed a maximum likelihood estimators of the proportion of the population in set A.

Let

• \(\pi\) : True probability of A in the population

• \(p\) : Probability that the spinner points to A

• \(X_{i}=\begin{cases}1,\text{ if the } i\text{th} \ \text{ subject says yes}\\0,\text{ if the } i\text{th} \ \text{ subject says no}\end{cases}\)

Index the sample set so that the first \(n_1\) report “yes”, while the second \(n-n_1\) report “no”.

The likelihood function of a sample set is

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

The log of the likelihood function is:

(21.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] \]

The first-order necessary condition for maximizing the log likelihood function with respect to \(\pi\) is:

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

or

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

If \(p \neq \frac{1}{2}\), then the maximum likelihood estimator (MLE) of \(\pi\) is:

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

We compute the mean and variance of the MLE estimator \(\hat \pi\) to be:

(21.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}\]

and

(21.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}\]

Equation (21.5) indicates that \(\hat{\pi}\) is an unbiased estimator of \(\pi\) while equation (21.6) tell us the variance of the estimator.

To compute a confidence interval, first rewrite (21.6) as:

(21.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} \]

This equation indicates that the variance of \(\hat{\pi}\) can be represented as a sum of the variance due to sampling plus the variance due to the random device.

From the expressions above we can find that:

• When \(p\) is \(\frac{1}{2}\), expression (21.1) degenerates to a constant.

• When \(p\) is \(1\) or \(0\), the randomized estimate degenerates to an estimator without randomized sampling.

We shall only discuss situations in which \(p \in (\frac{1}{2},1)\)

(a situation in which \(p \in (0,\frac{1}{2})\) is symmetric).

From expressions (21.5) and (21.7) we can deduce that:

• The MSE of \(\hat{\pi}\) decreases as \(p\) increases.

21.3. Comparing Two Survey Designs

Let’s compare the preceding randomized-response method with a stylized non-randomized response method.

In our non-randomized response method, we suppose that:

• Members of Group A tells the truth with probability \(T_a\) while the members of Group B tells the truth with probability \(T_b\)

• \(Y_i\) is \(1\) or \(0\) according to whether the sample’s \(i\text{th}\) member’s report is in Group A or not.

Then we can estimate \(\pi\) as:

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

We calculate the expectation, bias, and variance of the estimator to be:

(21.9)\[\begin{split} \begin{aligned} \mathbb{E}(\hat{\pi})&=\pi T_a + \left[ (1-\pi)(1-T_b)\right]\\ \end{aligned} \end{split}\]

(21.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}\]

(21.11)\[ \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} \]

It is useful to define a

\[ \text{MSE Ratio}=\frac{\text{Mean Square Error Randomized}}{\text{Mean Square Error Regular}} \]

We can compute MSE Ratios for different survey designs associated with different parameter values.

The following Python code computes objects we want to stare at in order to make comparisons under different values of \(\pi_A\) and \(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 Ratio: 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

Let’s put the code to work for parameter values

• \(\pi_A=0.6\)

• \(n=1000\)

We can generate MSE Ratios theoretically using the above formulas.

We can also perform Monte Carlo simulations of a MSE Ratio.

cp1 = Comparison(0.6, 1000)
df1_theoretical = cp1.theoretical()
df1_theoretical

			MSE Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: 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 Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: p = 0.9
T_a	T_b	Bias
0.95	1.00	-0.030060000000000087	5.76	1.36	0.63	0.35
0.90	1.00	-0.060045000000000015	1.73	0.41	0.19	0.1
0.70	1.00	-0.17952999999999997	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.019769999999999954	10.59	2.5	1.15	0.64
	0.90	0.04005000000000014	3.63	0.86	0.39	0.22
	0.70	0.12005199999999994	0.46	0.11	0.05	0.03
	0.50	0.1997460000000001	0.17	0.04	0.02	0.01
0.95	0.95	-0.010137000000000063	18.65	4.41	2.02	1.12
0.90	0.90	-0.020103000000000093	10.48	2.48	1.14	0.63
0.70	0.70	-0.060487999999999875	1.71	0.4	0.19	0.1
0.50	0.50	-0.09934100000000001	0.66	0.16	0.07	0.04

The theoretical calculations do a good job of predicting Monte Carlo results.

We see that in many situations, especially when the bias is not small, the MSE of the randomized-sampling methods is smaller than that of the non-randomized sampling method.

These differences become larger as \(p\) increases.

By adjusting parameters \(\pi_A\) and \(n\), we can study outcomes in different situations.

For example, for another situation described in Warner [War65]:

• \(\pi_A=0.5\)

• \(n=1000\)

we can use the code

cp2 = Comparison(0.5, 1000)
df2_theoretical = cp2.theoretical()
df2_theoretical

			MSE Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: 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 Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: p = 0.9
T_a	T_b	Bias
0.95	1.00	-0.02523000000000003	7.0	1.69	0.75	0.44
0.90	1.00	-0.05027900000000002	2.23	0.54	0.24	0.14
0.70	1.00	-0.14986600000000005	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.024410000000000043	7.38	1.78	0.79	0.46
	0.90	0.04983899999999997	2.26	0.54	0.24	0.14
	0.70	0.14976900000000004	0.27	0.07	0.03	0.02
	0.50	0.24985100000000005	0.1	0.02	0.01	0.01
0.95	0.95	-0.00025999999999998247	24.29	5.86	2.59	1.52
0.90	0.90	-0.00010899999999997023	25.73	6.2	2.74	1.61
0.70	0.70	-0.0004390000000000782	25.75	6.21	2.74	1.61
0.50	0.50	0.0007679999999999909	24.91	6.01	2.65	1.56

We can also revisit a calculation in the concluding section of Warner [War65] in which

• \(\pi_A=0.6\)

• \(n=2000\)

We use the code

cp3 = Comparison(0.6, 2000)
df3_theoretical = cp3.theoretical()
df3_theoretical

			MSE Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: 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 Ratio: p = 0.6	MSE Ratio: p = 0.7	MSE Ratio: p = 0.8	MSE Ratio: p = 0.9
T_a	T_b	Bias
0.95	1.00	-0.03031649999999997	3.27	0.8	0.34	0.19
0.90	1.00	-0.0603515	0.91	0.22	0.09	0.05
0.70	1.00	-0.18008749999999996	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.01973400000000003	6.7	1.64	0.69	0.39
	0.90	0.03976600000000008	2.01	0.49	0.21	0.12
	0.70	0.11978900000000003	0.24	0.06	0.02	0.01
	0.50	0.2001385	0.09	0.02	0.01	0.0
0.95	0.95	-0.010474500000000053	14.78	3.61	1.53	0.85
0.90	0.90	-0.020373499999999933	6.32	1.54	0.65	0.36
0.70	0.70	-0.059945499999999985	0.92	0.23	0.1	0.05
0.50	0.50	-0.10010349999999996	0.34	0.08	0.03	0.02

Evidently, as \(n\) increases, the randomized response method does better performance in more situations.

21.4. Concluding Remarks

This QuantEcon lecture describes some alternative randomized response surveys.

That lecture presents a utilitarian analysis of those alternatives conducted by Lars Ljungqvist [Lju93].