In studies in medical science or public health, the health outcome we want to predict or explain with independent variables is often binary. For example, the outcome could be whether a patient has a specific disease or not (e.g. diabetic or not), or it could be whether a patient develops a new condition within a certain period of time (e.g. recurrence of cancer within 3 years since surgery). This makes logistic regression one of the most commonly used methods among others in health research.

Logistic regression assumes that the response variable follows the Bernoulli distribution and the probability of “success” is equal to the logistic function of linear combination of independent variables. There is no distributional assumption required for the independent variables in logistic regression. However, when an independent variable has a very skewed distribution to the right, it is common to apply log transformation of such variables in medical science and public health (Ekwaru and Veugelers, 2018; Feng et al., 2014). A typical example of variables with right-skewed distribution includes lab values such as Creatinine, CRP (C-Reactive Protein), or CA19-9 (Carbohydrate antigen 19-9). These variables tend to have values within a “normal” range when the patient is healthy, but can have much higher values, sometimes multiple hundred times the normal value, when the patient is in an abnormal condition. This makes the distribution severely skewed to the right. Although logistic regression does not require independent variables to have symmetric distribution, these variables with skewed distribution tend to explain the response variable better when they are transformed into a more symmetric distribution. This is because extremely high, abnormal lab values work as influential points and thus tend to ruin the overall model fit. Also, it is unlikely that having multiple hundred times the normal lab value really has hundred times the impact that a patient with a normal lab value would have on a logit of “success” probability.

Log-transformation is simply replacing a variable X with log(X). Due to its simplicity and popularity, a lot of statistical software such as SAS and SPSS includes built-in log-transformation. It is a special case of the Box-Cox transformation (Box and Cox, 1964) and tends to make rightly skewed distribution less skewed. The downside of log-transformation is that it can only be applied to positive variables because logarithm is not defined at zero or negatives.

Although lab values are usually nonnegative, often times they have zeros. This is because in many cases, due to technical restrictions, there is a limit of detection. That is, if the amount of the substance of interest goes smaller than a certain value, it is not detected and is recorded as zero. Therefore, to log-transform a lab value, it is wise to check if there are zeros, and if so, preprocessing to handle zeros is necessary before log-transformation.

The most commonly used method to handle zeros before log-transformation is the so called “started logarithm” (Rocke and Durbin-Johnson, 2003). This method is adding a small positive constant before log-transformation. That is, replacing X with log(X + c), where c is a very small, positive constant. For c, 1, or half of the minimum non-zero values, or arbitrary numbers such as 0.01 or 0.001 are often used. There have been studies on how to select optimal value for c in the literature (Ekwaru and Veugelers, 2018; Rocke and Durbin-Johnson, 2001, 2003; Durbin and Rocke, 2004). Also, there have been suggestions for new methods of zero imputation as well, but it is challenging to use such methods in practice due to their complexity (Bellégo et al., 2006).

Park (2023) introduced five zero imputation methods commonly used in practice including “started logarithm” and compared their performances in linear regression setting. However, the impact of use of these methods in logistic regression remains unclear.

In this study, we investigate the effect of zero imputation methods for log-transformation of independent variables in logistic regression setting using a simulation study. The rest of this paper is organized as follows: In Chapter 2, we introduce the logistic regression setting with a right-skewed independent variable. In Chapter 3, the five most common zero imputation methods used in medical science are introduced. Chapter 4 presents the details and the results of the simulation study, and we conclude with recommendations for zero imputation methods in Chapter 5.

We are interested in the situation where the logistic regression model (2.1) explains the relationship between the binary response variable and the continuous, nonnegative independent variable as,

The model above reflects the common situation where X has right-skewed distribution and it explains P(Y = 1) better when X is log-transformed. We assume that log(X) has normal distribution (thus X follows log-normal distribution which is skewed to the right) and we cannot observe the exact value of X if X < d, where d is the limit of detection and its value is unknown. That is, we observe the realization of X̃ instead of X, where

Setting up X̃ this way, we are reflecting the nature of lab values which have skewed distribution and have zeros due to limit of detection.

Let x̃₁, . . . ,x̃_n and y₁, . . . ,y_n be the observed values of X̃ and Y in a sample with size n. Using these observations one would try to fit the model in the following:

However, log(x̃_j) cannot be computed for x̃_j = 0. Therefore, we need to use transformation z_i = f (x̃_i) to replace log(x̃_j) in (2.2). Here, f (x̃_i) needs to be close to log(x̃_j) and should be defined at zero.

We consider five different choices for z_i = f (x̃_i) to replace log(x̃_j) in model equation (2.2) as follows:

3.1. Shift by one

This method adds one to the whole observations, and then takes the logarithm. This can be written as follows:

The rationale behind this method is that it transforms zero to zero. That is, z_i becomes zero when x̃_i = 0. This method only involves arithmetic operations and does not require conditional statements, which makes the use of this method simple and easy. The downside of this method is that when X̃ is distributed too closely to 1, the impact of shifting by one can be too big so that z_i can be very different from log(x_i). In practice, sometimes any other arbitrary number such as 0.01 or 0.001 is added instead of 1, but then it does not transform zero to zero anymore and the rationale is lost.

3.2. Shift by half of the smallest nonzero

In this method, half of the minimum value of nonzero observations is added to x̃_i = 0 before log-transformation. In particular, we use z_i as follows:

This method is different from the ‘shift by one’ method in the sense that the added constant is determined based on the distribution of the observed data. Thus using this method prevents z_i from becoming too far from log(x_i). If we think of min_x̃i>0{x̃₁, . . . ,x̃_n} as an estimate of d, d̂, then this method is equivalent to shifting all the observations by arithmetic mean of zero and d̂ before log-transformation.

3.3. Shift by the square root of the smallest nonzero

We can consider adding the square root of min_x̃i>0{x̃₁, . . . ,x̃_n} instead of halving it, as follows:

Here, we can think of minx˜i>0{x˜1,…,x˜n} as the geometric mean of one and d̂. When min_x̃i>0{x̃₁, . . . ,x̃_n} < 4, this method adds a slightly bigger constant to the observed values than ‘shift by half of the smallest nonzero’ method, but otherwise it adds a smaller constant resulting in z_i being closer to log(x_i).

3.4. Replace zeros with half of the smallest nonzero

This method keeps nonzero observations the same and only replaces zeros with a half of the minimum of nonzero observations. That is, we use transformation as follows:

For positive observations, this method does not change anything before log-transformation, resulting in z_i being pretty close to log(x_i) overall. However, it requires a conditional statement in the analysis and this can be challenging in some software.

3.5. Replace zeros with the square root of the smallest nonzero

Similar to the previous method, this method only replaces zeros with a positive constant and keeps nonzeros as they are. The difference is that the square root of the minimum of nonzero observations is used as the positive constant rather than half of the smallest nonzero. This can be written as follows:

When the minimum of nonzero observations is less than 1, minx˜i>0{x˜1,…,x˜n}>minx˜i>0{x˜1,…,x˜n} holds. This means that zeros can be replaced by a constant that is greater than the minimum of nonzero observations, so that the order of the observations may not be preserved.

To evaluate the performance of the five zero imputation methods outlined in the previous chapter, we carried out a simulation study and report the results here. Here, we outline the steps to generate simulated data. First, to simulate a continuous, right-skewed independent variable, we generated random numbers that follow a log-normal distribution. We considered three different kinds of log-normal distribution: Lognormal(0, 1²), Lognormal(3, 1²), Lognormal(3, 3²), which are similar to the distributions of commonly used lab values CEA (carcinoembryonic antigen), CRP, and CA19-9, respectively. To help grasp the overall distribution, mean, median, 1, 3, 5, 10, and 99.5 percentile of X is listed in Table 1.

Second, the lowest 1, 3, 5, or 10% of the generated values were replaced by zero to reflect the measurement error due to limit of detection.

Third, the values of the response variable were generated from the Bernoulii distribution with the success probability calculated by equation (2.2) based on the simulated data of the independent variable. We set β₁ = 1 for all scenarios. For β₀, we set the values such that the “success” probability becomes 0.5 or 0.2 when X = E(X). That is, we simulated data for two types of outcome distribution: P(Y = 1|X = E(X)) = 0.5 and P(Y = 1|X = E(X)) = 0.2. The first type represents the case where the distribution of the binary response variable is balanced, and the second type shows the case where the response variable has an unbalanced distribution.

In summary, we considered 6 different scenarios. Scenario 1, 2, and 3 are for the cases where P(Y = 1|X = E(X)) = 0.5 and Scenario 4, 5, and 6 are for the cases where P(Y = 1|X = E(X)) = 0.2. In Scenario 1 and 4, the independent variable X follows Lognormal(0, 1²). Scenario 2 and 5 have Lognormal(3, 1²), and Scenario 3 and 6 have Lognormal(3, 3²) for the distribution of X.

We generated 1000 sets of samples, each with a sample size of 200. We repeated the same experiment with sample size of 100 and 300, and the results were similar, but are not reported here. For measure of performance, we calculated bias of MLE estimates of β₀, β₁, and exp(β₁) which is the OR (odds ratio) of X. We also generated a separate test set with a size of 1000 for each scenario so that we could evaluate the prediction performance of the resulting logistic regression model. The c-index of the model was calculated and reported. In addition, we calculated calibration in the large and calibration slope based on the test set so that we can understand the effect of the zero imputation method on calibration of the resulting logistic regression model. Calibration means the agreement between the observed outcomes and the predicted response. In the logistic regression context, calibration in the large stands for the mean difference between estimated vs. actual logit of the success probability, and the calibration slope refers to the slope of the recalibration model: Logit(p) = α + β · logit(p̂) where p = Pr(Y = 1|X) (Steyerberg, 2019). If the fitted model has perfect calibration, the calibration in the large would be 0 and the calibration slope would be 1. The performance measures were calculated for each of the 1000 simulated data sets, and we report their mean and standard deviation in Tables 2–7 and Figures 1–10.

First, let’s take a look at the scenario 1 to 3 where the response variable has balanced distribution. Figure 1 shows that the bias in the intercept term tends to be negative. This is because, all five zero imputation methods tend to make the value of independent variable greater than or equal to the originally observed value. That is, z_i ≥ log(x̃_i) holds for all cases except for rare cases where minx˜i>0{x˜1,…,x˜n}<minx˜i>0{x˜1,…,x˜n}<zi<1 holds and ‘replace zeros with square root of the smallest nonzero’ method is used. Hence, β₀+β₁z_i would be greater than β₀+β₁ log(x_i) for most cases when β₁ > 0. To compensate the increase in the value of β₀ + β₁X in the logistic regression model, the estimate of β̂₀ tends to become lower than β₀. Bias of β̂₀ is the greatest (in terms of absolute value) when the ‘shift by 1’ method was used in scenario 1 because using this method in scenario 1 the difference between the observed value log(x̃_i) and the imputed value z_i is the greatest.

Bias in the slope is presented in Figure 2. With all five zero imputation methods, the range of z_i’s becomes narrower than that of the original observations log(x_i)’s. Therefore, β̂₁ tends to be greater than β₁ to make β₀ + β₁z_i cover the range of β₀ + β₁ log(x_i), which is why the bias of β̂₁ tends to be positive. Shrinkage in the range is the greatest in scenario 1 using ‘shift by 1’ method, which is reflected in the magnitude of the bias of β̂₁. Bias in the OR is directly related with bias in the slope because OR = exp(β₁) holds. Figure 3 shows that the bias of OR can be quite dramatic because of exponentiation.

The c-indices in the training set and testing set were almost the same for all five methods because the different zero imputation methods do not change the order of observations of the independent variable most of the time.

Figure 4 shows calibration in the large is all positive, which means that the success probability is underestimated overall using z_i. Figure 5 shows that the calibration slopes are less than 1, which implies that the estimated success probabilities are too extreme, that is, they are close to either zero or one. This seems to be due to the fact that the bias of β̂₁ is positive. That is, because β₁ is overestimated, effect of independent variable is exaggerated.

Now, we consider the results of scenario 4 to 6, where the response variable has an unbalanced distribution. From Figures 6, 7, and 8, we can see similar patterns to those in scenario 1 to 3, where the bias of the intercept term tends to be negative and the bias of the slope is more likely to be positive.

Figure 9 shows calibration in the large is all positive, which is also similar to scenario 1 to 3, except that the magnitude of calibration in the large is much bigger in scenario 4. The calibration slope was less than 1 in scenario 5 and 6, but tends to be greater than 1 in scenario 4. These results imply that when the success probability is low and the range of the independent variable is narrow, the success probability tends to be greatly underestimated and the effect of the independent variable is also underestimated. Overall, shift by 1 method remains to be the worst method, similarly to scenario 1 to 3.

We investigated the impact of zero imputation methods for nonnegative independent variables with skewed distribution in logistic regression. We evaluated and compared the performance of five zero imputation methods commonly used in clinical research using a simulation study. Overall, bias in the estimated coefficients tends to be high when the range of independent variable is narrow, especially with the ‘shift by 1’ method. Bias can be dramatically high in OR estimate and thus, caution is necessary when selecting zero imputation method especially when the range of the observed values of independent variable is rather narrow. The effect of the zero imputation method in the discrimination of the resulting logistic regression model is minimal, but calibration suffers somehow, although there is not so much of a difference among the zero imputation methods. Overall, it is recommended to avoid the shift by 1 method because it can yield too high a bias in the parameter estimation despite its ease of use. Shift by half of the smallest nonzero, replace zeros with half of the smallest nonzero, and replace zeros with the square root of the smallest nonzero methods showed comparable and stable performance. In the future, investigation of the effect of the zero imputation methods in the Cox proportional hazards model, which is another frequently used method in clinical research, is warranted. Another possible research subject includes investigation of the effect of the zero imputation methods on agreement studies, where a new technique to measure a quantity is compared to the gold standard.

This research was supported by the Korea National Open University Research Fund.