{'fraction_non_alphanumeric': 0.0541261315403974, 'fraction_numerical': 0.022273973631822108, 'mean_word_length': 4.129865125240848, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining case-control studies for identifiability and efficiency improvement in logistic regression

Can two separate case-control studies, one about Hepatitis disease and the other about Fibrosis, for example, be combined together? It would be hugely beneficial if two or more separately conducted case-control studies, even for entirely irrelevant purposes, can be merged together with a unified analysis that produces better statistical properties, e.g., more accurate estimation of parameters. In this paper, we show that, when using the popular logistic regression model, the combined/integrative analysis produces a more accurate estimation of the slope parameters than the single case-control study. It is known that, in a single logistic case-control study, the intercept is not identifiable, contrary to prospective studies. In combined case-control studies, however, the intercepts are proved to be identifiable under mild conditions. The resulting maximum likelihood estimates of the intercepts and slopes are proved to be consistent and asymptotically normal, with asymptotic variances achieving the semiparametric efficiency lower bound.


Introduction
Logistic regression (Cox, 1958) is probably the most fundamental statistical tool to model categorical dependent variable. Let Y be a categorical response variable and let X be an observable d-dimensional vector of covariates. Without loss of generality, we focus on binary response typically coded as 0/1. A binary logistic regression model assumes where β ∈ R d is the slope parameter and α ∈ R is the intercept. The logistic regression is often regarded as a special case of the generalized linear models, which has been extremely popular in biomedical science for the study of the effect of certain exposure to possible disease or hazards. The response is thus often binary, e.g, Y = 0 represents the controls (non-disease) and Y = 1 represents the cases (disease), or there may be more than two categories of Y . In prospective studies, a sample of subjects or individuals is followed and their respective Y are recorded. For model (1) under prospective studies, the samples are simple random sampling from the underlying population. Thus α, β can be consistently estimated by the maximum likelihood estimation (MLE). However, many diseases are often fairly rare such that P (Y = 0|X = x) is close to 1. Large studies may produce very few diseased cases and thus very little information about the disease of interest.
Case-control study is a primary tool for the study of existing factors related to rare disease incidence, by taking samples separately from the case population and the control population, when the case population and the control population are clearly separated, e.g, through disease registry records. In econometrics, people are interested in the relationship of the covariates and the choices made by individual, e.g, the choice of transportation from Shanghai to Hong Kong. It would be easier and cheaper to take samples of individuals from the Hong Kong international airport, West Kowloon High-speed Train Station and Hung Hom train station, which is a choice-based sampling data (Manski and McFadden, 1981), than to take a single sample from the entire population. Generally speaking, according to Lawless(1997), sampling schemes that depend on the value of the outcome/response are called response-selective to response-biased sampling, among which case-control sampling is the most typical one; see Manski and Lerman (1977), Prentice and Pyke (1979), Cosslet (1981), Wild (1986, 1997), Manski (1993), Chen and Lo (1999), Chen (2001), Fithian and Hastie (2014), Liu et al. (2014) and Chen et al. (2017) etc. Another popular sampling design in epidemiology is the case-cohort sampling. The estimation and inference of case-cohort data with Cox's proportional hazards model, transformation models and other semiparametric models are studied by Prentice (1986), Lu and Tsiatis (2006), Cai andZeng (2004, 2007), Liu et al. (2010), Zeng and Lin (2014), Ni et al. (2016), Tao, Zeng and Lin (2017), etc. Unified estimations for semiparametric linear transformation models, the accelerated failure time model and quantile regression under general biased sampling schemes were studied by Kim, Lu, Sit and Ying (2013), Kim, Sit and Ying (2016) and Xu et al. (2017). In fact, earlier work on nonparametric inference were developed for recovering the distribution function in two-sample problem in the presence of selection bias under known selection bias weight function; see Vardi (1982Vardi ( , 1985, Qin (1993) and Qin and Zhang (1997). When the population case percentage is known, large sample theory of the nonparametric maximum likelihood estimate for semiparametric biased sampling data was established by Gilbert (2000). A comprehensive discussion on biased sampling and associated parameter problems can be found in Qin (2017).
For logistic regression under case-control sampling, a remarkable result is that, the prospective estimating equation derived from the MLE is valid for a consistent estimate of the slope parameter, except the intercept term (Prentice and Pyke, 1979). The beauty of this method is its simplicity and ease of inference. We first point out that for model (1) under a single case-control study as in the literature, the intercept α and f (·) are not identifiable, where f (·) is the probability density function of the covariate X in the population. The rationale behind is that, with a single case-control study, the score function of α lies in the linear space spanned by the score function of f (·). This is discussed in details in Remark 2 in section 2.1. As a result, contrary to prospective studies, the intercept cannot be consistently estimated in a single case-control study. Nonetheless, a consistent estimate of the intercept is important in many machine learning problems such as image classification or pattern recognition, and it is crucial for better understanding of the incidence of the disease, when prior knowledge of the event/disease is not available or not reliable. Despite some nice theory on the estimation of β is developed under various circumstances, no specific construction of consistent estimation of α in case-control logistic regression is available in the literature. Moreover, the consistent estimate for β in single case-control study is not semiparametric efficient.
In this paper, we propose an efficient estimation for the intercepts and slope parameters of the logistic regression models under multiple/distributed case-control studies. With K case-control studies of different purposes from the same population or collected from K heterogeneous subpopulations, the data can be pooled together for integrative analysis. The intercepts and slope parameters in the logistic regression models can vary across different studies/subpopulations. By combining K case-control studies, surprisingly, the intercepts become identifiable under mild conditions (see Proposition 1 in section 2.1), as the score function of each intercept term no longer lies in the linear space spanned by the score of f (·) (see section 2.2 for detailed discussions); most importantly, the resulting estimates of the slope parameters are proved to be semiparametric efficient, with asymptotic variance smaller than the counterparts in single case-control study (see Remark 4 in section 2.2 and Remark 5 in section 2.4). The proposed estimation is based on the maximisation of the nonparametric likelihood function of the integrative data. The resulting estimates for the intercepts and slope parameters are shown to be consistent, asymptotically normal and asymptotically efficient. An iterative algorithm is employed to compute the maximum likelihood estimates numerically.

Theory and Methods
Without loss of generality, we focus on binary outcome. Suppose that there are K independent studies that might be of different purposes, or there are K heterogeneous subpopulations. For k = 1, . . . , K, the d-dimenisonal explanatory variables X and the binary response Y k satisfy where φ(x) = e x /(1 + e x ) and θ = (α 1 , . . . , α K , β ⊤ 1 , . . . , β ⊤ K ) ⊤ of dimension p = K + Kd. Let f (·) and F (·) be the density function and the cumulative distribution of X respectively. We assume that f (·) and F (·) are the same for k = 1, . . . , K in model (2).
For each k = 1, . . . , K, let f k1 (·) be the density function of the case population, i.e, the conditional distribution of X given Y k = 1, and f k0 (·) be that of the control population, i.e, the conditional distribution of X given Y k = 0. Suppose that the k-th case-control study is conducted by taking a random sample of n k1 cases from its case population, and a random sample of n k0 controls from its control population, denoted by D k = {(y ik , x ik )} n k i=1 . Note that n k1 and n k0 are pre-specified in case-control studies. Let n k = n k1 + n k0 and N = K k=1 n k . Thus, the pooled data are D = {D k , k = 1, . . . , K} consisting of K case-control studies. Under model (2), by the Bayes' formula, for k = 1, · · · , K, where c k (θ) = +∞ −∞ φ k (θ, x)f (x)dx is the population percentage of cases for the k-th study or the k-th subpopulation. For notational simplicity, we simply write c k (θ) as c k , k = 1, . . . , K. Let β k,0 and α k,0 be the true values of β k and α k , k = 1, . . . , K, respectively.

Identifiability
Throughout the paper, parameters including f (·) are not identifiable if two different choices have the same resulting distribution of random variables. In classical logistic regression model (1) under single case-control study (Scott and Wild, 1986), that is K = 1 in model (2), for any two values x and x 0 in the support of X, the logarithmic odds ratio (OR) is where f (·|Y = 1) is f 1 (·), the conditional density of X in the case population defined earlier.
Heuristically, independent samples from the case population and control distribution are observed separately in case-control studies, making β identifiable according to (4); but α is not involved in (4) and thus it is not identifiable under single case-control study.
With multiple/distributed case-control studies, however, the data can be pooled together for unified analysis. Under regularity conditions (C1)-(C3) given in section 2.4, the idea of combining data enables us to consistently estimate the intercepts. We discuss the identifiability of θ under multiple case-control studies in the following proposition.
As a result, f (·) is identifiable if and only if one of the α k is identifiable.
(d) Assume f is continuous, and β k are not all same. Then, those α k with β k = 0 are not identifiable and those with β k = 0 are not identifiable.
Part(a) is consistent with the results in Prentice and Pyke (1979). Part(b) tells the relationship between the identifiability of the intercepts and f (·). Part(c) has an important implication that when each separate study/subpopulation shares the same/common intercept and slope parameters, the intercept is still not identifiable despite the data combining. In other words, when a single case-control study is sub-divided into multiple smaller studies, combining those smaller case-control studies cannot make the intercept identifiable. Part(d) implies that those sub-models with β k = 0, their respective α k is always not identifiable, regardless of the identifiability of f (·). For the rest of this paper, we assume that

The Score Functions and Maximum Likelihood Estimation
We next introduce our proposed maximum likelihood estimation. The likelihood function of and the log-likelihood function is According to the profile-likelihood method in Lin (2006, 2007), the nonparametric component f (·) in (6) can be profiled over its observed values. To avoid uninteresting discussions, we assume there is no tie in the realizations of X. With slight abuse of notation, we still use x ik , i = 1, . . . , n k , k = 1, . . . , K to denote the realizations. We define the estimator of F (·) as a step function with jumps only at the distinct observed values subject to Our proposed maximum likelihood estimator for θ and p ik , denoted byθ andp ik , is defined as the maximizer ofl n (θ, p; D) under the constraints in (8). Thus, the nonparametric maximum likelihood estimator (NPMLE) for F (x) is given byF n (x) = x ik ≤xp ik .
Remark 1 Let vec(A) be the vectorization of matrix A of order m × n, namely, where e t is a K-vector with the t-th element being 1 and other elements being 0. And the score function of f , denoted by S f [g 1 ], is the partial derivative of l(θ, f ; D) along the direction (θ, f ε = f + εg 1 ) with ε being a small constant, such that f ε ≥ 0 and Here, BV [D] is a class of functions on domain D with bounded total variation. Then, Given the observations x ik , i = 1, . . . , n k , k = 1, . . . , K, the log-likelihood l(θ, f ; D) can be approximated byl(θ, p; D); and the score function of p ik , i = 1, . . . , n k , k = 1, . . . , K, based on the profile likelihoodl(θ, p; D) is Moreover, the linear space spanned by S * f (p ik ), i = 1, . . . , n k , k = 1, . . . , K, is a N-dimensional surface on S f (Bickel et al. 1993). As a result, the linear space spanned by S f contains the linear space spanned by S where the equality holds by invoking (9). This observation reveals that S α and S * f are on the same space when K = 1; in other words, S α is on the linear space spanned by S f when K = 1. As a result, the intercept and the population percentage of cases are not identifiable in single case-control logistic regression.
Remark 2 When K ≥ 2, we focus on K = 2 without loss of generality.
• When β 1 = β 2 , β 1 = 0 and β 2 = 0, the score functions of α 1 and α 2 are respectively. Similar to Remark 1, direct calculations yield that the surface of the score function S f is Apparently, the score functions of α 1 and α 2 are not on any spanned linear space of S f when β 1 = β 2 , β 1 = 0 and β 2 = 0. This fact offers insights that the intercepts α 1 , α 2 can be identifiable by combining the two case-control studies, when β 1 and β 2 are different and nonzero. Thus, α 1 , α 2 , β 1 , β 2 and f (·) in model (2) can be consistently estimated by the proposed maximum likelihood estimation.
• When β 1 = 0 and β 2 = 0, the score function of α 1 is (11) and that of α 2 is indicating that S α 2 is in the linear space spanned by S f , but S α 1 is not. In other words, α 2 is unidentifiable but α 1 is identifiable. This is consistent with Proposition (1d) in section 2.1.
• Lastly, for the special case that β 1 = 0 and β 2 = 0, the score functions of α 1 and α 2 are respectively. It is clear that S α 1 and S α 2 are on the spanned linear space of S f . As a result, α 1 and α 2 are unidentifiable when β 1 = 0 and β 2 = 0.
Remark 3 We need to emphasize that, β k in model (2) can be always separately estimated by the celebrated estimation in Scott and Wild (1986) or Qin (2017, pages 210-212) with the k-th case-control study, k = 1, . . . , K, though α k is unidentifiable. In contrast to the classical separate estimation in the literature, our proposed method based on the integrative data is able to identify the intercepts and give more efficient estimation of β k , k = 1, . . . , K, which achieves the semiparametric efficiency lower bound, as evidenced in the simulation studies.

An Iterative Algorithm
Direct maximization of (7) is challenging and potentially unstable as the second constrain in (8) is intractable. To circumvent the difficulty, we first plug inc t = K k=1 n k i=1 φ t (θ, x ik )p ik , t = 1, . . . , K into (7), and consider to maximize a Lagrange function by incorporating the first constraint of (8) as follows where λ > 0 is the Lagrange multiplier. Taking first derivative of (12) with respect to p ik and θ, we obtain ∂l(θ, p; D) and Next, multiplying p ik on both sides of (13), one can easily verify that λ = 0 because of N = K t=1 n t . Consequently, We now devise an iterative algorithm to maximize (12). Set the initial value of the parameter ik , k = 1, . . . , K, i = 1, . . . , n k }. For the j-step, given p (j) , we update θ (j) by solving ik . Notice that the objective function in (12) is convex in θ, thus the solution to (16) is unique. Subsequently, given (θ (j) , p (j) ), we update p (j+1) according to the recursive formulâ We iterate between (16) and (17) until convergence. Denote the resulting estimator for θ and p byθ andp. Thus,F n (x) ≡ K k=1 n k i=1p ik I(x ik ≤ x). Moreover, the population percentage of casesc t is estimated by The iterative algorithm is summarized as follows:

Algorithm 1 The iterated algorithm for computing the MLEs
Set an initial value p (0) . repeat Step 1. Given p (j) , compute θ (j) by solving (16) with the lgfbs package in R.
Here κ 1 and κ 2 in the stopping criterions are set to be n −6 in our numerical studies.

Asymptotic properties
Let θ 0 ≡ (α 1,0 , . . . , α K,0 , β ⊤ 1,0 , . . . , β ⊤ K,0 ) and Conditions (C1)-(C2) are regularity conditions to ensure the identifiability of the parameters. The continuously distributed condition of X is imposed for technical convenience. It can be weaken to accommodate discrete predictors, as long as it is satisfied that there exist at least three different values x 1 , x 2 , x 3 in the support of X such that β ⊤ k,0 x 1 , β ⊤ k,0 x 2 and β ⊤ k,0 x 3 are also different for those nonzero β k,0 , k = 0, . . . , K. Condition (C3) assumes that θ 0 is an interior point of a compact set. We next present the asymptotic properties in three theorems.
Theorem 1 and Theorem 2 indicate that maximizing the profile log-likelihood (7) over (θ, p) leads to consistent, asymptotically normal and semiparametric efficient estimate for θ under mild conditions. Theorem 3 provides a simple and easy-to-implement variance estimation.
Remark 4 Our proposed estimator is semiparametric efficient under K case-control studies by Theorem 2. To have more insights into the efficiency gain, we consider model (2) with K = 2, and focus on the estimation of the parameter θ 1 ≡ (α 1 , β ⊤ 1 ) ⊤ of the first study as a toy example. Three estimators of θ 1 are considered: the maximum likelihood estimator for single case-control study with known f (·), denoted byθ (1)⊤ 1 ) ⊤ ; the maximum likelihood estimator for single case-control study with unknown f (·) (Scott and Wild, 1986), denoted bŷ θ When f (·) is known in a single case-control study, the score function of θ 1 is When f (·) is unknown in a single case-control study, according to Remark 2, the score of Then, the score function of β 1 where H(·) ∈ G is the projection of ψ 1 − Eψ 1 onto the linear space spanned by S f that minimizes Var{S (2) (β 1 )}. Apparently, Var{S (1) (β 1 )} is always larger than Var{S (2) (β 1 )}. Nonetheless, by combining the two case-control studies, the score function of θ 1 becomes where b 2 is some constant and g(·) ∈ G such that Var{S (3) (θ 1 )} is minimized and nonnegative definite.
Theoretically speaking, the maximum likelihood estimatorθ (1) 1 for single case-control study with known f (·) is most efficient for estimating θ 1 , thus combining another case-control study cannot improve efficiency. However, when f (·) is unknown, combining case-control studies would lead to efficiency improvement in the estimation of θ 1 . For instance, when the sample size of the second case-control study to be combined is much larger than that of the first one, our proposed estimatorθ can be nearly as efficient asθ (1) 1 asymptotically, implyingβ (3) 1 is more efficient thanβ (2) 1 for estimating the slope parameter β 1 . In other words, it can be shown that, where r k = n k /N, k = 1, 2.

Empirical results
Finite-sample studies are carried out to examine the performance of our proposed method under different scenarios. We consider binary outcome and generate independent data from the following models where β k ∈ R d , k = 1, . . . , K. The d-vector predictor X follows the same distribution F (·) across K logistic regression models. Thus, the true population percentage of cases in the k-th subpopulation can be calculated via p k ≡ P (Y k = 1) = φ(α k + β ⊤ k x)dF (x), k = 1, ..., K. The case-control sampling are then conducted by taking samples separately from the case population and the control population from each subpopulation respectively. Let q k = n k1 /n k , k = 1, ..., K, be the case proportion in the k-th case-control study.

Simulated data: verifying identifiability conditions
We investigate the identifiability of K intercepts under various settings. In cases (a1)-(a4), X 1 and X 2 are independent standard normal random variables. For simplicity, set n 1 = n 2 , n 11 = n 10 and n 21 = n 20 . The pooled sample size N = n 1 + n 2 = 500. The following scenarios (a1)-(a6) are tried for K = 2 and d = 2. The results of each case are based on 1000 replications.
In Table 1 and 2, we present the bias of the estimates of the regression parameters and p k , the empirical standard errors (SE), the average of the estimated standard errors (ESE) and the 95% coverage probabilities (CP) with our proposed method. It is seen that for all cases, β is always identifiable. For the identifiability of the intercepts, when the slope parameters are different across two models, α 1 and α 2 are both identifiable for Cases (a1)-(a2) in Table  1. We have to emphasise that in Case (a3), α 1 is identifiable but α 2 is not identifiable with the proposed method. On the other hand, in the case that the two subpopulations share the same β, α 1 and α 2 are identifiable in Case (a4), but α 1 and α 2 are not identifiable in Cases (a5)-(a6). In addition, the results in Table 1 and 2 show that the proposed pooled case-control studies method produces accurate estimation for the population percentages of cases and controls. All observations here are consistent with the discussions in section 2.1.

Simulated data: comparing three different methods
In the second part, we conduct simulations to examine the efficiency gain in estimating the slope parameters with our proposed method by comparing the variance of three estimators: the maximum likelihood estimator for single case-control study with known f (·) denoted byθ (1) , the maximum likelihood estimator for single case-control study with unknown f (·) (Scott and Wild, 1986) denoted byθ (2) and our proposed estimator denoted byθ (3) . The settings of cases (b1)-(b4) are given in Table 3. For the convenience of comparison, we set the dimension of β to be one in all cases in this subsection. The results of cases (b1)-(b4) summarized in Tables 4-5 are based on 1000 replications.
By comparing the ESEs of β i in Tables 4-5, one can observe thatβ (1) is most efficient andβ (2) is least efficient among the three estimators, which confirms our theory. Significant efficiency gain in estimating the slope parameter with our proposed method are observed when some r k is getting close to 0, as in the second case-control study in Cases (b2) and (b3). Furthermore, as shown in Table 4 Case (b3),β 2 (1) is nearly as efficient asβ 2 (3) and the estimated standard errors ofα (1) 2 are significantly smaller, when the sample sizes of the other two case-control studies are large enough. On the other hand, when the number of case-control studies increases, it can be seen from Cases (b1) and (b4) in Table 5 thatβ i . Lastly, we observe that when the proportion of cases q i in single case-control study sample is close to 1 or 0 and |p i − q i | is large, as in the 1st, 2nd, 4th and 5th case-control studies in Case (b4), combining several case-control studies can balance the proportions between cases and controls, and thus leads to more efficient estimation for β i .

Hepatitis C infection dataset
The proposed method is applied to real-life data from a published study on liver fibrosis and cirrhosis in patients with chronic hepatitis C infection, which is available at https:// archive.ics.uci.edu/ml/datasets/HCV+data (Lichtinghagen et al. 2013). It is known that progressive fibrosis is a major cause of morbidity and mortality in chronic liver disease. Factors having high correlations with fibrosis stages in chronic liver disease are included in this study.
Following the data pre-processing procedure in Hoffmann et al. (2018), 589 patients are selected aged between 23 to 77. They have been divided into 4 categories already: the healthy ones, the Hepatitis, the Fibrosis and the Cirrhosis, with sample size 533, 20, 12 and 24 respectively. Among the 12 predictors (age, sex and 10 other measurands), according to Hoffmann et al. (2018), 6 predictors including ALB, BIL, CHE, GGT, AST, ALT are used in the decision tree for clustering and performs well. Hence, we fit model (2) with the 6 factors. To illustrate the idea of combining case-control studies, we randomly partition the 533 healthy samples into three subsets with size 177, 177 and 179 respectively, and obtain 3 case-control studies on Hepatitis, Fibrosis and Cirrhosis as the mixture of the healthy subgroup (controls) with each diseased groups (cases) respectively. Especially, the casecontrol study on Hepatitis consists of 177 controls and 20 Hepatitis cases; the case-control study on Fibrosis consists of 177 controls and 12 Fibrosis cases; the case-control studies on Cirrhosis consists of 179 controls and 24 Cirrhosis cases.
We apply our combining case-control studies method to fit model (2) with 6 predictors and K = 3. We repeat our proposed method by randomly splitting the healthy group for 50 times. For comparison, the odds ratio method by Prentice and Pyke (1979) is applied to each case-control study. The results are presented in Table 6, from which one can see that each factor contributes to the severity of the disease. By combining the three case-control studies, our proposed estimator gives stable and more efficient estimates compared with the classical odds ratio estimator. Notes: "Par." stands for parameter; "True" means the true value of the parameter; the symbol * means the value is greater than 10 4 . Notes: "Par." stands for parameter; "True" means the true value of the parameter; The symbol * means the value is greater than 10 4 .   Notes: "Par." stands for parameter; "True" means the true value of the parameter.  Notes: "Par." stands for parameter; "True" means the true value of the parameter. Notes: "CC" stands for the proposed method by combining case-control studies; "ESE" represents the estimated standard errors for the parameter estimate; "OR" stands for the odds ratio estimator by Prentice and Pyke (1979).