Nonuniform Negative Sampling and Log Odds Correction with Rare Events Data
Year: 2021

Abstract:
We investigate the issue of parameter estimation with nonuniform negative sampling for imbalanced data. We first prove that, with imbalanced data, the available information about unknown parameters is only tied to the relatively small number of positive instances, which justifies the usage of negative sampling. However, if the negative instances are subsampled to the same level of the positive cases, there is information loss. To maintain more information, we derive the asymptotic distribution of a general inverse probability weighted (IPW) estimator and obtain the optimal sampling probability that minimizes its variance. To further improve the estimation efficiency over the IPW method, we propose a likelihood-based estimator by correcting log odds for the sampled data and prove that the improved estimator has the smallest asymptotic variance among a large class of estimators. It is also more robust to pilot misspecification. We validate our approach on simulated data as well as a real click-through rate dataset with more than 0.3 trillion instances, collected over a period of a month. Both theoretical and empirical results demonstrate the effectiveness of our method.

heading: 1 Introduction
text: Imbalanced data are ubiquitous in scientific fields and applications with binary response outputs, where the number of instances in the positive class is much smaller than that in the negative class. For example, in online search or recommendation systems where billions of impressions appear each day, non-click impressions usually dominate. Using these non-click impressions as negative instances is prevalent and proved to be useful especially in modern machine learning systems [21, 12, 4, 9, 14]. A common approach to address imbalanced data is to balance it by subsampling the negative class [6, 19] and/or up-sampling the positive class [2, 10, 20, 5], and a great deal of attentions have been drawn to this problem, see [15, 17, 3, 7, 23, 27, 1, 25, 8, 18, 30], and the references therein.
In this paper, we focus on subsampling the negative class, i.e., negative sampling, while keeping all rare positive instances. This is beneficial for modern online learning systems, where negative sampling significantly reduces data volume in distributed storage while saving training time for multiple downstream models. Rare events data and negative sampling are studied in [30], but it focused on linear logistic regression and only considered uniform sampling. We approach this problem under general binary response models from the perspective of optimal subsampling, which aims to minimize the asymptotic variance of the IPW estimator [32, 28, 31, 33]. This topic was not well investigated neither for imbalanced data nor for negative sampling. In addition, existing optimal
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
probability formulation minimizes the conditional asymptotic variance and the variations due to the data randomness are ignored. We fill these gaps by deriving the asymptotic distributions for both full data and subsample estimators under general binary response models, and we find the unconditional optimal probability for negative sampling. In addition, since IPW may inflate the variance [13, 22], we develop a more efficient likelihood-based estimator through nonuniform log odds correction to avoid IPW. Our main contributions are summarized as follows.
• Under a general binary response model with rare events data, we find that the difference between the full data estimator and the true parameter converges to a normal distribution at a rate that is tied to the number of rare positive instances only. This indicates the possibility to throw away the majority of negative instances without information loss, i.e, it justifies the usage of negative sampling.
• We show that there is no asymptotic information loss with aggressive negative sampling if the negative instances kept dominates the positive instances. However there is information loss expressed as a variance inflation if the negative instances are brought down to the same level of the positive instances. For this case, we obtain optimal subsampling probabilities by minimizing the unconditional asymptotic variance of a general IPW estimator.
• We develop a likelihood estimator through nonuniform log odds correction for sampled data, and prove that it has the smallest asymptotic variance among a large class of estimators.
• We apply the proposed method to a real online streaming click-through rate (CTR) dataset with more than 0.3 trillion instances and demonstrate its effectiveness.
The rest of paper is organized as follows. Section 2 defines the problem this paper focuses on and shows the full data results. Section 3 presents the optimal probability for the IPW estimator. Section 4 proposes the likelihood-based estimator and establishes its asymptotic optimality. Section 5 discusses practical implementation and the theoretical properties of the resulting estimators. Section 6 presents numerical experiments with simulated data, and applies the proposed method to a real online streaming CTR dataset with more than 0.3 trillion instances. Section 7 concludes the paper and points out limitations of our investigation. Proofs and additional experiments with misspecifications are in the supplementary materials.

heading: 2 Problem setup and full data results
text: Let {(xi, yi)}Ni=1 be training data of size N that satisfies the binary response model, Pr(y = 1 | x) = p(x;θ), (1)
where y ∈ {0, 1} is the binary class label (y = 1 for the positive and y = 0 for the negative), x is the vector of features, and θ is the unknown d-dimensional parameter.
Let N1 be the number of positive instances (N1 = ∑N i=1 yi) and N0 be the number of negative instances (N0 = N −N1). We consider the scenario of massive imbalanced data, i.e., N1 N0. For asymptotic investigations when N1 is much smaller than N0, it is more appropriate to assume that N1 increases in a slower rate compared with N0, i.e., N1/N0 → 0 in probability as N → ∞ [30]. This requires Pr(y = 1) = E{p(x;θ)} → 0 as N → ∞ on the model side. We assume that the parameter θ contains two components θ = (α,βT)T and the log odds can be written as
g(x;θ) := log { p(x;θ) 1− p(x;θ) } = α+ f(x;β),
where the true parameter, denoted as θ∗ = (α∗,β∗T)T, satisfies that α∗ → −∞ as N → ∞ and β∗ is fixed. Here f(x;β) is a smooth function of β, such as a linear function or a nonlinear neural network model. If it is linear, the model reduces to the logistic regression model. A diverging α∗ and a fixed β∗ indicate that both the marginal and conditional probabilities for a positive instance are small. Heuristically, this implies that a change in the feature value does make a rare event become a large probability event. We can also allow β∗ to change with N , but as long as β∗ has a finite limit, the situation is essentially the same as a fixed β∗. So here we assume β∗ to be fixed to simplify the presentation.
Throughout this paper, we use ġ(x;θ) and g̈(x;θ) to denote the gradient and Hessian of g(x;θ) with respect to (w.r.t.) θ, respectively. Let ‖v‖ be the Frobenius norm and v⊗2 = vvT for a vector or a
matrix v. We denote π(x, y) as the sampling probability for an instance (x, y). For negative cases, we sometimes use a shorthand notation π(x) to denote its sampling probability. We use oP (1) to represent a random quantity that converges to 0 in probability as N →∞.

heading: 2.1 It is the number of positive instances that matters
text: We show that for rare events data, the estimation error rate for θ is related to N1 instead of N . Using the full data, a commonly used estimator of θ is the maximum likelihood estimator (MLE)
θ̂f := argmax θ N∑ i=1 [ yig(xi;θ)− log { 1 + eg(xi;θ) }] .
To investigate the theoretical properties of the estimator θ̂f , we make the following assumptions.
Assumption 1. The first, second, and third derivatives of f(x;β) and ef(x;β)f(x;β) w.r.t. any components of β are bounded by a square intergrable random variable B(x).
Assumption 2. The matrix E{ġ⊗2(x;θ∗)} is finite and positive definite.
Assumption 1 imposes constraints on the smoothness of f(x;θ) w.r.t. θ and on the tightness of the distribution of x. Assumption 2 ensures that all components of θ are estimable.
Theorem 1. Let Vf = E{ef(x;β ∗)}M−1f where Mf = E{ef(x;β ∗)ġ⊗2(x;θ∗)}. Under Assumptions 1 and 2, as N →∞,√
N1(θ̂f − θ∗) −→ N(0, Vf), in distribution.
This result shows that even if all the N instances are used for model training, the resulting estimator θ̂f converges to the true parameter θ∗ at the rate of N −1/2 1 , which is much slower than the usual rate of N−1/2 for regular cases. This indicates that with rare events data, the available information about unknown parameters is at the scale of N1. Although N →∞ much faster than N1 →∞ (in terms of N1/N → 0), N does not explicitly show in the asymptotic normality result of Theorem 1. For the specific case of linear logistic regression with g(x;θ) = α + xTβ, our Theorem 1 reduces to Theorem 1 of [30].

heading: 2.2 Negative sampling
text: Since the available information ties to the number of positive instances instead of the full data size, one can keep all the positive instances and significantly subsmaple the negative instances to reduce the computational cost. For better estimation efficiency, we consider nonuniform sampling. Let ϕ(x) > 0 be an integrable function and ρ be the sampling rate on the negative class. Without loss of generality, assume E{ϕ(x)} = 1 so that π(xi) := ρϕ(xi) is the sampling probability for the i-th data point if yi = 1. We present the nonuniform negative sampling procedure in Algorithm 1.
Algorithm 1 Negative sampling For i = 1, ..., N : if yi = 1, include {xi, yi, π(xi, yi) = 1} in the sample;
if yi = 0, calculate ϕ(xi) and generate ui ∼ U(0, 1); if ui ≤ ρϕ(xi), include {xi, yi, π(xi, yi) = ρϕ(xi)} in the sample.

heading: 3 Weighted estimation and its optimal negative sampling probability
text: Let δi = 1 if the i-th data point is selected and δi = 0 otherwise. Given a subsample taken according to π(xi, yi) = yi + (1− yi)ρϕ(xi), i = 1, ..., N , the IPW estimator of θ is
θ̂w = argmax θ N∑ i=1 δiπ −1(xi, yi) [ yig(xi;θ)− log { 1 + eg(xi;θ) }] . (2)
We need an assumption on the subsampling rate to investigate the subsample estimators.
Assumption 3. The subsampling rate ρ satisfies that cN := eα ∗ /ρ → c for a constant c such that 0 ≤ c <∞.
Note that cN cannot be zero but its limit c can be exactly zero. It can be shown that cNE{ef(x;β ∗)} = N1(N0ρ) −1{1 + oP (1)}, so cE{ef(x;β
∗)} is the asymptotic positive/negative ratio for the sample. The theorem below shows the asymptotic distribution of θ̂w. Theorem 2. Under Assumptions 1-3, if E[{ϕ(x) + ϕ−1(x)}B2(x)] <∞ then as N →∞,√
N1(θ̂w − θ∗) −→ N(0, Vw), in distribution,
where Vw = Vf + Vsub and Vsub = cE{ef(x;β ∗)}M−1f E{ϕ−1(x)e2f(x;β ∗)ġ⊗2(x;θ∗)}M−1f .
The subsampled estimator θ̂w have the same convergence rate as the full data estimator. However, the asymptotic variance Vw may be inflated by the term Vsub from subsampling the negative cases. Here Vsub is zero if c = 0. If we keep much more negative instances than the positive instances in the subsample (c = 0), then there is no asymptotic information loss due to subsampling (Vw = Vf ). In this scenario, the number of negative instances can still be aggressively reduced (ρ→ 0) so that the computational efficiency is significantly improved, and the sampling function ϕ(x) does not play an significant role. If we have to reduce the negative instances to the same level of the positive instances (0 < c < ∞), then the variance inflation Vsub is not negligible. In this scenario, a well designed sampling function ϕ(x) it is more relevant and critical. Theorem 2 shows that the distribution of the estimation error err := √ N1(θ̂w − θ∗) is approximated by a normal distribution N(0,Vw). Thus for any > 0 the probability of excess error Pr(err > ) is approximated by Pr{‖N(0,Vw)‖ > } = Pr (∑d j=1 λjz 2 j > ) , where λj’s are eigenvalues of Vw and zj’s are independent standard normal random variables. Therefore, a smaller trace of Vw means a smaller probability of excess error of the same level since tr(Vw) = ∑d j=1 λj .
The following theorem gives the optimal negative sampling function for the IPW estimator. Theorem 3. For a given sampling rate ρ, the asymptotically optimal ϕ(x) that minimizes tr(Vw) is
ϕos(x) = min{t(x;θ∗), T}
E[min{t(x;θ∗), T}] , (3)
where t(x;θ) = p(x;θ)‖M−1f ġ(x;θ)‖ and T is the maximum number so that ρ[min{t(x;θ∗), T}] ≤ E[min{t(x;θ∗), T}] with probability one. Remark 1. If ρ→ 0 then ρt(x;θ∗) ≤ E{t(x;θ∗)} almost surely, so T can be dropped (i.e., T =∞). If limN→∞ ρ > 0 then c = 0, so the variance inflation due to subsampling is negligible. Thus, the truncation term T can be ignored in practice with imbalanced data; it only plays an role for not very imbalanced data when the sampling ratio is not very small.
In ϕos(x), t(x;θ∗) consists of two components, p(x;θ∗) and ‖M−1f ġ(x;θ∗)‖. The term p(x;θ∗) is from the binary response model structure, and it gives a higher preference for a data point with larger p(x;θ∗) in the negative class. The term ‖M−1f ġ(x;θ∗)‖ corresponds to the A optimality criterion [24]. For computational benefit in optimal sampling, the A optimality is often replaced by the L optimality [e.g., 32, 28]. Under the L optimality, ‖M−1f ġ(x;θ∗)‖ is replaced by ‖ġ(x;θ∗)‖, and this minimizes the asymptotic mean squared error (MSE) of Mf θ̂w. In case the gradient is difficult to obtain, one could ignore the gradient term and simply use p(x;θ∗) to replace t(x;θ∗). Although this does not give an optimal sampling probability, it outperforms the uniform subsampling because it takes into account the information of the model structure like the local case-control (LCC) sampling.

heading: 4 More efficient estimation based on the likelihood with log odds correction

heading: 4.1 Nonuniform log odds correction
text: The optimal function ϕos(x) assigns larger probabilities to more informative instances. However, the IPW estimator in (2) assigns smaller weights to more informative instances in estimation, so the resulting estimator can be improved to have higher estimation efficiency. A naive unweighted
estimator is biased, and the bias correction approach in [8, 29] does not work for the sampling procedure in Algorithm 1 even when the underlying model is the logistic regression. We adopted the idea of [11] and seek to define a more efficient estimator through finding the corrected likelihood of the sampled data based on any negative sampling probability π(x). For data included in the subsample (where δ = 1), the conditional probability of y = 1 is
pπ(x;θ) := Pr(y = 1 | x, δ = 1) = 1
1 + e−g(x;θ)−l , (4)
where l = − log{π(x)}. Please see the detailed derivation of (4) in the supplement. To avoid the IPW, we propose the estimator based on the log odds corrected likelihood in (4), namely,
θ̂lik = argmax θ `lik(θ), where `lik(θ) = N∑ i=1 δi [ yig(xi;θ)− log { 1 + eg(xi;θ)+li }] . (5)
With θ̂w in (2), π(xi) is in the inverse. If an instance with much smaller π(xi) is selected in the sample, it dominates the objective function, making the resulting estimator unstable. With θ̂lik, this problem is ameliorated because π(xi) is in the logarithm in the log-likelihood of (5) and log(v) goes to infinity much slower than v−1 as v ↓ 0.

heading: 4.2 Theoretical analysis of the likelihood-based estimator
text: The following shows the asymptotic distribution of θ̂lik.
Theorem 4. Under Assumptions 1-3, if E{ef(x;β∗)ϕ−1(x)B(x)} <∞, then as N →∞,√ N1(θ̂lik − θ) −→ N(0,Vlik) in distribution,
where Vlik = E{ef(x;β ∗)}Λ−1lik and Λlik = E
[ ef(x;β
∗)ġ⊗2(x;θ∗) 1+cϕ−1(x)ef(x;β∗)
] .
Theorem 3 shows that the proposed estimator θ̂lik is asymptotic normal with variance Vlik. We compare Vlik with Vw to show that θ̂lik has a higher estimation efficiency than θ̂w.
Theorem 5. If Vw and Vlik are finite and positive definite, then Vlik ≤ Vw, i.e., Vw −Vlik is non-negative definite. The equality holds when c = 0 and in this case Vlik = Vw = Vf .
Since the estimator θ̂lik in Eg. (5) is based on the conditional log-likelihood of the sampled data, it actually has the highest estimation efficiency among a class of asymptotically unbiased estimators.
Theorem 6. Let X = (x1, ...,xN )T be the feature matrix and Dδ denote the sampled data. Let θ̂U be a subsample estimator with the following asymptotic representation:
θ̂U = U(θ ∗;Dδ) +N−1/21 oP (1), (6)
where the variance V{U(θ∗;Dδ)} exists, and U(θ∗;Dδ) satisfies that E{U(θ∗;Dδ) | X} = θ∗ and E{U̇(θ∗;Dδ) | X} = 0 with U̇(θ∗;Dδ) = ∂U(θ∗;Dδ)/∂θ∗T. If N1V{U(θ∗;Dδ)} → VU in probability, then VU ≥ Vlik.
Remark 2. Theorem 6 tells us that θ̂lik is statistically the most efficient among a class of estimators, which includes both θ̂w and θ̂lik as special cases. The condition E{U(θ∗;Dδ) | X} = θ∗ ensures that the estimator is asymptotically unbiased. The condition E{U̇(θ∗;Dδ) | X} = 0 can be intuitively interpreted as the requirement that the derivative of the oP (1) term in (6) is also negligible. Clearly, this is satisfied by all unbiased estimators for which the oP (1) term in (6) is 0.

heading: 5 Practical considerations
text: In this section, we first show how we design a more practical estimator based on the previous results and then present the theoretical analysis behind these improvements.

heading: 5.1 Making estimators more practical
text: Like the LCC sampling [8] and the A-optimal sampling [29], the ϕos(x) in (3) depends on θ, and thus a pilot value of θ, say θ̃, is required in practice. Here, θ̃ can be constructed from a pilot sample, say (x̃i, ỹi), i = 1, ..., ñ, taken by uniform sampling from both classes with equal expected sample sizes. With θ̃ obtained, calculate t̃i = p(x̃i; θ̃)‖M̃−1f ġ(x̃i; θ̃)‖ for i = 1, ..., ñ, where M̃f is the Hessian matrix of the pilot sample objective function. In case that the Hessian matrix and gradients are hard to record, the gradient term can be dropped, i.e., use t̃i = p(x̃; θ̃). As mentioned in Remark 1, T can be ignored for very imbalanced data. The expectation in the denominator of (3) can be approximated by mimicking the method of moment estimator, namely by
ω̃ = 2N1(ñN) −1∑
ỹi=1 t̃i + 2N0(ñN) −1∑ ỹi=0 t̃i. (7)
In practice, it is common to set a lower threshold, say %, on sampling probabilities, to ensure that all instances have positive probabilities to be selected. This is more critical for the IPW estimator while the likelihood estimator is not very sensitive to very small sampling probabilities.
Denote the pilot value as ϑ̃ = (θ̃, ω̃). The following probabilities can be practically implemented
πos% (x; ϑ̃) = min[max{ρϕ̃os(x), %}, 1], where ϕ̃os(x) = ω̃−1t(x; θ̃), (8)
and we use this in our experiments in Section 6. The ω̃ in (7) is essentially a normalizing term so that the expected subsample size is ρ. For some practical systems such as online models, ω̃ is often treated as a tuning parameter. We will illustrated this in our real data example in Section 6.2.

heading: 5.2 Theoretical analysis of practical estimators
text: Denote θ̂ϑ̃w and θ̂ ϑ̃ lik the IPW estimator and the likelihood estimator, respectively, with practically estimated optimal probability in (8). We have the following results.
Theorem 7. Assume that the pilot estimator ϑ̃ is independent of the data, ϕ̃os(x) → ϕplt(x) in probability, and ρ−1%→ cl > 0. Under Assumptions 1-3, as N →∞, the θ̂ϑ̃w satisfies that,√
N1(θ̂ ϑ̃ w − θ)→ N(0,Vpltw ) in distribution,
where Vpltw has the same expression of Vw except that ϕ(x) is replaced by max{ϕplt(x), cl}. If Pr{ϕplt(x) ≥ cl} = 1 and the pilot estimator is consistent, then Vpltw achieves the optimal variance.
Theorem 8. Assume that the pilot estimator ϑ̃ is independent of the data, ϕ̃os(x) → ϕplt(x) in probability, and % = o(ρ). Under Assumptions 1-3, as N →∞, θ̂ϑ̃lik satisfies that,√
N1(θ̂ ϑ̃ lik − θ)→ N ( 0,Vpltlik ) in distribution,
where Vpltlik has the same expression of Vlik except that ϕ(x) is replaced by ϕplt(x). If the pilot estimator is consistent, then Vpltlik achieves the variance with the optimal probability.
Remark 3. Theorems 7 and 8 do not require ϑ̃ to be consistent, i.e., the pilot estimator can be misspecified. However ϑ̃ has to be consistent in order to achieve the asymptotic variances with the optimal probability. Compared with the likelihood estimator θ̂ϑ̃lik, the IPW estimator θ̂ ϑ̃ w requires a stronger condition on % to have asymptotic normality, because πos% (x; ϑ̃) is in the denominator of the objective function for θ̂ϑ̃w so the turbulence of ϑ̃ is amplified. The requirement cl > 0 means that % cannot be too small in practice. The likelihood estimator does not have this constraint.

heading: 6 Numerical experiments
text: We present simulated results and report the performance on a CTR dataset with more than 0.3 trillion instances. More experiments with model and pilot misspecifications are available in the supplement.

heading: 6.1 Experiments on simulated data
text: We implemented the following negative sampling methods for logistic regression. 1) uniW: uniform sampling with the IPW estimator. 2) uniLik: uniform sampling with the likelihood estimator. 3) optW: optimal sampling with the IPW estimator. 4) optLik: optimal sampling with the likelihood estimator. We also implement the full data MLE (Full) and the LCC for comparisons.
In each repetition of the simulation, we generate full data of size N = 5 × 105 from the logistic regression with g(x;θ) = α+ xTβ. We set the true β∗ to be a 6× 1 vector of ones and set different values for α for different feature distributions so that the positive/negative ratio is close to 1:400. We consider the following feature distributions: (a) Normal distribution which is symmetric with light tails; the true intercept is α = −7.65. (b) Log-normal distribution which is asymmetric and positively skewed; the true intercept is α = −0.5. (c) t3 distribution which is symmetric with heavier tails; the true intercept is α = −7. We set the sampling rate as ρ = 0.002, 0.004, 0.006, 0.01, and 0.02 for all sampling methods. In each repetition, uniform samples of average size 100 are selected from each class to calculate the pilot estimates, so the uncertainty due to the pilot estimates are taken into account. We also consider pilot misspecification by adding a uniform random number from U(0, 1.5) to the pilot so that the pilot is systematically different from the true parameter. We repeat the simulation for R = 1000 times to calculate the MSE as R−1 ∑R r=1 ‖θ̂(r) − θ‖2, where θ(r) is the estimate at the r-th repetition.
Results are presented in Figure 1. The MSEs of all estimators decrease as the sampling rate ρ increases. The optLik outperforms other sampling methods in general, and its performance is close to the full data MLE when ρ = 0.02, especially if the pilot is consistent. The advantage of optLik over optW is more evident for smaller ρ. If the pilot is misspecified, optLik is quite robust but optLik and LCC significantly deteriorate. Since uniLik and uniW do not require a pilot, they are not affected by pilot misspecification, but they have much larger MSE than optLik and optW when the pilot is consistent, especially uniW. Different feature distributions also affect the performances of different sampling methods.
To verify Theorem 1 numerically, we let the full data size N increase at a faster rate than the average number of positive instances Na1 = E(N1) so that the probability Pr(y = 1) decreases towards zero. As N increases, we fixed the value of β∗ and set decreasing values to α∗ to mimic our scaling regime. We considered two covariate distributions: 1) a multivariate normal distribution for which a logistic regression model is a correct model; and 2) a multivariate log-normal distribution for which a logistic regression model is misspecified. We simulated for 100 times to calculate the empirical variance V̂f of the full data MLE, and report the results in Table 1. It is seen that tr(V̂f ) is decreasing towards zero. According to Theorem 1, tr(V̂f ) should converge in a rate of 1/N1, so Na1 tr(V̂f ) should be relatively stable as N increases. This is indeed the case as seen in the third and sixth
columns of Table 1. On the other hand, Ntr(V̂f ) gets large dramatically as N increases. The aforementioned observations confirm the theoretical result in Theorem 1. Furthermore, Table 1 shows that the convergence rate in Theorem 1 may also be true for some misspecified models.

heading: 6.2 Experiments on real data
text: We conduct experiments on an internal CTR dataset from our products with more than 0.3 trillion instances. There are 299 input features, including user and ad profiles together with rich context features. We concatenate all features as a single high-dimensional vector x, whose length is greater than 5, 000. We use a 3-layer fully connected deep neural network as a feature extractor that maps x into a 256-dimensional dense features h(x). The downstream model is a linear logistic regression model taking h(x) as input to predict the binary output whether the user clicks the ad or not.
The positive/negative ratio is around 1:80. Due to limited storage, we first uniformly drop 80% negative instances and try various negative sampling algorithms on the rest 60 billion instances. On this pre-processed dataset, we then split out 1 percent of the data for testing, and do negative sampling on the rest data. The entire training procedure scans the subsampled data set from begin to end in one-pass according to the timestamp. The testing is done all the way through training. We calculate the area under the ROC curve (AUC) for testing instances. The entire model is trained on a distributed machine learning platform using 1,200 CPU cores and can finish within two days.
We adopt the sampling probability πos% (x; ϑ̃) as demonstrated in (8) and use the empirical sampling rate to calibrate the normalizing term ω̃. As mentioned before, we remove the gradient term for implementation tractability so t(x; θ̃) = p(x; θ̃). As the result will show, this approximate score produces consistent results with the simulation experiments. We compare four methods: 1) uniW and 2) uniLik are the same as in simulated experiments; and 3) optW and 4) optLik represent nonuniform sampling probability using πos% (x; ϑ̃) with the IPW estimator and with the likelihood estimator. We use “opt" to refer to both optW and optLik when focusing on the sampling.
We study two sets of negative sample rates (w.r.t. the original 0.3 trillion data): (i) [0.01, 0.05]. In this case, negative instances still outnumber positive instances. (ii) [0.001, 0.008]. In this case, negative subsampling is too aggressive that positive instances dominates negative instances. We demonstrate the result in Figure 2. When the sample rate is moderate (Figure 2.(b)), opt is nearly optimal and can significantly outperform uniform subsampling. There is still a small gap between the IPW estimator and the likelihood estimator. When the sample rate is extremely small (Figure 2.(a)), opt becomes sub-optimal and its gap w.r.t. uniform sampling is closer. Still we see a clear difference between the IPW estimator and the likelihood estimator. This is due to a larger asymptotic variance for uniform sampling when the sampling rate goes to zero. Note that on this huge data set, a small relative AUC gain (e.g., 0.0005) usually could bring about significant revenue gain in products.

heading: 6.3 Sensitivity analysis
text: In practice, we find it necessary to use cross-validation to find the optimal % in (8) given a fixed sample rate. As demonstrated in Figure 3, for each sampling rate we tune % from a very small value (in this case, 10−6) all the way to its largest value. When % achieves its largest value, ρω̃−1 = 0 and πos% (x; ϑ̃) reduces to uniform negative sampling. We observe that πos% (x; ϑ̃) with a moderate % achieves the best results for various sampling rates. Note that % is the lower truncation level for negative sampling probabilities. Thus, a larger % means more “simple" negative instances (ρω̃−1p(xi; θ̃) is smaller than
%) will have better chances to be selected. This coincides with recent empirical results for negative subsampling in large search systems [14], where they mix “simple" (smaller p(xi; θ̃)) and “hard" (larger p(xi; θ̃)) negative instances. We also observe an empirical variance of around 0.0001 for each setup, demonstrating that the relative improvement is consistent.

heading: 7 Conclusions, discussion, and limitations
text: In this paper, we have derived asymptotic distributions for full data and subsample estimators with rare events data. We have also found the optimal sampling probability that minimizes the unconditional asymptotic variance of the IPW estimator, and proposed the likelihood estimator through nonuniform log odds correction to further improve the estimation efficiency. The optimal probability depends on unknown parameters so we have discussed practical implementations and investigated the convergences of the resultant estimators. Experiments on both simulated and real big online streaming data confirm our theoretical findings. Heuristically, our nonuniform negative sampling method gives preference to data points that are harder to observe, so it would give preference to sub-populations that are difficult to observe. We do not think this would bring any negative societal impacts. We have also examined our real data application and did not notice discrimination against any sub-populations.
This paper has the following limitations: 1) We assume that the model is correctly specified. Theoretical properties with model misspecification is not considered, and this important question requires future studies. 2) Oversampling the rare positive instances is another common practice to deal with imbalanced data. Its application together with negative sampling is not considered. 3) The assumption of a fixed β∗ means that both the marginal and conditional probabilities for a positive instance to occur are small. If the proportion of the positive instances for a given feature in a dataset is large,
this assumption may not be appropriate. Further investigation is need to see if our results is still applicable. In the following, we present three scenarios that the scaling regime may or may not fit in:
Scenario 1 (phone call while driving) Car crashes occur with small probabilities, and making phone calls while driving significantly increases the probability of a car crash. However, the probability of car crashes among people making phone calls while driving is still small, so for these types of features, our scaling regime is appropriate to model the rare events.
Scenario 2 (anti-causal learning) Anti-causal learning [26, 16] assumes that label (y) causes observations (x). Thus Pr(x|y) represents the causal mechanism and is fixed. One standard example is that diseases (y) cause symptoms (x). Our scaling regime fits the framework of anti-causal learning. To see this, using Bayes Theorem we write the log odds as
α+ f(x;β) = log Pr(y = 1)
Pr(y = 0) + log Pr(x | y = 1) Pr(x | y = 0) .
In anti-causal learning, the marginal distribution of y changes while the conditional distribution of x given y is fixed. Thus only the scale factor α changes, and f(x;β) is fixed.
Scenario 3 (local transmission of the COVID-19) Our scaling regime has some limitations and does not apply to all types of rare events data. For example, although the COVID-19 infection rate for the whole population is low, the infection rate for people whose family members are in the same house with positive test results is high. This means the change of a family member’s test result converts a small-probability-event to a large-probability-event, and our scaling regime would not be appropriate.
Acknowledgments and Disclosure of Funding
HaiYing Wang’s research was partially supported by NSF grant CCF 2105571.

Review summaries:
This paper considers the estimation of the parameters of a binary response model when the positive instances are rare. Rareness is modeled via a scaling regime on the logits, with the key restriction being that the scaling occurs _uniformly over features_. The paper first shows that the fraction of positive instances bottlenecks estimation and uses this to motivate subsampling the negative instances, to save on computation. The paper considers two negative subsampling approaches: importance weighting the original likelihood and maximizing it (IPW), or maximizing the log-likelihood under subsampling (Lik). Under correct model specification, asymptotic analysis is provided for both methods, with the latter enjoying many advantages, including stability and smaller asymptotic variance. The approach is tested on synthetic and click-through-rate data. The experiments support the theory and go beyond it, by showing that the methods work well even under model misspecification.

The paper discusses the sampling strategy for imbalanced binary-label datasets. For a binary-label dataset where the number of negative samples $N_0$ is much larger than that of positive samples $N_1$, we often "subsample" a part of negative samples in order to balance two labels.  The paper discusses the asymptotic convergence rate of the optimal model parameter against the number of samples in the dataset, and improves the rate from conventional model parameter estimation method. -   Theorem 1 proved that the convergence rate is proportional to $N_1^{-1/2}$, that is, determined not by the number of all instances but of positive instances if $N_0\gg N_1$. -   Theorems 2 and 3 composes a function that provides the best strategy to subsample negative samples. -   Theorems 4 through 6 modifies the result of Theorems 2 and 3 for improving the efficiency by incorporating the sampling strategy. -   Theorems 7 and 8 discusses the effect of the "pilot" model parameter and provides the construction of it for better convergence rate.

This paper considers the question of parameter estimation in a binary response model (e.g., logistic regression), in the imbalanced case where the fraction of positive instances is small/vanishing. This question is formulated by considering a well-specified model, and adding an intercept to the log-odds ratio which diverges as the sample size goes to infinity (with the other parameters kept fixed).  The main results are the following:  - First, it is shown (in Theorem 1) that the convergence rate of the MLE (computed on the full dataset) only depends on the number N_1 of positive instance, rather than on the sample size (which also includes negative instances). MLE also shown to be asymptotically normal, and the asymptotic covariance is provided.  - Since the convergence rate of MLE only depends on N_1 (but not on the total number N of samples), it can make sense to perform "negative sampling", namely to reduce N by only retaining a small fraction of negative instances, in order to reduce computations. Of course, removing negative instances creates a bias, so the authors considers inverse probability weighted (IPW) estimators, where negative samples are re-weighted (inversely proportionally to the deletion probability). An asymptotic normality result is obtained, with a covariance larger than MLE (as expected), by an additional covariance depending on the rejection probability function (Theorem 2). The optimal rejection function (minimizing asymptotic variance) is deduced (Theorem 3).  - Finally, the authors consider a different estimator, based on log-odds correction. This estimator corresponds to MLE on the selected data, based on adjusting the likelihood of the selected data by accounting for the deletion mechanism. This leads to a correction which is less sensitive to small probability instances than the IPW method. The asymptotic normality and covariance of this estimator is provided (Theorem 4). Interestingly, it is shown (Theorem 5) that this scheme is always better than the IPW (with optimal rejection function), in that it leads to a smaller asymptotic covariance. An optimality result of this estimator among (asymptotically) unbiased estimators is then provided in Theorem 6.