Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,480 characters · 27 sections · 49 citation commands
Modelling with Sensitive Variables
\nocite{reguly2024doc_app}
Over the last decade, there has been an explosion in the amount of data collected on individuals, firms, the economy, and society. Governments and the private sector gather large amounts of information that researchers and analysts use on a daily basis. Privacy constraints to share these data or willingness to respond in case of surveys, often prevent access to sensitive data in its original form. We introduce here a novel sampling method that naturally protects individual privacy or increases the willingness to answer through discretization, while minimizes information loss from an estimation perspective. Income is perhaps one of the most commonly used sensitive variable in economics. Instead of providing the sensitive data itself or asking directly for the actual weekly personal income, we work with discrete categories such as `weekly personal income below $\$100$, between $\$100$ and $\$400$, or above $\$400$. It is well known in the econometric literature that this simple discretization leads to modeling problems (see, e.g., hsiao1983regression, manski2002inference, beresteanu2008asymptotic, beresteanu2011discretized, pacini2019interval, or abrevaya2020interval_fe), which we solve by using multiple discretization schemes. In the context of income, instead of using only one discretization scheme, we use multiple discretizations by changing the interval boundaries. For example, discretization scheme (2): below $\$50$, between $\$50$ and $\$350$ or above $\$350$; discretization scheme (3): below $\$150$, between $\$150$ and $\$450$ or above $\$450$, etc. We call the method split sampling,\footnote{The term split sampling in this paper is not related to the technique occasionally used in chromatography (see e.g., split1977gas1) or cross-validation methods in machine learning, which splits the initial sample into folds.} and show convergence in distribution for the discretized sensitive variable to the original unobserved variable as both the number of observations and the number of discretization schemes go to infinity. Our solution is simple in the sense that each discretization scheme allows (maintains) data protection or enhances the willingness to respond, while the properties of the underlying distribution are recovered by the combination of the resulting discretized variable. Building on our result of convergence in the distribution, we investigate three types of regression: i) the sensitive variable is an explanatory variable; ii) the sensitive variable is the outcome variable; and iii) sensitive variables are on both sides of the regression model. We use our split sampling method to point identify parameters in such linear models and derive the properties of the OLS estimator. The solution to identifying parameters in such context, is to condition the linear regression model on the appropriate known discretization intervals. The definition of the intervals, varies based on where the sensitive variable(s) are in the regression. Then split sampling allows to approximate the unknown distribution of this variable and obtain the conditional expectations needed to identify the parameter of interest without revealing their actual values. We show that this general method yields a consistent and asymptotically normal estimator for the parameters of interest in linear models. We then extend our discussion to nonlinear models and panel data as well. To demonstrate our method's finite sample properties, we provide some Monte Carlo evidence and apply our procedure to wage data from the Australian Tax Office (ATO) to investigate the gender wage gap in Australia. Our method complements the large literature on differential privacy (DP) to protect data privacy. DP has received much attention in the last decade and has become the industry standard.\footnote{It has been used by many technology companies, including Google erlingsson2014rappor, Microsoft ding2017collecting, LinkedIn kenthapadi2018pripearl, or Facebook. The United States Census Bureau employed DP in 2020 to safeguard individual confidentiality in the U.S. decennial census.} DP quantifies the notion of privacy for downstream tasks and aims to protect the most extreme observations as well. (see, e.g., dwork2006differentialprivacy, jordan2015machine, or bi2023distribution). There are two large branches of differential privacy; the first is “standard” or “central” DP, where data owners can only publish randomized statistics. Our paper relates to the second branch, which is called “local” DP, where data owners do not trust the “central server” or are not allowed/want to share the data in its original form. Local DP has multiple variants\footnote{dwork2014algorithmic uses a Laplace mechanism, rohde2020geometrizing uses local differential privacy, or avella2021privacy proposes parametric estimation, just to list a few variants.} and the major issue is that such data privatization may alter the analysis (estimation) results vis-a-vis the original data. As this problem is not new, Duchi02012018 and cai2021costofdp work out this trade-off between privacy and statistical accuracy in different setups, such as the mean estimator or least-squares parameters in linear regression. In the econometric context, bi2023distribution proposes a DP discretization that maintains statistical accuracy for a pre-specified model, while allowing for data privacy. Our approach to some extent is close to theirs in the sense of maintaining (asymptotic) consistency; however, this is the only link it shares with the differential privacy literature. For example, instead of adding noise (Laplace or any other) to the data, we employ discretization to “distort” the information.\footnote{Compared to the novel method introduced by bi2023distribution, in which they require specifying the model before the analysis, our discretization process is model-agnostic in the sense that it can be used with different model specifications. Furthermore, bi2023distribution achieves asymptotic convergence through splitting the sample into two parts and providing only one part to the end user. The other part is then used to transform the original data, add noise, and convert back to a representative distribution of the original data, resulting in an efficiency loss in the number of observations. Our method does not require such sample splitting.} Our approach also relates to the literature on partially identified parameters (e.g., see manski2003partial, manski2002inference, tamer2010partial, molinari2020microeconometrics, abrevaya2020interval_fe, pacini2019interval, and wang2022partial). This literature starts from the fact that the variable of interest is discretized (also called interval censored) for some reason (e.g., pacini2019interval discuss income surveys). manski2002inference show that the conditional expectation function and the regression parameters in such context cannot be point identified in general. They show how to derive set identification when the discretized variable is on the right-hand side. beresteanu2008asymptotic, bontemps2012set discuss the case when the outcome variable is discretized, while beresteanu2011discretized cover the case when both outcome and covariate(s) are discretized. The main toolkit of partially identified parameters in such context is moment (in)equalities.\footnote {Moment (in)equality models generalize the problems to multiple equations and/or inequalities, see for examples, chernozhukov2007estimation or andrews2010inference. beresteanu2008asymptotic show the asymptotic properties of such partially identified parameters. imbens2004confidence, chernozhukov2007estimation, and kaido2019confidence, among others, derive confidence intervals for these set-identified parameters. These methods are feasible ways to estimate parameter sets for a given conditional expectation function without any further assumptions.} The main advantage of these methods is that they do not require adding any noise to the variables or further distributional assumptions and allow valid statistical inferences on the parameters of interest. The magnitudes and signs of the estimated parameter vector can be interpreted in the same way as the classical regression coefficients. However, the main drawback is that the estimated parameters are not point identified but rather up to a set to which the parameter vector belongs. In many empirical applications, this set may be too wide to be economically meaningful. Our method circumvents the problem of using interval censored variables and shows that if a multiple discretization scheme is applicable, this additional information can be used to get point identification along with more precise estimates. The paper is organized as follows: Section (ref) introduces the discretization problem, data privacy, the parameters of interest, and also, the notation. Section (ref) describes the proposed split sampling approach. We discuss the shifting method in detail and show convergence in distribution to the unknown underlying distribution. In Section (ref), we identify parameters in a multivariate linear regressions and derive the asymptotic properties of the least square estimator for the cases of discretized variables on the right or left hand side. We also outline the case of discretization on both side of a regression model. Section (ref) discusses how to extend our method to nonlinear models and infer implications with panel data and fixed-effect type of estimators. Section (ref) presents some Monte Carlo evidence along with an empirical application on the Australian gender wage gap. Section (ref) concludes.
Consider $Z \sim f_Z(z;a_l, a_u)$ an i.i.d. random variable, where $f_Z(z;a_l,a_u)$ denotes the probability density function (pdf) with support $[a_l,a_u]$, where $a_l, a_u \in \mathbb{R}, \, a_l < a_u$. We assume that $f_Z(\cdot)$ is unknown and continuous. Let $Z_i$ be $i=1,\ldots, N$, realizations of $Z$. Instead of providing the sensitive variable $Z_i$, one observes a discretized version $Z^*_i$, through the following discretization process:
where $v_m \in \mathcal{C}_m$, $m=1 , \dots , M$ is the assigned value for each interval, and $\mathcal{M}^0$ is the discretization mechanism. The value $v_m$ can be a measure of centrality (e.g., midpoint) or an arbitrarily assigned value within its interval. $M$ denotes the number of intervals that are known. Also, the discretization intervals $\mathcal{C}_m$, are set by the data provider or the survey maker, and we take them as given. It is helpful to see the discretization intervals as independent from the random variable $Z$, but this is not a necessary requirement in our case. For simplicity, we use the terms interval and class interchangeably for $\mathcal{C}_m$.
In differential privacy, there is a well-established notion for data protection called $\varepsilon$-differential privacy (dwork2006differentialprivacy, dwork2006calibrating). The idea is the following: one changes a single observation in a dataset and wants to see how it impacts the publicly provided information. $\varepsilon$-differential privacy captures this notion as a ratio of probabilities. Let $\mathcal{Z}$ be an arbitrary sample taken from $f_Z(z;a_l,a_u)$, and $\mathcal{M}(\cdot)$ the privatization mechanism, which privatizes/discretizes $\mathcal{Z}$ into $\mathcal{Z}^*$ for public release. Consider two neighboring realizations of $\mathcal{Z}$; $\mathrm{z}$, and $\mathrm{z}'$ that are different in just one observation. In our case it means to drop one observation from $\mathrm{z}$ to get $\mathrm{z}'$.\footnote{In differential privacy it may also mean to add different noise to one observation.} Consequently, $\varepsilon$-differential privacy requires that the ratio of the probability of any privatized value given one sample to the probability of it given the other sample is upper bound by $e^\varepsilon$, that is:
where $\varepsilon \geq 0$ is a privacy factor. Note that the $\Pr[\cdot]$ refers to the randomness implied in the assignment mechanism $\mathcal{M}(\cdot)$. If $\varepsilon$ is small, we say privacy protection is high. Typically, the literature set $\varepsilon$ to 1 or 2. When the numerator and denominator are zero, the convention is to define them as 0. In the case when the denominator is zero, but the numerator is not, we face information leakage. To address this case, the literature uses the “approximate differential privacy”, given by:
where $\delta$ gives the probability of information leakage, typically a small value. Remark 1: The use of known fixed discretization intervals $\mathcal{C}_m$ allows to hide identity. Remark 2: In our setup, the discretization intervals' size $(M)$ mitigates the trade-off between privacy and accuracy. To illustrate this, the size of the discretization interval is $||\mathcal{C}_m|| = c_{m} - c_{m-1}$. One end of the spectrum is an infinite number of interval; $M\to \infty \implies ||\mathcal{C}_m|| \to 0$, so we observe each sensitive value in its original form without providing any privacy. On the other hand, if we use only one interval, $M=1 \implies ||\mathcal{C}_m|| = a_u - a_l$, which means $Z^*_i$ takes only one value. This case implies the numerator and the denominator are always the same, thus $\varepsilon = 0$. Remark 3: In our approach $\mathcal{C}_m$ are set and known, thus for a given discretization mechanism $\mathcal{M}(\cdot)$, one can calculate $\varepsilon$ and $\delta$ using Equation ((ref)).
We are interested in the regression parameter $\beta$, where $g(\cdot)$ is a known continuous function,
For the sake of simplicity, we introduce our notation in scalar terms; $\beta$ stands for a scalar parameter belonging to a subset of a compact finite-dimensional space $(\mathcal{B} \subset \R)$, while $Y, X$ are scalars. Facing the discretization of a variable can affect the identification of $\beta$ in three ways. In the first and simplest case, the explanatory variable $(X^*)$ is discretized; the second possibility is when the outcome variable $(Y^*)$ is discretized; and lastly, where both variables are discretized $(Y^*, X^*)$. To show an example for the identification problem, let us consider the following model:
Now, instead of observing $X$, we observe the discretized version of $X^*$, resulting in,
Under the usual assumptions, the naive OLS estimator for $\beta$ results in $\hat{\beta}^n = \left ( X^{*\prime} X^* \right )^{-1} X^{*\prime} Y = \left( X^{*\prime} X^* \right )^{-1} X^{*\prime} X \beta + \left( X^{*\prime} X^* \right )^{-1} X^{*\prime} u\,$, that is in general $\hat{\beta}^n \not\to \beta$ as $X^{*} \neq X$. Furthermore, the effect of the discretization is contagious; thus, if one uses further variable(s) $W$ with parameter(s) $\gamma$, that are correlated with $X$, then $\hat\gamma^n \not\to \gamma$. This result implies that the practice of assigning arbitrary values (e.g., midpoints) to the discretized variable leads to biases in the parameter estimates even when the discretized variables are controls. Note that the size of the bias depends on the underlying distribution of $X$ and the discretization process.\footnote{We provide more detailed discussion on this topic with proofs in \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 5.} In fact, this phenomenon is known in the literature, and manski2002inference shows that with discretized variables, only set identification is possible. However, in this paper, instead of applying set identification, we show that if the marginal distribution of the discretized variable is known, then under some mild conditions, we can point identify $\beta$ while preserving the discretization process.
The key to our approach is the use of split sampling, which recovers the marginal distribution of $Z$ even if it is discretized. Split sampling consists of two steps: 1) discretize the sensitive variable into multiple samples; 2) combine the samples by using the information on the known interval boundaries. Although the idea is simple, it requires some heavy notation for rigorous proofs. Thus, first, let us take a very simple illustrative example, and then we discuss the univariate case in general. As a last step, we extend our approach to the multivariate models. Let $Z$ be a random variable with the following pdf, with support boundaries $a_l = 0$ and $a_u=4$,
As the first step, let us discretize $Z$ with two intervals $M=2$ and let us use two discretization schemes or split samples. The first split sample discretizes the $Z$ as $[0, 2)$ and $[2, 4]$, while in the second split sample $[0, 1)$ and $[1, 4]$. The top section of Figure (ref) referred to as “split samples” visualizes this setup. By taking these two separate sets of information, we cannot recover the true probabilities for $Z$.
However, in the second step, we combine information from the two split samples into the “working sample”, which is represented in the bottom part of Figure (ref). In this simple example, the working sample has $3$ intervals, with boundaries $[0, 1)$, $[1, 2)$, and $[2, 4]$. By combining the information from the split samples to the working sample, we can deduce the probabilities of the non-observed intervals. Now, let us generalize our method. First, create multiple samples ($s=1,\dots,S$) and discretize them through different schemes while fixing the number of intervals ($M$) within each split sample. Discretization for the $s$-th split sample looks exactly as the problem introduced in Equation ((ref)), with the only difference being that the interval boundaries are different,\footnote{To simplify the notation, we use instead of $Z^{*(s)}$ simply $Z^{(s)}$.}
The number of observations across split samples $(N^{(s)})$ can be the same or, more likely, different. As a second step, we combine information into the “working sample”. The construction of the working sample's interval is simply given by the union of unique interval boundaries within each split sample,
One can generalize the approach with any type of distribution. Let us write the probability for a random observation to be in a given class of a split sample,
where $\mathcal{S}_s$ denotes the $s$-th split sample. We can express, the probability of an observation falling into the working sample's $\mathcal{C}^{WS}_b$ interval as
In practice, we recommend splitting the sample size into equal parts, thus $\Pr( Z \in \mathcal{S}_s ) = 1/S$ and use uniform probability for $\Pr \left( Z \in \mathcal{C}^{WS}_b \mid Z \in \mathcal{C}^{(s)}_m \right)$. Then Equation ((ref)) simplifies to
The shifting method is an application of split sampling. It uses an equal distance discretization scheme as the benchmark and shifts the boundaries of the intervals with a certain value. The interval widths for the different split samples remain the same, except for the boundary intervals. As one increases the number of split samples, the size of the shift becomes smaller and smaller. This enables the mapping of the marginal probability distribution at the limit. First, we introduce the discretization process proposed for the shifting method, and then, as a second step, we show how to use the realizations of this discretization to create a synthetic variable that maps the original variable's marginal distribution. Figure (ref) shows an illustrative example for split samples with $S=4$ and $M=4$ classes.
To derive the properties of the shifting method, we define the uniform class widths for the first split sample as $\frac{a_u-a_l}{M-1}$. We split these intervals into $S$ parts, thus we can shift the class boundaries $S$ times. This defines the size of the shift as
This implies that the number of intervals in the working sample is $B = S(M-1)$. The boundary points for each split sample are given by, \footnote{In the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement} under Section 1, Algorithm A1 summarizes how to construct split samples using the shifting method.}
As the second step, we introduce a synthetic random variable $Z^\dagger$ that maps the underlying marginal distribution for $Z$. Observations from a particular interval in the split sample $s$ can end up in several candidate intervals from the working sample $\mathcal{C}^{WS}_b$. To make it more transparent, we define the discretization intervals $\mathcal{C}^{(s)}_b$ as sets of $\mathcal{C}^{WS}_b$,
We create the synthetic variable $Z^\dagger$ by assigning each discretized observation from its split sample $\mathcal{C}^{(s)}_m$ to one of its corresponding working sample's intervals $\mathcal{C}^{WS}_b$ with uniform probability. Note that this assignment mechanism does not assume any distributional form for $Z$ itself, but it randomly assigns each value of $Z^{(s)}$ to the working sample's related interval with uniform probability to get $Z^\dagger$. The random assignment mechanism can be written as
The unconditional probability of $Z^\dagger \in \mathcal{C}^{WS}_b$ for the shifting method, following Equation ((ref)) and assuming ${\Pr( Z \in \mathcal{S}_s )=1/S}$ now is
Next, we outline the assumptions under which $Z^\dagger \overset{d}{\rightarrow} Z$ as the number of split samples $(S)$ and observations $(N)$ go to infinity.
Assumption (ref).a) ensures that the number of observations is always higher than the number of split samples.\footnote{Note, that we can decrease $\mathrm{c}$ to be arbitrarily close to $0$.} This is required to identify each point of $f_Z(z;a_l,a_u)$ through the mapping of split samples to the working set. Assumption (ref).b) imposes a mild assumption on the underlying distribution: the support of the random variable is not disjoint within the working samples' interval, thus $\int_{c^{WS}_{b-1}}^{c^{WS}_b} f_Z(z)\mathrm dz > 0 \,,\,\, \forall c^{WS}_b$.
Let us generalize our results to the multivariate case. Let $\bfZ = (\bfz_1,\dots,\bfz_p,$ $\dots,\bfz_P)$ distributed as $\bfZ \sim \bff_{\bfZ}(\bfz;\bfa_l,\bfa_u)$, i.i.d. random variables with $P$ dimensions, where $\bff_{\bfZ}(\cdot)$ denotes the (joint) pdf with $\bfz = (z_1,\dots,z_p,\dots,z_P)$,\footnote{Note here that $\bfz$ is not a row from $\bfZ$, but represents a vector that is used with the distribution $\bff_{\bfZ}(\cdot)$} and known support $\mathbb{A}^P = [a_{1,l}, a_{1,u}] \times \ldots \times [a_{P,l}, a_{P,u}] \subseteq \R^P$. We assume $\bff_{\bfZ}(\cdot)$ is unknown and continuous. Instead of observing $\bfZ$, we observe $\bfZ^*$ through a (multivariate) discretization process similar to Equation ((ref)). The difference is that instead of intervals, we use a multi-dimensional grid,\footnote{One-by-one discretization does not work in our case as it will not be able to recover the joint distribution of $\bff_{\bfZ}(\cdot)$. See more in \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 1.3.} where each grid is defined as $\bf\mathcal{C}_m = [c_{1,m-1}, c_{1,m}) \times \ldots \times [c_{P,m-1}, c_{P,m})$ with $c_{m,p}$ being the $m$'th grid point for variable $p$. Note that $c_{p,0} = a_{p,l}$ and $c_{p,M} = a_{p,u}$, for $p=1,\ldots, P$, thus the first boundary vector and last boundary vector contains the boundaries of the distribution's support. $\bfv_m \in \bf\mathcal{C}_m$ is the assigned vector for each grid, that can be a measure of centrality (e.g., focus point) or other point within the grid. In the case of $P=2$, we have a 2D grid that defines the discretization scheme, and each $\bf\mathcal{C}_m$ stands for a rectangle with $M^2$ different partitions. The shifting method in the multivariate case utilizes the initial discretization scheme with $\bf\mathcal{C}_m$ and shifts the boundary points of the grids with $h_p = \frac{a_{p,u} - a_{p,l}}{S(M-1)}$ for each $p=1,\ldots, P$. This results in $(S(M-1))^P$ different grids in the working sample, which is much larger than the initial $M^P$. Changing the grids' boundary points allows for mapping the underlying distribution using the same inductive logic to prove convergence in distribution. More formally, first consider the following assumptions:
See the complete proof in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement} under Section 1.1.
We investigated the speed of convergence for $F_{Z^\dagger}(\cdot) \to F_Z(\cdot)$ as a function of the number of split samples $(S)$ to get insight into the choice of $S$ in finite samples for the univariate case. In general, the rate of convergence is $1/S$, except around the boundaries of the support. Boundaries of the support are defined as: $c^{(1)}_1 + (c^{(S)}_2-c^{(1)}_1)$ for the lower bound, and $c^{(1)}_M + (c^{(1)}_{M-1}-c^{(1)}_M)$ for the upper bound. In these regions, the speed of convergence is $\frac{\log S}{S}$. This result means that even with small values of $S$, one can get a decent mapping of the underlying marginal distribution of $Z$.\footnote{We show the detailed derivations for the speed of convergence in \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, under Section 1.2.}
Let us consider the case when $Z$ has unbounded support: $a_l = -\infty$ or $a_u = \infty$ or both. If we set the first and last intervals' boundary points to infinity, we are facing censoring. As will be discussed in Section (ref), one can extend our approach for censoring by using the continuous mapping theorem. However, here we propose a simpler solution: truncate observations that fall into the first or last interval: $\mathcal{C}^{(s)}_1, \mathcal{C}^{(s)}_M, \,\, \forall s$. Although this method loses observations and important units, the flip side is that for the truncated distribution now all previous results apply.
The $\varepsilon$-differential privacy can also be used for split sampling methods. The shifting method introduces a random assignment of each unit to a specific split sample (with a uniform probability of $1/S$). The user receives both the discretized value and the discretization scheme.\footnote{Creating the synthetic variable $Z^\dagger$ can be done by the user; no additional information is gained or provided from that step.} Realized interval sizes are specific to the shifting method. In general, $||\mathcal{C}^{(s)}_m|| = \frac{a_u-a_l}{M-1},$ $\forall m = 2,\dots, M-1$, thus it does not depend on the number of split samples, only on the number of intervals $M$. However, at the lower and higher bounds, the interval sizes are different: $||\mathcal{C}^{(s)}_m|| = (s-1) \frac{a_u-a_l}{S(M-1)},\,\, \forall s,\, \text{ with } m=1 \text{ or } m=M$. This implies that around the boundaries, the number of split samples influences the size of the intervals, thus $\varepsilon$ depends on both $M$ and $S$, that causes an overall weaker privacy protection. As we show in our empirical application, this does not have real empirical relevance. Note, however, that truncation of the boundary values could provide a simple fix for this issue.
It is possible to work out other discretization schemes that can map the underlying marginal distribution of $Z$. In fact, in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, under Section 4, we introduce the magnifying method, which offers an alternative method to learn $f_Z(z;a_l,a_u)$. The key to all possible split sampling methods is to show convergence in distribution. Different discretizations may yield different speeds of convergence and other $\varepsilon$-differential privacy properties.
This section discusses how to use the synthetic $Z^\dagger$ variable to identify and estimate $\beta$ in a linear regression model. As linear models and OLS have well-known properties, we directly address modelling with multiple regressors. First, we discuss the identification $\beta$ for each model, depending on where the discretization happens. This step helps us to pin down which theoretical expressions to quantify. In all cases, we highlight the discretized variable with an asterisk in the superscript. Let $\bfy = (Y_1,\ldots,Y_N)'$ a $(N\times 1)$ vector, $\bfX = (\bfx_1,\dots,\bfx_k,\dots,\bfx_K)$ with dimensions $(N \times K)$ where $\bfx_k = (\bfx_{k,1},\ldots,\bfx_{k,N})'$. The $(K \times 1)$ parameter vector of interest is $\bfbeta = (\beta_1,\ldots,\beta_K)'$. We allow further conditioning variables that are not discretized denote by $\bfW = (\bfw_1,\dots,\bfw_j,\dots,\bfw_J)$ i.e., a $(N \times J)$ matrix with $\bfw_j = (\bfw_{j,1},\ldots,\bfx_{j,N})'$ and parameter vector $\bfgamma = (\gamma_1,\ldots,\gamma_J)'$, as this mimics standard practices better. $\bfbeta$ and $\bfgamma$ belong to a subset of a compact finite-dimensional space $(\mathcal{B},\mathcal{G})$. Let the unknown data-generating process (DGP) be
where $\bfX$ and $\bfW$ are linearly separable. $\bfu$ is homoskedastic (for simplicity), and we use the usual OLS assumptions: $\mathbb{E}[\bfu|\bfX,\bfW] = 0$ and $\operatorname*{plim}_{N\to\infty}\left[\frac{1}{N}\bfX\bfM_{\bfW}\bfX\right]^{-1}=\bfQ$, positive definite matrix, with the usual residual maker $\bfM_{\bfW} = \bfI - \bfW (\bfW'\bfW)^{-1}\bfW'$. In all three discretization cases, we follow a three-step procedure:
First, let us discuss the case when we observe $\bfX^*$, instead of $\bfX$. We utilize the result of Proposition (ref), hence $\bfX^\dagger \overset{d}{\rightarrow} \bfX$, where $\bfX^\dagger = \left(\bfx^\dagger_1,\dots,\bfx^\dagger_k,\dots,\bfx^\dagger_K\right)$ and $\bfx^\dagger_k = (\bfx_{k,1},\ldots,\bfx_{k,N})'$.
To identify $\bfbeta$, from the DGP defined by Equation ((ref)), we need to use the conditional expectations given $\bfX^* \in\bf\mathcal{C}^{(s)}_m$. Applying the conditional expectation operator on our model yields,
Note that the classes, $\bf\mathcal{C}^{(s)}_m$ are mutually exclusive for given $s$ along $m$ and $k$. Let us define $\kappa(s,\bfm)$ a function\footnote{Note that $\bfm = (m_1,\dots,m_K)$ is a vector, defining the grid points in $\bf\mathcal{C}^{(s)}_m$ for each variable $k$.} that maps the conditional expectations of $\bfX$ for given grid $\bf\mathcal{C}^{(s)}_m$. Proposition (ref) ensures equality between the conditional expectation of the observable $\bfX^\dagger$, and the conditional expectation of the unknown $\bfX$ given $\bfX^* \in \bf\mathcal{C}^{(s)}_m$, which is needed to identify $\bfbeta$,
Let $\bfkappa = (\kappa(1,{\bf1}),\dots,\kappa(S,{\bf1}),\dots,\kappa(S,\bfm),\dots,\kappa(S,\bfM))'$ be a vectorized version of $\kappa(s,\bfm)$ with $(SM^K \times 1)$ dimension. We propose a joint estimation for $\bfkappa$ using OLS\footnote{We can estimate all conditional expectations jointly via OLS. Otherwise simple conditional means for a given interval $\bf\mathcal{C}^{(s)}_m$ would yield the same result. Joint estimation however allows to establish asymptotic properties more easily.},
where $\mathbf{1}_{\{\cdot\}}$ stands for the indicator function and results in a matrix with dimensions $(N\times SM^K)$ with the same ordering as $\bfkappa$. Let us note, $\hat{\bfkappa}$ does not require the actual values of $\bfX$, only $\bfX^\dagger$ and $\bfX^\dagger \in \bf\mathcal{C}^{(s)}_m$. Under standard OLS assumptions, $\sqrt{N}\left(\hat{\bfkappa} - \bfkappa \right) \mathrel{\overset{\makebox[0pt]{\mbox{\scriptsize a}}}{\sim}}{} \mathcal{N}\left(\mathbf{0},\bfOmega_{\bfkappa} \right)$, where $\bfOmega_{\bfkappa}$ is the variance-covariance matrix. See the derivations for the proposed estimator in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, under Section 2.1.
To get a consistent estimator for $\bfbeta$, let us follow Equation ((ref)), and define $\ddot\bfX$, which replaces $\bfX^*$ with the corresponding conditional expectations via $\hat\bfkappa$. Furthermore, let $\ddot\bfy = \hat{\mathbb{E}}\left[\bfy|\bfX^* \in \bf\mathcal{C}_m\right]$ and $\ddot\bfW = \hat{\mathbb{E}}\left[\bfW|\bfX^* \in \bf\mathcal{C}_m\right]$ the respective conditional averages.
while $\bfe = (e_1,\dots,e_N)$. Now the OLS estimator for $\bfbeta$ is
where $\bfM_{\ddot\bfW}$ is the usual residual maker using $\ddot\bfW$. Note that $\E[e_i]=0$ for all $i$ since $\ddot{\bfX}$ is a consistent estimate. Moreover, $\E[e_i e_j] = 0$ for $i\neq j$ due to $\bf\mathcal{C}^{(s)}_m$ being mutually exclusive in $m$ and $k$. As long as the discretization is mean independent from the error term: $\mathbb{E}\left[\bfu|\bfX,\bfW,\bfX^*\in\bf\mathcal{C}_m\right] = \mathbb{E}\left[\bfu|\bfX,\bfW\right] = 0$, the OLS is a consistent estimator.
The asymptotic distribution of $\hat\bfbeta$ can be derived similarly. As $\bfu \sim iid \left ( \bfzero, \sigma^2_{\bfu} \bfI \right )$ and as $\ddot{\bfX}$ and $\ddot{\bfW}$ are consistent estimates of the conditional means, $\mathbb{E}\left[\bfe\bfe'|\ddot{\bfX},\ddot{\bfW}\right] = \sigma_{e}^2\bfI$. Under assumptions $N,S \to \infty$ and $\mathbb{E}\left[\bfu|\bfX,\bfW,\bfX^*\in\bf\mathcal{C}_m\right] = \mathbb{E}\left[\bfu|\bfX,\bfW\right]$, we get $\sqrt{N}\left(\hat\bfbeta-\bfbeta\right) \mathrel{\overset{\makebox[0pt]{\mbox{\scriptsize a}}}{\sim}}{} \mathcal{N}\left(\mathbf{0}, \sigma_{e}^2\left(\ddot{\bfX}'\bfM_{\ddot\bfW}\ddot{\bfX}\right)^{-1} \right)$.
The second case is when the outcome variable is discretized, thus $\bfy$ is not observed, only $\bfy^*$. To simplify our analysis let us neglect $\bfW$ along with the parameter vector $\bfgamma$ as they can be seen as part of $\bfX$ and $\bfbeta$ respectively.
To identify $\bfbeta$ with discretized $\bfy^*$, let us partition the domain of $\bfX$'s into mutually exclusive partitions, denoted by $\bf\mathcal{D}_l$ and defined by the researcher, such that it is (mean) independent from the error term $\bfu$. We use the same logic as derived in Equation ((ref)), but instead of conditioning on $\bfX^* \in \mathcal{C}_m$, we use $\bfX \in \bf\mathcal{D}_l$. Similar derivation can be done,
$\Pr[\bfy^* \in \mathcal{C}_m | \bfX \in \bf\mathcal{D}_l], \mathbb{E} \left[ \bfX | \bfX \in \bf\mathcal{D}_l \right]$ are known quantities and we aim to learn $\mathbb{E}[\bfy | \bfy^* \in \mathcal{C}_m, \bfX \in \bf\mathcal{D}_l]$. Note that $\mathbb{E} \left[ \bfu | \bfX \in \bf\mathcal{D}_l \right] = 0$ by the definition of $\bf\mathcal{D}_l$. Let us define $\pi(s,m,\bfl)$ a function that maps the conditional expectations of $\bfy$ for a given interval $\bfy \in \mathcal{C}^{(s)}_m$ and partition $\bfX \in \bf\mathcal{D}_l$.
where $\bfy^\dagger = (Y^\dagger_1,\ldots,Y^\dagger_N)$ and we used the results of the shifting method to show convergence.
Let $\bfpi = \left(\pi(1,1,{\bf1}),\dots,\pi(S,1,{\bf1}),\dots,\pi(S,M,{\bf1}),\dots,\pi(S,M,\bfl),\dots,\pi(S,M,\bfL)\right)'$ be the vectorized version of $\pi(s,m,\bfl)$. Estimating $\bfpi$, via OLS yields,
The vector $\hat{\bfpi}$ has dimensions of $(SML^K \times 1)$, $ \mathbf{1}_{\{\cdot\}}$ is the corresponding indicator matrix, and estimation of $\bfpi$ does not require the actual values of $\bfy$, only $\bfy^\dagger$, and observing $\bfy^\dagger \in \mathcal{C}^{(s)}_m$ and $\bfX \in \bf\mathcal{D}_l$. Under the standard OLS assumption, $\sqrt{N}\left(\hat{\bfpi} - \bfpi \right) \mathrel{\overset{\makebox[0pt]{\mbox{\scriptsize a}}}{\sim}}{} \mathcal{N}\left(\mathbf{0},\bfOmega_{\bfpi} \right)$, where $\bfOmega_{\bfpi}$ is the corresponding variance-covariance matrix. See the derivations for the proposed estimator in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, under Section 2.2.
As we have shown in our identification strategy, the standard regression model can be rewritten in the form of Equation ((ref)). Let $\tilde{\bfy}$ be a $(N\times 1)$ vector, where we replace the discretized elements of $\bfy^*$ with the consistent estimates of $\sum_l\sum_m\mathbb{E}[\bfy | \bfy^* \in \mathcal{C}_m, \bfX \in \bf\mathcal{D}_l] \Pr[\bfy^* \in \mathcal{C}_m | \bfX \in \bf\mathcal{D}_l]$, using $\hat\bfpi$ and sample analogues for $\Pr[\bfy^* \in \mathcal{C}_m | \bfX \in \bf\mathcal{D}_l]$. Let us redefine $\tilde{\bfX} = \hat\mathbb{E} \left[ \bfX | \bfX\in \bf\mathcal{D}_l \right]$, as the conditional sample means for $\bfX$ in partition $\bf\mathcal{D}_l$. With the replaced measures, we get
where $\bfnu = \left ( \nu_1, \ldots, \nu_N \right )'$. Note that $\E[\nu_i]=0$ for all $i$ since $\tilde{\bfy}$ and $\tilde{\bfX}$ are consistent estimates. Moreover, $\E[\nu_i \nu_j] = 0$ for $i\neq j$ due to that $\bf\mathcal{D}_l$ are mutually exclusive $\forall l,$ and $\E [\nu_i | \tilde{\bfX} ] = 0$ since the partitioning is (mean) independent from the error term and does not affect the sampling error. Let $\hat{\bfbeta} = \left ( \tilde{\bfX}' \tilde{\bfX} \right )^{-1} \tilde{\bfX}' \tilde{\bfy}$, then under similar argument, $\hat{\bfbeta} - \bfbeta = o_p(1)$. The asymptotic distribution of $\hat\bfbeta$ can be derived in the same spirit. As $\bfu \sim iid \left( \bfzero, \sigma^2_{\bfu} \bfI \right )$ and $\tilde{\bfX}$ is a consistent estimate of the conditional means, $\mathbb{E}\left[\bfnu\bfnu'|\tilde{\bfX}\right] = \sigma_{\nu}^2\bfI$. The proposed estimator is asymptotically distributed as $\sqrt{N}\left(\hat\bfbeta-\bfbeta\right) \mathrel{\overset{\makebox[0pt]{\mbox{\scriptsize a}}}{\sim}}{} \mathcal{N}\left(\mathbf{0}, \sigma_{\nu}^2\left(\tilde{\bfX}'\tilde{\bfX}\right)^{-1} \right)$.
Our last case is when the discretization happens with both outcome $\bfy^*$, and with one or more explanatory variables $\bfX^*$. In this case, we do not need to partition the domain of $\bfX$, but can use the discretization grids $\bf\mathcal{C}^{(s)}_m$ for the explanatory variable. As $\bfy$ is also discretized, we need to partition $\bfW$ with $\bf\mathcal{D}_l$, similarly to the case discussed in Section (ref). Identification is the same as with discretized outcome variable and the sample estimator for the conditional expectations follows the same logic as the OLS estimator for the $\bfbeta$ parameter. We discuss this case in details in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, under Section 3.
A possible extension is the use of nonlinear models with one or more variables discretized. Let us consider the following general model:
where $g(\cdot)$ denotes a known continuous function. The identification procedure is similar to Equations ((ref), (ref)), with the exception that $g(\mathbb{E}[\cdot;\bfbeta]) \neq \bfbeta g(\mathbb{E}[\cdot])$ in general. As shown with split sampling $\bfZ^\dagger \overset{p}{\rightarrow} \bfZ$. If so, with discretized explanatory variable(s) by continuous mapping theorem $g(\bfX^\dagger;\bfbeta)~\overset{p} {\rightarrow}~g(\bfX;\bfbeta)$. Similar argument applies when discretization happens on the left hand side $(\bfy^*)$ or on both sides $(\bfy^*,\bfX^*)$. As the estimation, let $\hat{\bfbeta}(\bfy, \bfX,\bfW)$ denote a consistent estimator of $\bfbeta$ with $\rho (\bfX,\bfW) =$\\ $\sqrt{N} \left [ \hat{\bfbeta} (\bfy, \bfX, \bfW) - \bfbeta \right ]$ such that $\rho(\bfy, \bfX,\bfW) \overset{d}{\rightarrow} \bff(\bf0, \bfOmega)$. Under the assumptions made earlier, for all cases, the nonlinear quantities should converge too. When the regressors are discretized, we need $\bfX^\dagger \overset{d}{\rightarrow} \bfX \implies \rho (\bfy, \bfX^\dagger,\bfW) \overset{d}{\rightarrow} \rho (\bfy, \bfX,\bfW) $. In case of discretized outcome, $\bfy^\dagger \overset{d}{\rightarrow} \bfy \implies \rho (\bfy^\dagger, \bfX,\bfW) \overset{d}{\rightarrow} \rho (\bfy, \bfX,\bfW)$; while for both variables $(\bfy^\dagger, \bfX^\dagger) \overset{d}{\rightarrow} (\bfy,\bfX) \implies \rho (\bfy^\dagger, \bfX^\dagger,\bfW) \overset{d}{\rightarrow} \rho (\bfy, \bfX,\bfW)$. By the continuous mapping theorem under appropriate regularity conditions these convergences hold. The technical details of these conditions, however, could be an interesting subject of future research. In case of panel data, the extension of this methodology is relatively straightforward. The most important difference is in the identification. If the interval of an individual does not change over the time periods covered, the individual effects in the panel and the parameter associated with the choice variable cannot be identified separately. The within transformation would wipe out the interval variable as well. When the interval does change over time, but not much, then we are facing weak identification, i.e., in fact very little information is available for identification, so the parameter estimates are going to be highly unreliable.\footnote{See more on this issue in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 5.6} In the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 6, we extend our method towards fixed effect type of estimators. This is a good opportunity to raise awareness through a phenomenon that we call the {\it perception effect}. The perception effect is relevant if the discretization happens by surveys. There is much evidence in the behavioral literature that the answers to a question may depend on the way the question is asked (see, e.g., Haisley2008). Note, that this is present regardless of whether split sampling has been performed or not. However, with split sampling, there is a way to tackle this issue, much akin to the approach a similar problem has been dealt with in the panel data literature.
For the simulation experiments, we consider a simple univariate linear regression model
where we discretize $Y_i$, or $X_i$ or both. We set the parameter of interest $\beta=0.5$ and report the Monte Carlo average bias and its standard deviation in parenthesis with $1,000$ repetitions and $N=10,000$ observations. As the size of the bias depends on mapping the conditional expectations, we use multiple distributions for $\epsilon_i$ or $X_i$. To be more specific, we use “Normal” (standard normal distribution truncated at $-1$ and $3$), “Logistic” (standard logistic distribution truncated at $-1$ and $3$), “Log-Normal” (standard log-normal distribution truncated at $4$ and subtracted $1$ to adjust the support), “Uniform” (uniform distribution between $-1$ and $3$), “Exponential” (exponential distribution with rate parameter $0.5$, truncated at $4$ and subtracted $1$) and “Weibull” (weibull distribution with shape parameter $1.5$ and scale parameter $1$, truncated at $4$ and subtracted $1$). All of these distributions have known finite support between $-1$ and $3$. We use these distributions to generate $\epsilon_i$ when the outcome or both variables are discretized. In these cases $X_i \sim \mathcal{N}(0,0.25,-1,1)$ (truncated normal distribution between -1 and 1). When discretization happens with the explanatory variable $X_i$, we use the different distributions to generate $X_i$ and use the same truncated normal distribution for $\epsilon_i$. This setup allows us to control for the domain of $Y_i$ in all cases. For the discretization of the variable(s) in all cases we use $M=5$ and for split sampling $S=10$. We have experimented with different setups; see more in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 7. To estimate $\beta$, our method is shown by the “\textit{Shifting method}”.\footnote{We used mid-values as observations for the working sample's values $(v^{WS}_b)$.}. As alternatives, we use “\textit{Mid-point regression}”, a commonly used naive model. This method uses mid-points for $v_m$ with OLS for estimation. When the outcome variable is discretized, we use further popular estimation methods: “\textit{Set identification}” which provides estimates for the lower and upper bounds of the parameter set;\footnote{Estimation is based on beresteanu2008asymptotic, we use the package by beresteanu2010stata.} “\textit{Ordered probit and logit}” models\footnote{Note that the estimated maximum likelihood “naive” parameters reported here are not designed to recover $\beta$ and to be interpreted in the linear regression sense. Therefore, we call the difference from $\beta$ {\sl distortion} rather than bias.} and “\textit{Interval regression}”\footnote{We assume a Gaussian conditional distribution to model the censored interval outcome.} as alternative estimators. We specify our models when applicable to estimate as $Y_i = \alpha + X_i'\beta + \eta_i$, where $\epsilon_i = \alpha + \eta_i$ with $\mathbb{E} (\epsilon_i) = \alpha$, to adjust for those distributions, which does not have zero mean.\footnote{Ordered choice models' implementation in Stata removes the intercept parameter to identify $\beta$.}
Table (ref) shows our results. The shifting method provides smaller average bias than the mid-point regression by magnitudes between 5-50 almost everywhere. The two exceptions can be found in the case of discretized explanatory variable. When we use a uniform distribution, there is no bias for mid-point regression, as the conditional expectation is the mid-point of the interval. The second case is the exponential distribution, where the case is similar, as the curvature of the used distribution is rather flat.\footnote{This result is in line with the theoretical results shown in the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 5.} The shifting method also outperforms the other competing methods, when discretization happens on the left hand side and provides significant reduction in the bias when both variables are discretized. We have run several other Monte Carlo experiments and the results are similar: the shifting method outperforms all alternatives. An important result is consistency, as we increase $N$ along with $S$, we get smaller biases. This consistency does not hold for the competing methods. For a detailed discussion see the \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 7.
We illustrate our approach with a short study of the Australian gender wage gap from 2017. We obtained wage and socio-economic data from the Australian Tax Office (ATO). The individuals' sample files record a $2\%$ sample of the whole population,\footnote{For details, see ATO's website: \url{https://www.ato.gov.au}.} including the actual wage values. This enabled us to estimate $\beta$ when no discretization happens. In practice, researchers use interval censored data coming from surveys or privatized data. To mitigate this fact, we employ three different equally distanced discretization methods with $M=3,5$, and $10$ and a privatization method that adds Laplacian noise to the data.\footnote{We used the python package diffprivlib by diffprivlib for differential privacy implementation.} Table (ref) shows regression estimates of the gender dummy (1 is female) on log yearly wage while controlling for age, age square, marital status, occupation, and regional variables. The “Directly observed” row shows the estimated parameter when the actual values are used for yearly wages. In the following rows we show the estimates when the discretization process is done in equal-sized intervals with $M=3,5,10$. Mid-point regression does not use split samples, while with the shifting method we use $S=10$. In all cases, the yearly wage was discretized, and then we took the log of the discretized mid values as this represents the practice. We report the point estimates along with the standard errors for $\hat\beta$. Furthermore, we report the calculated $\varepsilon$-differential privacy measures based on Equation ((ref)). Finally, we report an estimate for $\hat\beta$ using differential privacy method, while setting $\varepsilon = 1$. We need to note that the estimates using differential privacy are highly volatile depending on the randomization.
The results show clearly that the shifting method provides closer estimates than mid-point regression everywhere, whereas the latter statistically produces different results from the “directly observed” value.\footnote{ For more discussion on the empirical example and for other model setups, see \href{https://regulyagoston.github.io/papers/MSV_online_supplement.pdf}{online supplement}, Section 8.} Observe that with the shifting method, $\varepsilon$-differential values are larger but still much lower than the commonly used value of 1 in the literature. The differential privacy method, while using a larger value for $\varepsilon$-differential than any of the realized $\varepsilon$-differential value for the discretization method (which means less data protection), produces a worse and imprecise measure. This is not surprising, as it is documented in the differential privacy literature that these methods comes with cost of accuracy (see, e.g., bowen2020comparative, cai2021costofdp, or bi2023distribution).
This paper deals with linear models using sensitive variables. We propose to use the discretization process to protect data privacy or increase response rates in surveys. This results in discretized (also called interval-censored) variables, which makes econometric modeling difficult as the conditional expectation cannot be point identified in general. We propose the split sampling method that introduces multiple discretization schemes; thus, instead of using one set of intervals, the sensitive variable is discretized through multiple versions. We combine these discretized realizations in a way that, under some mild conditions, they converge in distribution to the original (unknown) variable. We introduce the shifting method, an easy to implement algorithm for split sampling that preserves data privacy while also satisfying the conditions for convergence in distribution. With the help of the shifting method, we can point identify parameters of interest in linear models. The necessary point identification assumptions depend on where the discretization happens: i) the sensitive variable is an explanatory variable; ii) the sensitive variable is the outcome; or iii) discretization happens on both sides. We examine each case and derive the appropriate OLS estimators in a multivariate regression setup. We show the asymptotic properties of these estimators and discuss extensions to nonlinear models and panel data. We provide some Monte Carlo evidence to show that our methods have superior finite sample properties compared to the “usual” ones. Finally, we apply our method to estimate the Australian gender wage gap. We achieve not only smaller $\varepsilon$-differential values that show better data protection properties, but consistent parameter estimates even when the number of intervals is small.