Should a Normal Imputation Model be Modified to Impute Skewed Variables?

Abstract

Researchers often impute continuous variables under an assumption of normality–yet many incomplete variables are skewed. We find that imputing skewed continuous variables under a normal model can lead to bias. The bias is usually mild for popular estimands such as means, standard deviations, and linear regression coefficients, but the bias can be severe for more shape-dependent estimands such as percentiles or the coefficient of skewness. We test several methods for adapting a normal imputation model to accommodate skewness, including methods that transform, truncate, or censor (round) normally imputed values as well as methods that impute values from a quadratic or truncated regression. None of these modifications reliably reduces the biases of the normal model, and some modifications can make the biases much worse. We conclude that, if one has to impute a skewed variable under a normal model, it is usually safest to do so without modifications–unless you are more interested in estimating percentiles and shape than in estimating means, variances, and regressions. In the conclusion, we briefly discuss promising developments in the area of continuous imputation models that do not assume normality.

Keywords

missing data missing values incomplete data regression transformation normalization multiple imputation imputation

Introduction

Imputation is an increasingly popular method for handling data with missing values. When using imputation, analysts fill in missing values with random draws from an imputation model and then fit the imputed data to an analysis model. In multiple imputation (MI), the process of imputation and analysis is repeated several times and the results of the several analyses are combined (Rubin 1987; Allison 2002; Kenward and Carpenter 2007).

In an ideal world, the imputation model would perfectly represent the distribution of the data. But perfect fidelity can be very difficult to achieve and in practice it is often unnecessary. All that is necessary is that the imputation model preserves those aspects of the distribution that are relevant to the analysis model. For example, if the analyst plans only to estimates means, standard deviations, and the parameters of a linear regression model, then the imputation model needs only to preserve the means, variances, and covariances among the variables that will be analyzed.

When the goals of analysis are limited in this way, as they often are, a crude imputation model can yield usable results. For example, in some settings, it can be acceptable to impute dummy variables, squared terms, and interactions as though they were normal when conditioned on other variables (Horton, Lipsitz, and Parzen 2003; Allison 2005; Bernaards, Belin, and Schafer 2007; von Hippel 2009). The resulting imputed values will look implausible if inspected closely but that often has little effect on the analytic results. In fact, attempts to edit the imputed values to improve their plausibility can introduce bias by changing the variables’ means, variance, and covariances (Horton et al. 2003; Allison 2005; Bernaards et al. 2007; von Hippel 2009).

The point of imputation is not that the imputed values should look like observed values. The point is that the imputed variable should act like the observed variable when used in analysis.

This article considers methods for imputing a different class of nonnormal variables—namely, variables with skew. Like other nonnormal variables, skewed variables are often imputed as though they were conditionally normal. The word conditionally is important here; later, we will encounter situations where a variable is skewed and yet the residuals are approximately normal when the skewed variable is conditioned on other variables. In that situation, a conditionally normal imputation model may be perfectly specified even though the imputed variable is skewed. In many other situations, though, imputing a skewed variable from a normal model entails some degree of misspecification.

In calculations and simulations, we find that conditionally normal imputation of a skewed variable can often produce acceptable estimates if the quantities being estimated are means, variances, and regressions. The estimates do have biases under some circumstances but the biases are typically small. However, the biases grow much larger if we estimate quantities that depend more strongly on distributional shape—quantities such as percentiles or the coefficient of skewness.

To increase the skew of a normally imputed variable, popular references recommend modifications to the normal imputation model. Modifications include rounding (censoring) the imputed values, truncating the imputed values, imputing from a truncated regression model, or transforming the incomplete variable to better approximate normality. We evaluate all these methods as well as a new method that adds quadratic terms to the imputation model.

In univariate data, we find that such modifications can work well if very carefully applied. But in bivariate and trivariate data, we find that even careful modification of the normal imputation model is problematic. At best, modifications reduce the biases only a little; at worst, modifications make the biases much worse. We conclude that, if you have to impute a skewed variable as though it were normal, it is safest to do so without modification.

In the conclusion, we discuss nonnormal imputation models that are currently in development. Early evaluations of these methods look promising and we hope that they will soon become more widely available.

Our presentation proceeds in order of complexity, progressing from univariate to bivariate to multivariate data.

Univariate Data

Imputation is rarely useful in univariate data since, if univariate values are missing at random (MAR), the observed values are a random sample from the population and can be analyzed just as though the data were complete. Yet the simplicity of the univariate setting can help clarify the advantages and disadvantages of different imputation methods.

Suppose we have a simple random sample of n values from a nonnormal variable X with finite mean μ and finite variance σ². A randomly selected n_mis of the X values are missing, so that n_obs = n–n_mis values are observed. We assume that values are MAR (Rubin 1976; Heitjan and Basu 1996), which in the univariate setting means that the n_obs observed X values are a random sample from the n sampled cases and therefore a random sample from the population.

How well can normal imputation, with or without modification, impute this nonnormal variable?

Fully Normal (FN) Imputation

The first and simplest technique is fully normal (FN) imputation (Rubin and Schenker 1986), under which a nonnormal variable X is imputed as though it were normal.

Since the observed values are a random sample from the population, we can obtain consistent estimates ${\hat{μ}}_{o b s}, {\hat{σ}}_{o b s}^{2}$ from the observed values alone, then use those estimates to impute the missing values with random draws from a normal distribution $N ({\hat{μ}}_{o b s}, {\hat{σ}}_{o b s}^{2})$ . We then analyze the mix of imputed and observed values as though it were complete.

The consistency of the observed-value estimators ${\hat{μ}}_{o b s}, {\hat{σ}}_{o b s}^{2}$ does not depend on X being normal. The most obvious observed-value estimator are the mean ${\overset{ˉ}{X}}_{n_{o b s}}$ and variance $s_{n_{o b s} - 1}^{2}$ of the observed values, which are minimum variance unbiased (MVU) estimators whether X is normal or not (Halmos 1946). And the estimators that are commonly used in MI are obtained by supplementing ${\overset{ˉ}{X}}_{n_{o b s}}$ and $s_{n_{o b s} - 1}^{2}$ with random variation to obtain Bayesian posterior draw (PD) estimators ${\hat{σ}}_{o b s, P D}^{2}$ and ${\hat{μ}}_{o b s, P D}$ —

{\hat{σ}}_{o b s, P D}^{2} = s_{n_{o b s} - 1}^{2} \frac{n_{o b s - 1}}{U}, w h e r e U ~ χ_{n_{o b s} - 1 + ν_{p r i o r}}^{2}

{\hat{μ}}_{o b s, P D}^{} ~ N ({\overset{ˉ}{X}}_{n_{o b s}}, \frac{{\hat{σ}}_{o b s, P D}^{2}}{n_{o b s}})

—which simulate random draws from the posterior density of

μ, σ^{2}

(Rubin and Schenker 1986). Note that

ν_{p r i o r}

in equation (1) is the degrees of freedom in the Bayesian prior (von Hippel 2012).

The PD estimators, like the MVU estimators, are consistent.¹ However, the PD estimators are not particularly efficient, and ${\hat{σ}}_{o b s, P D}^{2}$ can be biased in small samples (von Hippel 2012).

Small-sample properties of the MVU and PD estimators are discussed elsewhere (von Hippel 2012). We can avoid some detail by assuming that we have a very large sample where ${\hat{μ}}_{o b s}, {\hat{σ}}_{o b s}^{2}$ are practically indistinguishable from $μ, σ^{2}$ . Under that assumption, we are effectively imputing missing values from $N (μ, σ^{2})$ , and any biases we encounter are asymptotic—that is, any biases we encounter cannot be eliminated by increasing the sample size.

Figure 1 illustrates the large-sample properties of FN imputation in a setting where the observed data are skewed. Here the observed values come from a standard exponential density $X_{o b s} ~ E x p (μ)$ , with mean $μ = 1$ and variance $μ^{2} = 1$ , while the imputed values $X_{i m p} ~ N (μ, μ^{2})$ come from a normal density with the same mean and variance as the observed exponential variable.

Figure 1.

Distribution of a standard exponential observed variable and a normal imputed variable that has the same mean and variance.

Although the observed and imputed variables have the same mean and variance, they do not have the same distributional shape, and that implies that, in a mix of observed and imputed values, some quantities will be estimated with bias. For example, estimates of skewness will be biased toward zero since the imputed values have no skew. And the first sixteen percentiles, at least, will be negatively biased since 16 percent of the imputed distribution is negative while the observed distribution is strictly nonnegative. The lower panel of Figure 1 examines all the percentiles by comparing the cumulative distribution function (CDF) of the observed and imputed variables. Percentiles 24–89.5 are positively biased since the imputed CDF is right of the observed CDF, but all the other percentiles are negatively biased since the imputed CDF is left of the observed CDF. In some distributions, FN imputation of nonnormal variables causes more bias in the extreme percentiles than in percentiles near the median (Demirtas, Freels, and Yucel 2008) but in this skewed example the 50th percentile has substantial bias while the 90th percentile has almost none.

In sum, when nonnormal data are imputed under a normal model, estimates of the mean and variance will be consistent, but there can be considerable bias in estimating the shape and percentiles. Simulations confirm that FN imputation yields consistent estimates for $μ$ and $σ^{2}$ under a variety of nonnormal distributions (Rubin and Schenker 1986; He and Raghunathan 2006; Demirtas and Hedeker 2008) but FN imputation yields biased estimates of many percentiles (He and Raghunathan 2006; Demirtas et al. 2008; Demirtas and Hedeker 2008).

Imputing Within Bounds: Censoring and Truncation

What can make the observed and imputed distributions more similar? A popular approach is to bound normally imputed values within a plausible range. For example, in imputing a positively skewed variable like body weight, we might require the imputed values to be larger than the smallest observed body weight, or not so small as to be “biologically implausible” (Centers for Disease Control and Prevention 2011). Options for bounding imputed variables are available in most popular imputation software, including IVEware, the MI procedure in SAS 9.2, the mi impute command in Stata 12, and the Missing Values Analysis package in SPSS 16.0.

There are two general approaches to bounding: censoring and truncation. For example, in Figure 1, where we imputed a standard exponential variable as though it were normal, we could have used censoring or truncation to ensure that the normally imputed values were nonnegative. To do this, we would first generate normal imputations $X_{i m p} ~ N (μ_{p r e}, σ_{p r e}^{2})$ where $μ_{p r e} {, σ}_{p r e}^{2}$ are the mean and variance of the imputed variable before bounds have been imposed. We would then round negative imputed values up to zero (censoring), or reimpute negative values repeatedly until a positive value occurs (truncation). Truncation results in a distribution of imputed values that is strictly positive, while censoring results in a distribution that is strictly positive except for a point mass at $X_{i m p, c e n s o r} = 0$ .

Bounding imputed values changes the mean and variance of the imputed variable. After bounding, the mean of the imputed variable is higher than $μ_{p r e}$ because negative values have been shifted upward, and the variance of the imputed variable is smaller than $σ_{p r e}^{2}$ because out-of-bounds values have been shifted inward. More specifically, if an imputed normal variable $X_{i m p} ~ N (μ_{p r e} {, σ}_{p r e}^{2})$ is truncated to the left of c, the truncated variable has mean and variance

\begin{aligned} E (X_{i m p, t r u n c}) = μ_{p r e} + λ_{c} σ_{p r e} \\ V (X_{i m p, t r u n c}) = (1 - δ_{c}) σ_{p r e}^{2} . \end{aligned}

And if the same imputed normal variable is censored to the left of c, the mean and variance are

E (X_{i m p, c e n s o r}) = π_{c} c + (1 - π_{c}) E (X_{i m p, t r u n c}) .

V (X_{i m p, c e n s o r}) = (1 - π_{c}) (V (X_{i m p, t r u n c}) + {(λ_{c} - Z_{c})}^{2} π_{c} σ_{p r e}^{2}) .

Here

Z_{c} = (c - μ_{p r e}) / σ_{p r e}

is the standardized boundary c;

ϕ (Z_{c})

and

π_{c} = Φ (Z_{c})

are the standard normal PDF and CDF evaluated at

Z_{c}

;

λ_{c} = ϕ (Z_{c}) / (1 - Φ (Z_{c}))

is the Mills ratio; and

δ_{c} = λ_{c} (λ_{c} - Z_{c})

(e.g., Johnson, Kotz, and Balakrishnan 1994, section 10.3; Greene 1999:907; Jöreskog 2002).

Bias can occur if we fail to anticipate the effect of bounding on the mean and variance of the imputed variable. To return to our running example, if the observed variable is standard exponential $X_{o b s} ~ E x p (μ)$ , with mean $μ = 1$ and variance $μ^{2} = 1$ , then we might naively impute a normal variable with the same mean and variance $μ_{p r e} = σ_{p r e}^{2} = 1$ . But if we censor this normal variable to the left of $c = 0$ , the mean of the censored variable will be $E (X_{i m p, c e n s o r}) \approx 1.08$ and the variance will be $V (X_{i m p, c e n s o r}) \approx 0.75$ —that is, the mean will have a bias of 8 percent and the variance will have a bias of −25 percent. The biases will be even larger if we truncate the normal imputed variable to the left of $c = 0$ ; then the mean of the truncated variable will be $E (X_{i m p, t r u n c}) \approx 1.29$ , which is biased by 29 percent, and the variance will be $V (X_{i m p, t r u n c}) \approx 0.63$ , which is biased by –37 percent.

These biases can be avoided, or at least reduced, if we anticipate the effect of censoring on the mean and variance of the imputed variable. By inverting equations (2) and (3), we can choose $μ_{p r e}, σ_{p r e}^{2}$ so that, after bounds are imposed, the mean and variance of the bound imputations will be the same as those of $X_{o b s}$ , or at least close. Let us return to our example where the observed variable $X_{o b s} ~ E x p (μ)$ is standard exponential with mean $μ = 1$ and variance $μ^{2} = 1$ . If we want to get the same mean and variance from a normal variable that is censored on the left at c = 0, we can do so by choosing $μ_{p r e} \approx .785$ and $σ_{p r e}^{} \approx 1.29$ ; then, after censoring, we will have a censored variable with a mean and variance of $E (X_{i m p, c e n s o r}) = V (X_{i m p, c e n s o r}) = 1$ .

Similarly, if we want to get a similar mean and variance from a normal variable that is truncated on the left at c = 0, we can get close by choosing $μ_{p r e} = - 6.28$ and $σ_{p r e}^{} = 3.02$ . Then after truncation we will have a truncated variable with a mean of $E (X_{i m p, t r u n c}) = 1.1$ and a variance of $V (X_{i m p, t r u n c}) = 1$ .

Notice that, although the censored normal imputations have the same mean as the observed variable, the truncated normal imputations do not. In fact, with c = 0, it is impossible to get the mean and variance of a censored normal variable to equal $E (X_{i m p, t r u n c}) = V (X_{i m p, t r u n c}) = 1$ exactly. In fact, as we approach $E (X_{i m p, t r u n c}) = V (X_{i m p, t r u n c}) = 1$ , $μ_{p r e} {/ σ}_{p r e}^{}$ approaches –∞. In the area near $E (X_{i m p, t r u n c}) = V (X_{i m p, t r u n c}) = 1$ it is very hard to estimate the parameters $μ_{p r e} {, σ}_{p r e}^{}$ with any precision.

We have just seen that, even in a simple univariate setting, a truncated normal model can have trouble matching the moments of a variable that is not in fact truncated normal. When the truncated normal model is used for this purpose, the parameter estimates can be sensitive or even infinite. Later, we will encounter similar problems when we use the truncated model in a multivariate setting.

Figure 2 illustrates our attempts to match the distribution of a standard exponential variable with a truncated or censored normal variable. The figure illustrates a naïve approach, where the mean and variance of the imputed variable are unbiased before censoring or truncation, but biased after. The figure also illustrates less-biased approaches, where the effect of truncation or censoring is anticipated before the imputations are drawn. Notice that, if the effects of truncation are anticipated, the truncated normal variable achieves a close match to the exponential distribution (lower left)—notwithstanding some bias in the mean.

Figure 2.

Biased and less-biased ways of imputing a censored or truncated normal variable to match the distribution of an observed exponential variable.

Transformation

Instead of bounding imputed X values, some experts recommend transforming an incomplete skewed X variable to better approximate normality (Schafer and Olsen 1998:550; Schafer and Graham 2002:167; Allison 2002:39; Raghunathan et al. 2001:82-83). Following this advice, we would subject the observed values $X_{o b s}$ to a skew-reducing transformation such as a log or root, which yields a transformed variable $X_{o b s, t} = t (X_{o b s})$ with a mean and variance of $μ_{t}, σ_{t}^{2}$ . We then apply FN imputation to the transformed variable to obtain normal imputed values $X_{i m p, t} ~ N (μ_{t}, σ_{t}^{2})$ . Finally, we invert the transformation to recover the original observed values $X_{o b s} = t^{- 1} (X_{o b s, t})$ and put the imputed values $X_{i m p} = t^{- 1} (X_{i m p, t})$ on a comparable scale. SAS automates this process in the MI procedure, where the TRANSFORM command offers a choice among a variety of nonlinear transformations including the log, power, Box-Cox (1964), logit, and exponential transformations (SAS Institute 2001).

In the univariate setting, transformation works well if the transformed variable is in fact normal or close to it. For example, if we apply a log transformation to a lognormal variable, the transformed variable is exactly normal.

But transformation can yield substantial bias if the transformed variable is not close to normal. To return to our earlier example, if we sample an exponential variable $X_{o b s} ~ E x p (μ)$ and transform it with a log, the logged variable $X_{o b s, t} = l n (X_{o b s})$ is far from normal (Hawkins and Wixley 1986). If we impute a normal variable $X_{i m p, t} ~ N (μ_{t}, σ_{t}^{2})$ on the log scale, inverting the transformation will yield an imputed $X_{i m p} = e x p (X_{i m p, t})$ whose mean and variance are far from the mean and variance of the original exponential variable (He and Raghunathan 2006).

The challenge of transformation is that skew-reducing transformations are nonlinear, and nonlinear transformations do complicated things to the mean and variance. Although on the transformed scale the observed variable $X_{o b s, t}$ and the imputed variable $X_{i m p, t}$ have approximately the same mean and variance $μ_{t}, σ_{t}^{2}$ , after the transformation is inverted there is no guarantee that $X_{o b s}$ and $X_{i m p}$ , unless they have the same distribution, will also have the same mean and variance $μ, σ^{2}$ .

The transformation method can yield better results if the transformation is carefully chosen. In imputing an exponential variable $X_{o b s} ~ E x p (μ)$ , for example, it is not hard to see why the logarithm is a poor choice. The log of an exponential variable is undefined at $X_{o b s} = 0$ , and although a value of exactly zero will never occur in an exponential sample, values close to zero occur frequently, and the logs of these near-zero values will be large negative numbers. So the log of a standard exponential sample has a long left tail and is far from normal.²

Figure 3 illustrates two transformations that do a better job of normalizing an exponential variable $X_{o b s}$ : the square root transformation $\sqrt{X_{o b s}}$ and the fourth root transformation $\sqrt[4]{X_{o b s}}$ . Unlike the log, both the square root and the fourth root are defined at $X_{o b s} = 0$ . The square root transformation is common in data analysis while the fourth root is less common but comes as close to normalizing the exponential distribution as nearly any transformation can (Hawkins and Wixley 1986).

Figure 3.

Normal imputation of an observed exponential variable under square-root and fourth-root transformation.

Both the square-root and fourth-root transformations yield an imputed variable whose mean is unbiased, or nearly unbiased, when compared to the mean of an exponential variable. The variance of the imputed variable is biased, though less biased for the fourth root than for the square root.

More specifically,³ if we subject an observed standard exponential variable $X_{o b s} ~ E x p (1)$ to a square root transformation, the result is a transformed variable $\sqrt{X_{o b s}} = X_{o b s, t} ~ R a y (1 / \sqrt{2})$ that has a Rayleigh distribution with mean $μ_{t} = \sqrt{π} / 2$ and variance $σ_{t}^{2} = 1 - π / 4$ . On the transformed scale, we can generate a normal imputed variable $X_{i m p, t}$ with the same mean $μ_{t}$ and variance $σ_{t}^{2}$ . Finally, we invert the transformation, squaring all the values to recover the original observed variable $X_{o b s} = X_{o b s, t}^{2} ~ E x p (1)$ and produce an imputed variable $X_{i m p} = X_{i m p, t}^{2}$ that, as the square of a normal variable, has a scaled noncentral chi-square distribution, with a mean of $E (X_{i m p}) = 1$ and a variance of $V (X_{i m p}) = 2 - π^{2} / 8 \approx .77$ . The mean of $X_{i m p}$ is unbiased, but the variance of $X_{i m p}$ , when compared to the variance $V (X_{o b s}) = 1$ of the observed variable, is negatively biased by 23 percent.

We get better results with a fourth-root transformation, although there is still a small bias. If the observed variable $X_{o b s} ~ E x p (1)$ has a standard exponential distribution, then the fourth-root transformed variable $X_{o b s, t} = \sqrt[4]{X_{o b s}}$ is very close to normal with a mean of $μ_{t} \approx .91$ and a variance of $σ_{t}^{2} \approx .065$ . If we impute a normal variable on the transformed scale $X_{i m p, t} ~ N (μ_{t}, σ_{t}^{2})$ with the same mean and variance as $X_{o b s, t}$ , then when we invert the transformation the imputed variable $X_{i m p} = X_{i m p, t}^{4}$ has a mean of 1.006 and a variance of 1.135. The mean of $X_{i m p}$ has less than 1 percent bias while the variance of $X_{i m p}$ is positively biased by 13.5 percent. The shape of the imputed distribution closely matches the observed variable except in the range very close to X = 0 (Figure 3).

In short, in univariate data, the use of normalizing transformation to impute skewed variables can sometimes yield good results if the transformation is carefully chosen. However, we still cannot recommend transformation since it has disadvantages in bivariate settings, which we will discuss next.

Bivariate Data

Bivariate data present new challenges for the imputation of skewed variables. In bivariate data, the imputation model should preserve not just the marginal distribution of the skewed variable; the imputation must also preserve the relationship between the skewed variable and other variables. This can be difficult.

To set the stage, suppose that we have, again, a standard exponential variable $X ~ E x p (1)$ , but now suppose that we also have a second variable Y that fits a linear regression on X with normal residuals:

Y = α_{Y | X} + β_{Y | X} X + e_{Y | X}, where e_{Y | X} ~ N (0, σ_{Y | X}^{2}) .

Figure 4 displays n = 500 simulated observations on (X, Y) with

α_{Y | X} = β_{Y | X} = 1

and

R_{Y | X}^{2} = .7

. The figure shows a scatterplot for Y regressed on X with an ordinary least squares (OLS) line (left) and a scatterplot for X regressed on Y (right) with an OLS line and an OLS quadratic.

Figure 4.

X is standard exponential and Y is conditionally normal. The left panel shows a scatterplot for the regression of Y on X, with an OLS line. The right panel shows a scatterplot for the regression of X on Y, with an OLS line and an OLS quadratic.

In evaluating imputation methods, it will be important to realize that the regression of Y on X satisfies the OLS assumptions, but the regression of X on Y does not. In particular, the regression of X on Y violates the OLS assumption that the conditional expectation of the residuals $E (e_{X | Y} | Y)$ is zero for all Y. To the contrary, $E (e_{X | Y} | Y) > 0$ when Y is small and when Y is large; in the scatterplot this is evident from the fact that all the X values are above the OLS line when Y is small or large. The residuals are nonnormal as well, but this is less important since residual normality is not required for OLS estimates to be consistent (Wooldridge 2001).

In fact, the regression of X on Y is not even linear. When X is nonnormal, a linear regression of Y on X does not necessarily imply a linear regression of X on Y. As Figure 4 show, the fit of the regression of X on Y can be improved, though not perfected, by adding a quadratic term $Y^{2}$ to the regression. With the quadratic term added, the regression comes closer to satisfying the OLS assumption that $E (e_{X | Y} | Y) = 0$ .

Although the OLS estimates are not perfect fits for the regression of X on Y, over most of the range of X the linear fit is a decent approximation, and the quadratic fit is very good. We will find that that OLS estimates, though approximate, can be often be serviceable when we impute X conditionally on Y.

Having looked at complete data, we now suppose that some values are MAR. In the bivariate setting, the MAR assumption means the probability that a value is missing depends only on values that are observed. That is, the probability that Y is missing cannot depend on Y, but can depend on X in cases where X is observed. Likewise, the probability that X is missing can depend on observed Y values in cases where Y is observed (Rubin 1976; Heitjan and Basu 1996).

We can impute missing values using several different methods.

Linear Regression Imputation

The simplest parametric imputation model assumes that the missing values fit a linear regression with normal residuals. We call this linear regression imputation.

To understand the technique, it is helpful to start with a situation where the assumptions of the imputation model are met. Imagine that X is complete and Y is MAR. Then the missing values fit the linear regression in equation (4), and the MAR assumption means that the probability Y is missing is independent of $e_{Y | X}$ . In other words, the distribution of $e_{Y | X}$ is the same in the cases with Y observed as in the cases with Y missing, and therefore we can obtain consistent estimates of the regression parameters ${\overset{}{\hat{α}}}_{Y obs | X}, {\overset{}{\hat{β}}}_{Y obs | X}, {\overset{}{σ^{2}}}_{Y obs | X}^{}$ by simply regressing the Y on X in the cases with Y observed (Anderson 1957; Little 1992; von Hippel 2007). Then to impute missing Y values, we can simulate draws from the regression model in equation (4), with the observed Y estimates substituted for the parameters:

Y_{i m p} = {\hat{α}}_{Y o b s | X} + {\hat{β}}_{Y o b s | X} Y + e_{Y | X, i m p}, where e_{Y | X, i m p} ~ N (0, {\hat{σ}}_{Y o b s | X}^{2}) .

An important detail is what observed Y estimates should be used to run the imputation model. The simplest choice is the OLS estimates:

{\hat{α}}_{Y o b s | X, O L S}, {\hat{β}}_{Y o b s | X, O L S}, {\hat{σ}}_{Y o b s | X, O L S}^{2}

. And the estimates that are usually used in imputation models are obtained by adding random variation to the OLS estimates to obtain Bayesian posterior draw (PD) estimates:

{\hat{σ}}_{Y | X, P D}^{2} = {\hat{σ}}_{Y obs | X, OLS}^{2} \frac{n_{Y obs} - 2}{U}, where U ~ χ_{n_{Y obs} - 2 + ν_{prior}}^{2}

[\begin{matrix} {\hat{α}}_{Y | X, P D} \\ {\hat{β}}_{Y | X, P D} \end{matrix}] ~ N ([\begin{matrix} {\hat{α}}_{Y o b s | X, O L S} \\ {\hat{β}}_{Y o b s | X, O L S} \end{matrix}], \hat{V} [\begin{matrix} {\hat{α}}_{Y o b s | X, O L S} \\ {\hat{β}}_{Y o b s | X, O L S} \end{matrix}]),

where

n_{Y o b s}

is the number of cases with Y observed,

ν_{p r i o r}

is the prior degrees of freedom and

\hat{V} [{\hat{α}}_{Y o b s | X, O L S}, {\hat{β}}_{Y o b s | X, O L S}]

is the usual covariance matrix of the OLS estimates, with the estimate

{\hat{σ}}_{Y | X, P D}^{2}

substituted for the parameter

σ_{Y | X}^{2}

(Little and Rubin 1989; Kim 2004).

Like the OLS estimates, the PD estimates are consistent,⁴ and if Y is regressed on X in the imputed data, the resulting regression estimates are consistent as well, although ${\hat{σ}}_{Y | X, P D}^{2}$ can be biased in small samples (Kim 2004).

In sum, if Y is the incomplete variable, linear regression imputation yields consistent estimates.

If X is the incomplete variable, however, linear regression imputation can be inconsistent. With X incomplete, linear regression imputation assumes that X fits the following equation:

X = α_{X | Y} + β_{X | Y} Y + e_{X | Y}, w h e r e e_{X | Y} ~ N (0, σ_{X | Y}^{2}) .

But in discussing Figure 4 we saw that this equation is incorrect; instead,

e_{X | Y}

is nonnormal, the conditional expectation of

e_{X | Y}

is positive for small and large values of Y, and the conditional expectation of

X

is not a linear function of Y. Therefore, the parameters

α_{X | Y}

and

β_{X | Y}

, to the degree that they are meaningful, cannot be estimated consistently using OLS. And PD estimates that are based on the OLS estimates will be inconsistent as well.

Despite its inconsistency for incomplete X, linear regression imputation is a convenient approximation, and we will find out later that its biases are fairly small.

Quadratic Regression Imputation

Quadratic regression imputation is an attempt to improve on linear regression imputation by adding a squared term $Y^{2}$ to the regression equation for X. As we saw in Figure 4, a quadratic term can improve the fit of the regression. We should recognize, however, that the quadratic model is not a perfect specification and will not fit all data as well as it fits the data in Figure 4. A misspecified quadratic model is vulnerable to extrapolation error. That is, the curve that best fits the observed X values may not be the best fit for the missing X values if the observed and missing X values are in different parts of the distribution.

Imputing Within Bounds: Censoring, Truncation, and Truncated Regression

As in the univariate case, in the bivariate case, we can increase the skew of a normally imputed X by keeping the imputed values within bounds.

The simplest way to bound imputed values is censoring. We impute X values under a linear regression imputation model, and then round (censor) out-of-bounds X values up to the boundary value c. In our running example, where we impute an exponential variable as though it were conditionally normal, we would round negative imputed X values up to zero, which is the lower bound of the exponential distribution.

An alternative approach is truncation. Under a simple truncation algorithm, we impute X values under a linear regression imputation model, then reject any out-of-bounds (e.g., negative) X values and reimpute them until an in-bounds (e.g., positive) X value randomly occurs. This is the approach used by the MI procedure in SAS 9.2, which will terminate if the wait for an in-bounds value is too long. A faster approach is to impute from a truncated distribution using sampling importance resampling (Robert 1995; Raghunathan et al. 2001).

As in the univariate setting, in bivariate data censoring and truncation can improve the shape of the imputed distribution, but censoring and truncation also bias the mean and variance of X. And in the bivariate setting, censoring and truncation do not solve the problem that X has been imputed from a misspecified regression model whose parameters are inconsistently estimated. So bounding linear regression imputations has several sources of bias: bias from bounding imputed values, bias from a misspecified regression model, and bias from inconsistent parameter estimates. If we are lucky, these biases will offset each other; if we are unlucky, they will reinforce one another.

A more sophisticated approach is to impute values with a truncated regression model (Goldberger 1981), which is implemented by IVEware’s BOUNDS statement and by Stata’s mi impute truncreg command (Raghunathan, Solenberger, and Van Hoewyk 2002; Stata Corp. 2011). Under truncated regression imputation, we do not just truncate the imputed X values. We also estimate the regression parameters under the assumption that the observed X values come from a conditionally normal distribution that has been truncated in the same way. This is akin to the univariate approach that we discussed earlier, where the effect of truncation was anticipated in estimating the model parameters.

Like the normal regression model, the truncated regression model makes an assumption about the conditional distribution of X, and the assumption is not correct for many skewed variables. In Figure 4, for example, X follows an exponential distribution, not a truncated conditional normal distribution. The truncated regression model is only an approximation, and the approximation is not necessarily better than the approximation offered by a model without truncation. In fact, our discussion of univariate data showed that, when the assumptions of a truncated model are violated, the model parameters can be infinite and estimates can be very sensitive. Later, a simulation will show that similar problems can occur in bivariate data.

Transformation

As in the univariate setting, in the bivariate setting it is commonly recommended that skewed incomplete variables be transformed to approximate normality before imputation (Schafer and Olsen 1998:550; Raghunathan et al. 2001:82-83; Allison 2002:39; Schafer and Graham 2002:167). But transformation of incomplete variables is often unnecessary and can lead to bias.

Again, it is helpful to return to a simple example where X is standard exponential and Y is conditionally normal (see equation (4)). What should we do if Y is MAR? Y has some skew which it inherits from X, so we might be tempted to transform Y to a better approximation of normality. Yet transformation could only hurt the imputations, since the optimal imputation model is a linear regression that imputes Y on its original scale. Linear regression imputation is appropriate for Y because, although the marginal distribution of Y is skewed, the X−Y relationship is linear and the conditional distribution of Y is normal—that is, the residual $e_{Y | X}$ is normal in the regression of Y on X (see equation (4)). Transformation can only hurt the imputations by introducing nonlinearity and residual nonnormality.

So transformation can only hurt the imputation of Y. Transformation can also hurt the imputation of X, because transformation introduces curvature into the relationship between X and Y. In Figure 4, for example, we have an exponentially distributed X, which, as we saw earlier, can be transformed to a very good approximation of normality by a fourth root transformation $t (X) = \sqrt[4]{X}$ (Hawkins and Wixley 1986). Yet, if we use this transformation in a bivariate setting, the imputation model will be misspecified and incompatible with the analysis model. The analysis model is a linear regression that correctly assumes a linear relationship between Y and X, yet under transformation the imputation model assumes a linear relationship between Y and $\sqrt[4]{X}$ .

As an alternative to transforming X alone, we can transform all the variables in the imputation model. For example, if X is exponential, we can transform X to $\sqrt[4]{X}$ and Y to $\sqrt[4]{Y - M i n (Y)}$ where $M i n (Y)$ is the smallest observed Y value, subtracted to ensure that the value under the radical is zero or more. Here the purpose of transforming Y is not to normalize Y, but to put the complete variable Y on a scale that is more nearly linear with respect to the transformed incomplete variable X. If both X and Y were incomplete, the choice of transformation would be less straightforward since different transformations would likely be optimal for each variable. And if there were more than two incomplete variables, choosing a suitable transformation could become unwieldy.

Simulation Experiment

In this section, we carry out a simulation experiment to test the bias and efficiency of different methods for imputing a skewed X variable in a bivariate (X, Y) data set. We then extend the simulation to accommodate a third variable Z.

Bivariate Design

In each simulated bivariate data set, Y fits a normal linear regression on X:

Y = α_{Y | X} + β_{Y | X} X + e_{Y | X}, w h e r e e_{Y | X} ~ N (0, σ_{Y | X}^{2}) .

We hold the slope and intercept constant at

α_{Y | X} = β_{Y | X} = 1

since changing these values to any nonzero value would be equivalent to a trivial shifting or rescaling of the axes. We also hold the sample size constant at n = 100 observations, so that we can avoid being distracted by small sample biases⁵ and focus on the asymptotic properties of the imputation methods.

The experiment independently manipulates four factors:

The first factor is the distribution of X. We let $X ~ χ_{ν}^{2}$ be a chi-square variable where the degree of freedom $ν$ determines the mean $E (X) = ν$ , the variance $V (X) = 2 ν$ , and the coefficient of skewness $S k e w (X) = \sqrt{8 / ν}$ . We set the degrees of freedom to four different levels $ν = 1, 2, 4, 8$ so that the skew ranges from mild (1) to severe ( $\sqrt{8} \approx 2.8)$ . Note that with $ν = 2$ the chi-square distribution reduces to the standard exponential distribution that we have used in all our examples up to this point.

The second factor is the strength of the X–Y relationship as measured by the regression’s coefficient of determination $R_{Y | X}^{2}$ . We set the coefficient of determination to four different levels $R_{Y | X}^{2} = .1, .3, .5, .7, .9$ by first choosing the desired value of $R_{Y | X}^{2}$ and then setting $σ_{Y | X}^{2} = β_{Y | X}^{2} V a r (X) (1 - R_{Y | X}^{2}) / R_{Y | X}^{2}$ .

The third factor is the location of the missing X values. From each complete data set, we delete half the X values in three MAR patterns:

To let values be missing completely at random (MCAR), we delete X values with a constant probability of ½.

To favor deletion of large X values—that is, values in the tail of the X distribution—we delete X values with probability $\hat{F} (Y)$ where $\hat{F} (Y)$ is the empirical cumulative distribution of Y.

To favor deletion of small X values—that is, values near the peak of the X distribution—we delete X values with probability $1 - \hat{F} (Y)$ .

Note that we do not vary the fraction of missing X values, since it is clear that changing the fraction of missing values would not change the relative performance of the imputation methods. Deleting a larger or smaller fraction of values would just make the differences between the imputation methods more or less consequential.

The fourth factor is the method of imputation. We use the seven different methods described earlier:

Linear regression imputation where X is regressed on Y.

Linear regression imputation with imputed X values censored below $c = 0$ .

Linear regression imputation with imputed X values truncated below $c = 0$ .

Quadratic regression imputation, where X is regressed on Y and $Y^{2}$ .

The transform X method: linear regression imputation where $\sqrt[4]{X}$ is regressed on Y. For the chi-square variable $X$ , $\sqrt[4]{X}$ is very close to normal for any degrees of freedom $ν$ (Hawkins and Wixley 1986).

The transform all method: linear regression imputation where $\sqrt[4]{X}$ is regressed on $\sqrt[4]{Y - M i n (Y)}$ . Here $M i n (Y)$ is the smallest Y value in each simulated data set, subtracted to ensure that the value under the radical is nonnegative.

Truncated regression imputation, with the lower truncation point at $c = 0$ .

For each imputation method, we imputed each incomplete data set five times. In each of the five imputed data sets, we calculated the mean, standard deviation, and skew of X, and we fit an OLS regression of Y on X. We then averaged estimates from the five imputed data sets to yield a single set of MI estimates for each incomplete data set and imputation method.

All the imputation methods used the MI procedure in SAS 9.2, except for the truncated regression method, which used the mi impute truncreg command in Stata 12. Both SAS and Stata displayed some technical limitations. These limitations are not central to the results, but are worth describing briefly, as follows:

In truncating normal imputations, SAS’s MI procedure uses a simple rejection algorithm,⁶ rejecting negative values and reimputing them iteratively until a positive value randomly occurs. The wait for a positive imputation can be long, and in about 8 percent of data sets SAS terminated the attempt after 100 iterations. We reimputed these data sets with a lower truncation point of $c = - 6$ , which means that the truncation was not as severe as we wanted it to be. We could have maintained a truncation point of $c = 0$ if SAS truncated imputed values by sampling importance resampling, which is used by IVEware (Raghunathan et al. 2001, 2002). Sampling importance resampling is faster than simple rejection and ensures that imputed values are within bounds (Robert 1995).

In transforming imputed variable, SAS’s MI procedure⁷ refuses to inverse-transform any values that could not have been obtained by transformation. In imputing X, for example, we transformed the observed X values to $X_{o b s, t} = \sqrt[4]{X_{o b s}}$ and then supplemented the observed values with normal imputations $X_{i m p, t}$ which included some negative values. When we tried to invert the transformation by calculating $X_{i m p} = X_{i m p, t}^{4}$ , SAS terminated on encountering negative values of $X_{i m p, t}$ because a fourth root cannot be negative. We got around this problem by transforming and inverse-transforming the variable outside of the MI procedure.

Stata’s implementation of truncated regression imputation is quite slow, taking about twenty-four hours to impute all the experimental data sets (SAS imputed them in minutes). It is not clear whether this is an inherent problem with truncated regression or a problem in Stata’s implementation of it. In addition to being slow, the truncated regression model failed to converge in about 1 percent of incomplete data sets—typically data sets with low $ν$ and low $R_{Y X}^{2}$ . While again this could be an implementation issue, we suspect that failure to converge resulted from the truncated regression model’s inherently poor fit to a chi-square variable X with low degrees of freedom $ν$ . Recall that, even in the univariate setting, we found that fitting a truncated normal model to a chi-square variable with $ν = 2$ could lead to infinite parameters and unstable parameter estimates. After failures to converge, we refit the truncated regression successfully by moving the lower truncation point to $c = - 1$ .

Illustrative Results

Figure 5 uses all the imputation methods to impute a simulated data set with $ν = 2$ , $R_{Y X}^{2} = .3$ , and X values MCAR. The true regression line—that is, the line that fits the complete-data population—is shown for reference.

Figure 5.

Different imputation methods when the incomplete variable X is standard exponential with values missing completely at random, and the complete variable Y is conditionally normal.

Clearly, some methods handle these data better than others. Under linear regression imputation, the imputed values, although more symmetrically distributed than the observed values, nevertheless fit well around the true regression line. Censoring or truncating the linear regression imputations does not change the fit very much; censoring and truncating simply move the leftmost values a little to the right, and this has little influence since the points with the most influence on the regression line are on the far right, in the tail of the X distribution. Quadratic regression imputation produces fairly similar results, although the imputed X values are a little more dispersed than they are under linear regression imputation.

The transform X method and the transform all method have serious trouble with these data, because transformation curves the relationship between X and Y. In the tail, which has the most influence on the regression, all the imputed points lie below the true regression line. These points will negatively bias the estimated slope.

The truncated regression has even worse problems, with nearly all the imputed points in influential positions and below the true regression line. Evidently, the truncated regression model can be a poor fit to observed data that are not actually truncated, and very unrealistic parameter estimates can result. Recall that the truncated model can produce infinite parameter values even in a simple univariate setting.

Figure 6 gives even more serious examples of bad imputations obtained from the truncated regression model. In one example, the imputed X values have a negative relationship to Y, even though the relationship in the observed data is clearly positive. Some imputed X values have values of 200 or greater, even though all the observed X values are less than 11. In one example, nearly all the imputed values are negative, even though the purpose of truncated regression here is to produce positive imputations. The presence of negative imputations, though extremely rare, suggests implementation problems that go beyond the basic difficulty of fitting a truncated model to nontruncated data.

Figure 6.

Anomalous imputed data sets generated by truncated regression imputation in Stata 12. Note: Main effect of imputation method, averaged across other factors.

Comprehensive Results

The full bivariate simulation experiment yielded 36,000 sets of MI estimates—100 estimates for each of 360 different experimental conditions (4 $ν$ values × 5 $R_{Y X}^{2}$ values × 3 missing value patterns × 6 imputation methods). Within each experimental condition, we calculated the mean and standard deviation of the 100 MI estimates in order to estimate the expectation and standard error of each imputed-based estimator under each experimental condition. From the expectation and standard error, we calculated the bias and root mean square error (RMSE). Finally, we divided the bias and RMSE by the true value of the estimand in order to get the relative bias and relative RMSE—that is, the bias and RMSE expressed as a percentage of the true value of the estimand.

Table 1 summarizes the relative bias and relative RMSE of each imputation method, averaged across all the experimental conditions. Results are shown for seven different estimands. The first three estimands describe the marginal distribution of X in terms of the mean, standard deviation, and coefficient of skewness. The remaining four estimands describe the regression of Y on X in terms of the intercept, slope, residual standard deviation, and coefficient of determination.

Table 1.

Results of Two-Variable Simulation: Main effect of imputation method.

	Marginal Distribution of Incomplete X						Regression of Complete Y on Incomplete X
	Mean		SD		Skew		Slope		Intercept		Residual SD		R ² _Y\|X
Imputation method	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)
Linear regression	−6	15	−3	22	−52	57	1	26	15	95	0	12	3	37
Linear regression, censored	−1	13	−8	20	−37	46	5	27	−1	93	1	12	0	35
Linear regression, truncated	2	16	−9	21	−38	46	4	28	−6	96	2	12	−5	35
Quadratic regression	0	16	6	25	−31	57	−11	31	37	113	1	13	3	41
Truncated regression	⁴⁰⁴⁹³⁶	^4064963	⁴¹⁷³⁹³	^4181515	−29	43	−8	28	18	102	0	11	5	39
Transform X	5	21	27	62	21	48	−17	30	47	106	14	21	−11	37
Transform all	−5	16	−1	30	−5	29	−20	31	74	118	25	32	−31	50

Note: RMSE= root mean square error.

For all but one of the estimands, linear regression imputation yields estimates whose average relative bias is fairly small, often close to zero. Censoring or truncating the linear regression imputations yields similar results, with a little more bias for some estimands and a little less for others.⁸ All the other methods have more serious biases. The transformation methods have the worst biases for the regression slope and intercept, and the truncated regression method has the worst biases for the mean and standard deviation of X. The biases of the truncated regression method can be enormous, exceeding 400,000 percent for some parameters, because the method occasionally imputes wild outliers like those in Figure 6.

In short, for the estimands that are emphasized in most social research—means, standard deviations, and regression parameters—linear regression imputation, with or without censoring or truncation, can often produce reasonable results. The other imputation methods are no better for these estimands, and sometimes much worse.

This is not to say that linear regression imputation is good for every estimand—it is not. In fact, linear regression imputation does a very poor job of estimating parameters that reflect distributional shape. In estimating the coefficient of skewness, for example, linear regression imputation had an average relative bias of –52 percent—that is, the estimated skew was less than half the true skew, on average. With one exception, the other methods do not estimate skew very well either; they are less biased than linear regression imputation but still have positive or negative biases exceeding 20 percent. The one method that estimates skew with little bias is the transform all method, but we cannot recommend that method because it has serious biases in estimating the regression parameters.

The online Appendix Table A1 (which can be found at http://smr.sagepub.com/supplemental/) summarizes the simulation results in more detail, breaking them down across different levels of the manipulated factors. The most striking result is that the differences of the methods are largest when missing values are in the tail. This makes sense because the tail values have the most influence in estimating the regression line. Another striking result is that the biases of the transformation methods are worst when X is highly skewed—which is exactly when one would be most tempted to use transformation.

The biases of the different methods vary from one simulated condition to another, and there are circumstances where transformation or quadratic regression or truncated regression gives the best results. It is tempting to imagine that you could obtain good results by picking and choosing different methods to suit different circumstances, but that it is a dangerous game with only small benefits for guessing right, and large penalties for guessing wrong.

The safest approach is to use a method that tends to have small biases in a wide variety of settings, and by that criterion the best choice for bivariate data is linear regression imputation, with or without truncation or censoring.

Trivariate Extension

We now extend the simulation to a regression where there are three variables: a complete dependent variable Z and two independent variables, one complete (Y), and one incomplete (X). To do this, we simply keep the X and Y variables from the two-variable simulation, and add a Z variable that fits a linear regression on X and Y, with normal residuals:

Z = α_{Z | X Y} + β_{Z | X} X + β_{Z | Y} Y + e_{Z | X Y}, w h e r e e_{Z | X Y} ~ N (0, σ_{Z | X Y}^{2}) .

The manipulated factors are the same as before, with only minor adjustments:

The first factor is the distribution of X. Again $X ~ χ_{ν}^{2}$ is a chi-square variable with degrees of freedom $ν = 1$ ,2,4,8.

The second factor is the squared correlation $ρ_{X Y}^{2}$ between X and Y. This is the factor that we manipulated in the bivariate simulation, except that in the bivariate simulation $ρ_{X Y}^{2}$ was interpreted as the coefficient of determination $R_{Y | X}^{2}$ for the regression of Y on X.

The third factor is the pattern of missing X values. These patterns are defined in exactly the same way as in the bivariate setting. In the MCAR scenario, X values are deleted with probability ½; in the MAR tail pattern, X values are deleted with probability $\hat{F} (Y)$ ; and in the MAR peak pattern, X values are deleted with probability $1 - \hat{F} (Y)$ . Note that Z has no effect on whether X is missing, and that both Y and Z are complete.

The fourth factor is the imputation method. The methods are the same as in the bivariate simulation, except that X is imputed conditionally on Z as well as Y. In the transform all method, X is transformed to $\sqrt[4]{X}$ , Y is transformed to $\sqrt[4]{Y - M i n (Y)}$ , and Z is transformed to $\sqrt[4]{Z - M i n (Z)}$ , where Min(Y) and Min(Z) are the minimum Y and Z values in the imputed data set, subtracted to ensure that the value under the radical is nonnegative.

To avoid further complicating the experiment, we hold constant certain parameters of the regression of Z on X and Y. The constant parameters are the slopes and intercept, which are held constant at $α_{Z | X Y} = β_{Z | X} = β_{Z | Y} = 1$ , and the coefficient of determination, which is held constant at $R_{Z | X Y}^{2} = .5$ by setting $σ_{Z | X Y}^{2} = V (β_{Z | X} X + β_{Z | Y} Y)$ .

Table 2 summarizes the relative bias and relative RMSE of each imputation method for estimands relating to the marginal distribution of X and the regression of Z on X and Y. Table 2 presents average results; in the online appendix (which can be found at http://smr.sagepub.com/supplemental/) Table A2 breaks the results down by each manipulated factor.

Table 2.

Results of Three-Variable Simulation: Main effect of imputation method.

	Marginal distribution of incomplete X						Regression of complete Z on incomplete X and complete Y
	Mean		SD		Skew		Slope of X		Slope of Y		Intercept		Residual SD		R²_Y\|X
Imputation method	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)	Rel. bias (%)	Rel. RMSE (%)
Linear regression	−5	14	−2	22	−51	57	−2	65	4	42	−1	157	−1	8	2	17
Linear regression, censored	−1	13	−7	19	−37	45	2	66	4	40	−18	155	−1	8	2	17
Linear regression, truncated	2	15	−8	20	−37	46	−1	65	6	41	−25	161	0	8	1	17
Quadratic regression	0	15	7	25	−29	57	−13	64	7	41	9	155	−1	8	2	17
Truncated regression	^78933932	^607613905	^71427514	^546826216	−27	42	−6	64	2	40	−4	157	−1	8	2	17
Transform X	6	20	30	63	24	51	−27	57	18	37	3	142	0	7	0	17
Transform all	−4	16	0	32	−5	30	−31	60	22	41	11	146	1	8	−1	17

Note: RMSE= root mean square error.

The basic results of the three-variable simulation are similar to those of the two-variable simulation.

In estimating the regression parameters, and in estimating the mean and standard deviation of X, linear regression imputation gives perhaps the best results overall, with relative biases of 5 percent or less. Quadratic regression imputation gives similar results, as does censoring or truncating the linear regression imputations. The transformation methods give more biased regression parameters, and the truncated regression method occasionally imputes very large outliers and so gives seriously biased estimates of the mean and standard deviation.

In estimating the skew, all the methods give highly biased results except for the transform all method, which cannot be recommended since it gives the most biased estimates for the regression parameters. Again, the differences among the imputation methods are much more consequential if values are imputed in the tail rather than the peak of the X distribution (see Table A 2b, which can be found in the online Appendix at http://smr.sagepub.com/supplemental/).

In the three-variable simulation, the slope of Z on Y has bias, which is opposite in direction to the bias of the slope of Z on X. Bias in the slope of Z on Y occurs despite the fact that neither Z nor Y has any imputed values. The reason for this bias is that Y is correlated with X, so any bias in the slope of X engenders a compensating bias in the slope of Y. Table A2b in the online appendix (which can be found at http://smr.sagepub.com/supplemental/) confirms that the bias in the slope of Y is greatest when the squared correlation $ρ_{X Y}^{2}$ between X and Y is large. Conceptually, the bias is similar to the bias that would occur if X were omitted or measured poorly.

An initially surprising result of the three-variable simulation is that the regression slopes estimated by the transformation methods, despite having the most bias, sometimes have smaller RMSE than the regression slopes estimated by other methods. The presumed reason for this is that the transformation methods, because they use misspecified imputation models, reduce the correlation between Y and the imputed X, and this reduced correlation yields smaller standard errors when X and Y are used to predict Z. Under some circumstances, the reduction in standard error can more than make up for the increase in bias.

Is the potential for reduction in RMSE a reason to use the transformation methods? No. First, reduced RMSE in the regression slope is limited to the simulations where the correlation between X and Y is very, very high. Table A 2b in the online appendix (which can be found at http://smr.sagepub.com/supplemental/) shows reduced RMSE in the simulation where $ρ_{X Y}^{2} = .9$ (i.e., $ρ_{X Y} \approx .95$ ) but no reduction in RMSE if $ρ_{X Y}^{2} = .7$ (i.e., $ρ_{X Y} \approx .84$ ). Second, because RMSE was reduced by accident, we have no theory for exactly when a similar reduction in RMSE might be expected outside the simulation. Analysts who use the transformation methods are risking severe bias, with only a speculative hope of reduced RMSE. Analysts who are willing to trade some bias for a reduction in RMSE would be better advised to use a more principled method with stronger theoretical foundations, such as ridge regression (Hoerl and Kennard 1970; Muniz and Kibria 2009).

Applied Examples

Authors often use incomplete skewed variables in applied research. In this section, we discuss two examples and evaluate the decisions that were made by the authors.

Analysis of Female Legislative Candidates

In a cross-national study, Kunovich and Paxton (2005) pointed out a strong bivariate relationship between two percentages that vary across the world’s n = 171 countries: the percentage of parliamentary legislators who are female (Y) and the percentage of parliamentary candidates who were female (X). The regression of Y on X is important because a slope of less than one suggests that female candidates lose more often than they win. Y was complete and skewed to the right, with a coefficient of skewness of 1.4. The skew of X is similar but harder to estimate from the observed data, since X was missing for more than half (ninety-nine) of the countries.

Following my advice, the authors imputed missing X values using the transform X method. It now appears, in light of results in this article, that my advice was misguided and could have caused substantial bias if missing X values had been primarily in the tail. Fortunately, X was missing primarily in the peak, and in this situation—with a skew of 1.4, values missing from the peak, and $R_{Y X}^{2} = .63$ —simulation suggests that the transformation method is nearly unbiased, while imputation without transformation has a bias of about –12 percent (see Online Appendix A found at http://smr.sagepub.com/supplemental/). Note that the authors used a milder transformation than the fourth-root transformation evaluated in our simulation. In particular, the authors used the arcsine transformation ${sin}^{- 1} \sqrt{X}$ , which is practically indistinguishable from $\sqrt{X}$ for the range of X values in their data.

Reanalyzing the authors’ data without control variables, I found that the estimated slope of Y on X was .60 with the authors’ transformation, and .55 without transformation. My simulations suggest that the authors’ estimate is closer to the truth, but the difference is not large either way.

Analysis of Body Mass Index (BMI)

In a longitudinal study of n = 358 children, von Hippel, Nahhas, and Czerwinski (2012) estimated percentiles for change in BMI from age 3½ to age 18 years. BMI was measured every six months, but some measurements were missing. Our longitudinal imputation model was more complicated than the models considered in this article, but the model still assumed that the incomplete variables were conditionally normal. Skew in BMI was an important challenge, especially since the skew of the BMI distribution increases as children grow older.

As one of the authors, I knew that linear regression imputation could estimate the conditional mean of a skewed variable with little bias, but this was not reassuring since in this study we sought to estimate the percentiles instead of the mean. With or without transformation, it proved difficult to find an imputation model that gave plausible estimates for extreme tail percentiles such as the 90th or 95th.

In the end, we sidestepped the issue by imputing BMI increments—that is, changes in BMI from one measurement occasion to the next. Although BMI is skewed, the distribution of BMI increments is approximately symmetric and has only slightly more kurtosis than a normal distribution.

Conclusion

On the whole, the simulation results suggest that when an incomplete variable has skew, linear regression often gives reasonable, though not unbiased, estimates for the quantities most commonly estimated in social science—namely, means, standard deviations, and linear regressions (see also Demirtas et al. 2008). Ad hoc modifications of linear regression—through censoring, truncation, or transformation—rarely do much to improve the estimates, and in fact can make the estimates much worse.

Although the normal regression method is fairly good for estimating means, variances, and regressions, it can do a poor job of estimating quantities that depend on distributional shape. Such quantities include the coefficient of skewness, the percentiles, and quantities based on percentiles such as the Gini coefficient.

Although our discussion has been confined to normal imputation, we note that similar biases would be expected if estimates were obtained, without imputation, by applying maximum likelihood to the incomplete data under an assumption of normality. Maximum likelihood can be asymptotically biased when its distributional assumptions are violated (Yuan, Wallentin, & Bentler, in press). For example, in the Bivariate Data section, we discussed a situation where OLS estimates for the regression of X on Y—which are maximum likelihood estimates if X is conditionally normal—are biased for a skewed X.

To improve estimation from incomplete skewed variables, it would be helpful to have imputation methods that do not assume normality. One option is the distribution-free approach of imputing missing values by resampling values from similar cases. Variants of this idea are called hot-deck imputation, the approximate Bayesian bootstrap, or predictive mean matching (for a review, see Andridge and Little 2010). For example, in bivariate data (X, Y) with Y complete and X MAR, one would fill in missing X values with observed X values resampled from cases with similar values of Y.

While imputation by resampling is initially attractive, it can work poorly when observed values are sparse in one part of the distribution. For example, suppose that X is missing if and only if Y < 0. If we wish to impute the missing X values in cases with Y < 0, we have to resample X values from cases with Y ≥ 0—and this inevitably leads to bias. (See Allison 2000 for a similar but slightly more complicated example.) To take a less-artificial example from applied research, in our study of BMI growth among children, the very heaviest children had a number of missing measurements. These measurements could not be imputed by resampling, unless we were willing to resample BMIs from much lighter children.

A perhaps more promising idea is to model and impute nonnormal variables using flexible nonnormal distributions that can take a variety of shapes—such as Tukey’s gh distribution (He and Raghunathan 2006), the Weibull or beta density (Demirtas and Hedeker 2008), the generalized lambda distribution, or Fleishman power polynomials (Demirtas 2010). These approaches are currently in the development stage. Initial evaluations suggest that they can mimic very well the shape of many nonnormal distributions and that they preserve the relationships among variables at least as well as normal imputation methods (He and Raghunathan 2006, 2009, 2012; Bondarenko and Raghunathan 2007; Demirtas 2010; Demirtas and Hedeker 2008). We look forward to these flexible methods being available in software so that they can be used in applied research and evaluated further.

Footnotes

Acknowledgments

I thank Ray Koopman for help with the calculations in the subsection on Imputing Within Bounds: Censoring and Truncation. I also thank Pamela Paxton and Sheri Kunovich for sharing the data in the subsection on Analysis of Female Legislative Candidates.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

References

Allison

Paul D.

2000. “Multiple Imputation for Missing Data: A Cautionary Tale.” Sociological Methods & Research 28:301–09.

Allison

Paul D.

2002. Missing Data. Thousand Oaks, CA: Sage.

Allison

Paul D.

2005. “Imputation of Categorical Variables With PROC MI.” Prsented at the SAS Users Group International, 30th Meeting (SUGI 30), Philadelphia , PA.

Anderson

T. W.

1957. “Maximum Likelihood Estimates for a Multivariate Normal Distribution When Some Observations are Missing.” Journal of the American Statistical Association 52:200–03.

Andridge

Rebecca R.

Little

Roderick J. A.

. 2010. “A Review of Hot Deck Imputation for Survey Non-response.” International Statistical Review 78:40–64.

Bernaards

Coen A.

Belin

Thomas R.

Schafer

Joseph L.

. 2007. “Robustness of a Multivariate Normal Approximation for Imputation of Incomplete Binary Data.” Statistics in Medicine 26:1368–82.

Bondarenko

Irina

Raghunathan

Trivellore E.

. 2007. “Multiple Imputations Using Sequential Semi and Nonparametric Regressions.” Proceedings of the Survey Research Methods Section. Alexandria, VA: American Statistical Association.

Box

G. E. P.

Cox

D. R.

. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society Series B 26:211–52.

Centers for Disease Control and Prevention. 2011. “A SAS Program for the CDC Growth Charts.” Retrieved July 13, 2011 (http://www.cdc.gov/nccdphp/dnpao/growthcharts/resources/sas.htm).

10.

Demirtas

Hakan

. 2010. “An Application of Multiple Imputation under the Two Generalized Parametric Families.” Journal of Data Science 8:443–55.

11.

Demirtas

Hakan

Freels

Sally A.

Yucel

Recai M.

. 2008. “Plausibility of Multivariate Normality Assumption When Multiply Imputing Non-Gaussian Continuous Outcomes: A simulation Assessment.” Journal of Statistical Computation & Simulation 78:69–84.

12.

Demirtas

Hakan

Hedeker

Donald

. 2008. “Imputing Continuous Data Under Some Non-Gaussian Distributions.” Statistica Neerlandica 62:193–205.

13.

Goldberger

Arthur S.

1981. “Linear Regression after Selection.” Journal of Econometrics 15:357–66.

14.

Greene

William H.

1999. Econometric Analysis. 4th ed. New York: Prentice Hall.

15.

Halmos

Paul R.

1946. “The Theory of Unbiased Estimation.” The Annals of Mathematical Statistics 17:34–43.

16.

Hawkins

Douglas M.

Wixley

R. A. J.

. 1986. “A Note on the Transformation of Chi-Squared Variables to Normality.” American Statistician 40:296–98.

17.

Yulei

Raghunathan

Trivellore E

. 2006. “Tukey’s gh Distribution for Multiple Imputation.” The American Statistician 60:251–56.

18.

Yulei

Raghunathan

Trivellore E.

. 2009. “On the Performance of Sequential Regression Multiple Imputation Methods with Non Normal Error Distributions.” Communications in Statistics - Simulation and Computation 38:856.

19.

Yulei

Raghunathan

Trivellore E.

. 2012. “Multiple Imputation Using Multivariate Gh Transformations.” Journal of Applied Statistics 39:2177–98.

20.

Heitjan

Daniel F.

Basu

Srabashi

. 1996. “Distinguishing ‘Missing at Random’ and ‘Missing Completely at Random’.” The American Statistician 50:207–13.

21.

Hoerl

Arthur E.

Kennard

Robert W.

. 1970. “Ridge Regression: Biased Estimation for Nonorthogonal Problems.” Technometrics 12:55–67.

22.

Horton

Nicholas J.

Lipsitz

Stuart R.

Parzen

Michael

. 2003. “A Potential for Bias When Rounding in Multiple Imputation.” The American Statistician 57:229–32.

23.

Johnson

Norman L.

Kotz

Samuel

Balakrishnan

. 1994. Continuous Univariate Distributions, Vol. 1. 2nd ed. New York: Wiley-Interscience.

24.

Jöreskog

Karl G.

. 2002. “Censored Variables and Censored Regression.” Retrieved November 10, 2011 (http://www.ssicentral.com/lisrel/techdocs/censor.pdf).

25.

Kenward

Michael G.

Carpenter

James

. 2007. “Multiple Imputation: Current Perspectives.” Statistical Methods in Medical Research 16:199–218.

26.

Kim

Jae Kwang

. 2004. “Finite Sample Properties of Multiple Imputation Estimators.” The Annals of Statistics 32:766–83.

27.

Kunovich

Paxton

Pamela

. 2005. “Pathways to Power: The Role of Political Parties in Women’s National Political Representation.” American Journal of Sociology 111:505–52.

28.

Little

Roderick J. A.

1992. “Regression With Missing X’s: A Review.” Journal of the American Statistical Association 87:1227–37.

29.

Little

Roderick J. A.

Rubin

Donald B.

. 1989. “The Analysis of Social Science Data with Missing Values.” Sociological Methods & Research 18:292–326.

30.

Muniz

Gisela

Golam Kibria

B. M.

. 2009. “On Some Ridge Regression Estimators: An Empirical Comparisons.” Communications in Statistics - Simulation and Computation 38:621–30.

31.

Raghunathan

Trivellore E.

Lepkowski

James M.

Hoewyk

John Van

Solenberger

Peter W.

. 2001. “A Multivariate Technique for Multiply Imputing Missing Values Using a Sequence of Regression Models.” Survey Methodology 27:85–95.

32.

Raghunathan

Trivellore E.

Solenberger

Peter W.

Hoewyk

John Van

. 2002. IVEware: Imputation and Variance Estimation Software. Ann Arbor, MI.

33.

Robert

Christian P.

1995. “Simulation of Truncated Normal Variables.” Statistics and Computing 5:121–25.

34.

Rubin

Donald B.

1976. “Inference and Missing Data.” Biometrika 63:581–92.

35.

Rubin

Donald B.

1987. Multiple Imputation for Nonresponse in Surveys. New York, NY: Wiley.

36.

Rubin

Donald B.

Schenker

Nathaniel

. 1986. “Multiple Imputation for Interval Estimation From Simple Random Samples With Ignorable Nonresponse.” Journal of the American Statistical Association 81:366–74.

37.

SAS Institute. 2001. “The MI Procedure.” SAS/STAT Software: Changes and Enhancements, Release 8.2. Cary, NC: SAS Institute. Retrieved May 23, 2011 (http://support.sas.com/rnd/app/da/new/802ce/stat/chap9/index.htm).

38.

Schafer

Joseph L.

Graham

John W.

. 2002. “Missing Data: Our View of the State of the Art.” Psychological Methods 7:147–77.

39.

Schafer

Joseph L.

Olsen

Maren K.

. 1998. “Multiple Imputation for Multivariate Missing-Data Problems: A Data Analyst’s Perspective.” Multivariate Behavioral Research 33:545.

40.

Stata Corp. 2011. “Multiple Imputation for Missing Data.” Retrieved July 2, 2011 (http://www.stata.com/stata11/mi.html).

41.

von

Hippel

Paul

. 2007. “Regression With Missing Ys: An Improved Strategy For Analyzing Multiply Imputed Data.” Sociological Methodology 37:83–117.

42.

von

Hippel

Paul

. 2009. “How To Impute Interactions, Squares, and Other Transformed Variables.” Sociological Methodology 39:265–91.

43.

von

Hippel

Paul

. in press. “The Bias and Efficiency of Incomplete-Data Estimators in Small Univariate Normal Samples.” Sociological Methods & Research.

44.

von

Hippel

Paul

Nahhas

Ramzi W.

Czerwinski

Stefan A.

. 2012. “Percentiles for Change in Body Mass Index (BMI) from Age 3½ to 18 Years.” Unpublished manuscript.

45.

Wooldridge

Jeffrey M.

2001. Econometric Analysis of Cross Section and Panel Data. 1st ed. Cambridge, MA: MIT Press.

46.

Yuan

K.–H.

Wallentin

Bentler

P. M.

(in press). ML versus MI for missing data with violation of distribution conditions. Sociological Methods & Research.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.09 MB