Modern Multiple Imputation with Functional Data

Electronic Health Record (EHR) or Clinical Data often require longitudinal statistical methods, which account for the correlation between repeated measurements on the same subject. If one also assumes that these are taken from a smooth curve or data generating process, then we can exploit tools from FDA, which can produce gains in terms of flexibility and statistical power (He et al.,, 2011; Carnahan-Craig et al.,, 2018; Goldsmith and Schwartz,, 2017). Since hospital visits can occur both infrequently and irregularly, we can’t directly apply many FDA techniques to them. They pose a challenge to the current methods as well as for imputation. To address these challenges and illustrate the effectiveness of our approach over the current methods, we wish to apply them to an EHR data set where we predict if a smoker will relapse or not based on their Blood Pressure (BP) recordings over a span of 18 months.

The data provided by the Penn State PaTH to Health, assists research that uses patient data from multiple sources to further scientific discoveries. The data set describes patient-level data variables in a standardized manner (i.e., with the same variable name, attributes, and other metadata) along with information on demographics, encounters, diagnoses, and procedures. More information can be found on their website.

While there is a wealth of research related to smoking and blood pressure (BP) (Wang et al.,, 2018; Hansson et al.,, 1996; Primatesta et al.,, 2001), our goal here is not to make deep scientific statements, but rather to illustrate the utility of our methods, which we hope will prove useful to practitioners. The highest risk of relapse for smokers is during the end of their first and second year after quitting (Herd et al.,, 2009; García-Rodríguez et al.,, 2013). We, therefore, focus on modeling the relapse of the patients based on monitoring their BP within the first two years, which may be useful to the practitioners designing interventions for patients at risk of relapse.

In general, EHR data sets vary significantly with the timing and regularity of the appointments, and clinical measurements are affected by errors of varying types and degrees (Daymont et al.,, 2017). Similar challenges apply to measurements in the PaTH data set where we have various kinds of clinical measurements and information recorded. The ability to characterize trajectories of sparse irregular data has potential applicability to many clinical questions. Though the term sparsity is somewhat subjective in the context of functional/longitudinal data, many of the patients in the Path data set have just 2 measurements, while the most number of measurements for any patient is 17 (after cleaning and implementing the exclusion criteria). We can see from Figure 1 the modal number of measurements is 2, while relatively few patients had more than 7 clinical visits and almost none had more than 11. Also, from the cumulative observation Figure (appendix Figure 6), we observe that 94% of the patients had 10 or fewer measurements, 72% had at most 5 measurements, and 28% had no more than 2 measurements. On an average, we have around 4 measurements per patient. Sparsity arises due to many reasons in an EHR setting. Some patients never come back, some patients are not observed with any uniformity or regularity, etc. Having identified the BP trajectories as both sparse and irregular, we move to introduce the imputation methods that account for these conditions in a Functional framework before revisiting this data in Section 3.

The rest of the paper is organized as follows. In Section 2, we briefly go through PACE and Multivariate Imputation methods (MICE and MF) before introducing our proposed method using LLF, and modifications to the Multivariate Imputation methods using bins and careful initialization, to better deal with Functional Data. We present multiple simulations for the linear and non-linear cases under different values of sparsity in Section 3. This section also includes the EHR data, where we fit a scalar-on-function regression model to determine if a patient (smoker) will relapse or not at the end of 18 months using BP as the functional predictor. These examples help us to illustrate the limitations of previous approaches and demonstrate the usefulness of our methodology, which overcomes many of the issues discussed earlier. In the final section, we present our concluding remarks and future research directions, which pertain to better understanding of the bins and deeper statistical theory.

We assume the data is collected from trajectories that are independent realizations of a smooth random function, with unknown mean function $E(X(t))$=$\mu(t)$ and covariance function $C_{X}(t,s):=cov(X(s),X(t))$. We define the underlying functional covariates as $\left\{X_{i}(t):t\in[0,1];1\leq i\leq n\right\},$ where $t$ denotes the argument of the functions, usually time, and $i$ denotes the subject or unit. We assume that these curves are only observed at times $t_{ij}$ ($j=1,\ldots,m_{i}$) with some error:

$$x_{ij}=X_{i}\left(t_{ij}\right)+\delta_{ij}$$

Let $\mathbf{x}_{i}=\left(x_{i1},\ldots,x_{im}\right)^{\top}$ denote the vector of observed values on the function $X_{i}.$ Let $Y_{i}$, be the outcome, which is a function of $X_{i}$ and some error. Examples of such relations can be found in the next subsection.

Generally, when integration is written without limits, it is implied to be over the entire domain, usually standardized to $[0,1]$ for simplicity. The main focus of this work is to develop tools for consistently estimating the parameters in functional regression models.

The scalar-on-function regression model in the linear case is defined as:

$$Y_{i}=\alpha+\int\beta(t)X_{i}(t)dt+\varepsilon_{i}$$

(1)

A common way of estimating the model components is by Basis or Functional Principal Component (FPC) expansions, where we simplify the problem of estimating the parameters by projecting the functions to a finite dimension and then using multiple regression and least squares (Kokoszka and Reimherr,, 2018; Ramsay and Silverman,, 1997).

Non-linear modeling becomes very challenging with functional data due to the curse of dimensionality. One popular way of simplifying the problem is by using the generalized additive model, also known as continuously additive model (McLean et al.,, 2014; Müller et al.,, 2013; Wang and Ruppert,, 2013; Fan et al.,, 2015; Ma and Zhu,, 2016):

$$Y_{i}=\int f(X_{i}(t),t)dt+\epsilon_{i}$$

(2)

where the bivariate function $(x,t)\rightarrow f(x,t)$ is smooth, but unknown. It is commonly estimated using either basis expansions (Yao et al.,, 2005; Greven et al.,, 2010) (often with a tensor product basis) or using Reproducing Kernel Hilbert Spaces (Wang and Ruppert,, 2013; Reimherr et al.,, 2017). McLean et al., (2014) developed Functional Generalized Additive Models that enabled non-scalar response mapping.

Most methods for fitting the models discussed above require densely sampled functional data. For irregular and sparsely sampled data that are observed with error, estimating the modeling parameters becomes much more challenging. Directly smoothing the $x_{ij}$ to plug into a dense estimation framework seems like a straightforward idea but can result in substantial bias.

PACE (Yao et al.,, 2005) uses functional principal components (FPC) analysis, in which the FPC scores are imputed using conditional expectations. In addition to the requirements discussed previously, PACE also relies heavily on the data being Gaussian. The Karhunen-Loève or Principal Component expansion of $X_{i}(t)$ is given by:

$$X_{i}(t)=\mu_{X}(t)+\sum_{j=1}^{\infty}\xi_{ij}v_{j}(t),$$

(3)

where $v_{j}(t)$ are the eigenfunctions of $C_{X}$ with eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$. The scores are computed as:

$$\xi_{ij}=\left\langle X_{i}-\mu_{X},v_{j}\right\rangle$$

(4)

PACE proceeds by computing the conditional expectation of the scores given the observed data. This conditioning method is straightforward and tends to work much better than direct smoothing of $x_{ij}$. It provides the Best Linear Unbiased Predictors (BLUPs) under Gaussian assumptions and works in the presence of both measurement errors and sparsity. We can plug in the BLUPs values into a dense estimation framework to model the response.

PACE still suffers from a few major problems. One issue is that the imputation procedure of PACE does not consider the response $Y_{i}$ nor does it have any consideration for subsequent models that will be fit. This results in a bias while estimating model parameters (Petrovich et al.,, 2018). In addition, PACE is just a single imputation method and hence the uncertainty in the imputation is not properly propagated when forming confidence intervals, prediction intervals, or p-values. For this reason, the PACE software (Chen et al.,, 2019) uses an alternative approach for fitting linear models, which does not extend to non-linear models.

After understanding PACE, we now look into some of the standard methods used for imputation in the multivariate case. All of these methods have proven to work well in the multivariate setting but have never been tested in the functional setting, where the sample sparsity can be very high.

Multivariate Imputation by Chained Equations (MICE) (van Buuren,, 2007), also known as “fully conditional specification” or “sequential regression multiple imputation”, has emerged in the statistical literature as one of the principal methods for addressing missing data. MICE performs multiple imputations rather than single imputation and hence, it can account for the statistical uncertainty. In addition, the chained equations approach is very flexible and can handle variables of varying types (e.g. continuous or categorical) as well as complexities such as bounds or survey skip patterns. At a high level, the MICE procedure is a series of models whereby each variable with missing data is modeled conditional upon the other variables in the data. The MICE procedure is as follows: 1. We start by initialization, wherein we fill in all the missing values with mean imputation. 2. Next, we select the first variable with missing entries. A model is fit with this variable as the outcome and the other variables as predictors. 3. The missing values of the current variable are then replaced with predicted (imputed) values from the model in step 2. 4. Step 2 and step 3 are repeated while rotating through the variables with missing values sequentially. The cycling through each of the variables constitutes one iteration or cycle. 5. At the end of one cycle, all of the missing values have been replaced with predictions from the model that reflect the relationships observed in the data. The cycles are repeated a few times and after each cycle the imputed values are updated.

The number of cycles to be performed is pre-specified and after the last cycle, the final imputations are retained, resulting in one imputed data set. Generally, 10-15 cycles are performed. The idea is that, by the end of the cycles, the distribution of the parameters governing the imputations (e.g. the coefficients in the regression models) should have converged, in the sense of becoming stable. Different MICE software packages vary somewhat in their exact implementation of this algorithm (e.g. in the order in which the variables are imputed), but the general strategy is the same. Here, we have used the MICE (Buuren and Groothuis-Oudshoorn,, 2011) package in R.

A key advantage of MICE is its flexibility in using different models. Generally, the modeling techniques included in MICE are Predictive Mean Matching, Linear Regression, Generalized Linear Models, Bayesian Methods, Random Forest, Linear Discriminant Analysis, and many more. Its primary disadvantage is that it does not have the same theoretical justification as other imputation methods. In particular, fitting a series of conditional distributions, as is done using the series of regression models, may not be consistent with proper joint distribution, though some research suggests that this may not be a large issue in applied settings (Schafer and Graham,, 2002).

Linear case:

For the Linear case, we simulate $n$ iid random curves $\left\{X_{1}(t),\cdots,X_{n}(t)\right\}$ from a Gaussian process with mean 0 and covariance

$$C_{X}(t,s)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right),$$

which is the Matérn covariance function and $K_{\nu}$ is the modified Bessel function of the second kind. We set $\rho=0.5,$ $\nu=5/2$ and $\sigma^{2}=1.$ These curves are evaluated at $m$ equally-spaced time points from $[0,1]$. We assume that each observed point contains a normal measurement error with mean zero and variance $\sigma_{\delta}^{2}=0.3.$ We set $\beta(t)=w\times\sin(2\pi t),$ where $w$ is a weight coefficient used to adjust the signal. The response, $Y_{i}$ $(i=1,\cdots,n)$ is computed using the model in Equation (1), where $\alpha=0$ and $\sigma_{\epsilon}^{2}=1.$ In the binary response case, we define $Y_{i}$ $(i=1,\cdots,n)$ using Bernoulli and logit link function in Equation (1). Finally, for each curve, we assume a percentage ($s$) of the $m$ time points are unobserved. After the scores are imputed, we fit a scalar-on-function regression model using these imputed curves.

For the linear case, we simulate the data sets of different sample sizes, $n\in\{200,500,1000\}$ (results for n=200 and n=1000 are included in the appendix); different number of observations per curve, $m\in\{32,52\}$; and with different sparsity level, $s\in\{Medium,High\}$. For sparsity levels, Medium means $50\%$ of the points are missing for each curve and High means more than $85\%$ of the points are missing for each curve. Also, the values of $m$ are taken such that it helps with the process of binning. Each of these settings is simulated 10 times. Since we are primarily interested in the accuracy of the final estimates $\hat{\beta}(t)$, $\hat{Y}(t)$ and $\hat{X}(t),$ we report the Root Mean Square Error (RMSE) or Prediction Error for each of them.

Table 1 and Table 2 indicates that, in general, irrespective of the number of points, all the binned methods perform better compared to other methods for Prediction Error or RMSE of Prediction, $\beta$ coefficients and Imputation when the sparsity is Medium and the sample size is 500. Also, we can see that the RMSE of imputation for PACE is the same for both the tables. This is because we are using the same sample curves to generate scalar and binary response. There is no clear winner between MF and MLLF within the binned methods with or without PACE initialization. Again, for High sparsity, we notice a similar behavior as before, irrespective of the number of points, all the binned methods perform better. As we increase the sparsity, it seems like MICE and MLLF perform worse. This happens mainly because they do a poor job of imputing the curves, the effects of which get compounded when estimating the parameters and modeling. Again, binned methods are the best with no clear winner. The results for the other cases with different sample sizes and scalar response yield similar performance to these tables and can be found in the appendix Table 6 for $n=200$ and appendix Table 7 for $n=1000$.

We can see from Figure 2, how each method does in estimating the $\beta$ coefficient. We can observe that most of the MissForest extensions are catching the right shape and doing a good job. The plot does not contain PACE, MICE and MLLF as their estimates were very poor, which is also reflected in the RMSE for the $\beta$ coefficients from Table 1.

Figure 3 shows an example of imputed curves under different methods for one random sample curve. Here, the main focus is to see how binning helps in not only doing better imputation as seen from Table 1 but also the fact that it gives much smoother results as compared to the MissForest methods without the binning. The same effect can be seen with Local Linear Forest methods as well. This plot is included in the appendix Figure 7.

All simulations parameters are the same as before, except the response $Y_{i}$ $(i=1,\cdots,n)$ is computed using the model in Equation (2), where $f(X_{i}(t),t)=5*sin(X(t)^{2}*t^{2})$. For the Non-linear case, we also simulate data sets of sample size $n$ as 500 with different number of observations per curve, $m\in\{32,52\}$; and with different sparsity, $s\in\{Medium,High\}.$ Each of these settings is simulated 10 times. Since we are primarily interested in the accuracy of the final output $\hat{Y}(t)$ and $\hat{X}(t),$ we report the Root Mean Square Error (RMSE) or Prediction Error for each of them. Another Non-Linear model result can be found in the appendix Table8.

From Table 3 and Table 4, we observe that our proposed approaches are outperforming PACE and MICE for Imputation irrespective of the number of points and sparsity. When it comes to prediction, our methods are still better than PACE and MICE but the gap is not as large compared to the linear case. Overall, we see the same trend as in the Linear case, with the binned methods outperforming every other method under various simulation settings.

The major takeaway from all the simulations is that our methods perform the best under various settings. This is because our methods impute smoother curves resulting in better modeling and smoother estimates of the beta coefficients irrespective of the relation between the response and the functional predictors.

New statistical tools are vital for data such as PATH, which are very large and have a great deal of underlying structure. We see the performance of the developed tools for imputation with this longitudinal/functional data. The electronic medical records contain information about smokers (patients) who irregularly come for a check-up at the hospital. We have the blood pressure (BP) readings along with some other measurement values at each check-up of the patients. From previous studies, we know that majority of the relapse among smokers occur within the first two years. For cleaning the data, the exclusion criteria were based on the number of longitudinal measurements (time points). Patients who had a smoking history (at least smoked for a year) with less than 2 measurements were excluded. After cleaning, we are left with 122 patients, of which 61 smokers relapsed and 61 smokers didn’t, where each smoker is under observation for 18 months. Hence, here the sample size (n) is 122 and the number of time points (m) is 18. The data is sparse naturally as the patients don’t come in for check-up regularly and sometimes the measurements are missing even though there was a visit due to unknown reasons. We build a model to predict if a patient will relapse or not using the Blood Pressure (BP) measurement over an 18 month period.

We can infer from Table 5 that the binned methods again are outperforming the other methods in prediction irrespective of linear or non-linear modeling. Also, since both the model results are so close, the true relation looks linear. This is further supported by Figure 4, where we see only the binned methods have smooth imputed curves with the black points denoting the observed values for that patient. Although all the estimated $\beta$ coefficients seem to follow the same trend in Figure 5, the methods with binning have much smoother results, leading to better interpretability. Smoothness is inherent to functional data and that is why binned methods are able to perform so well.

Also, Figure 5 suggests that patients with low BP or sudden changes in their BP have a higher risk of relapse. Also, the curve isn’t constant suggesting that acceleration/velocity of the BP curve is important. We did check and found out that the average BP was higher in the control group (no relapse) than for the cases (relapse). We feel there might be confounding variables and further analysis is needed, which is out of the scope of the project as we are only interested in demonstrating the efficacy of our methods for imputation and training the model, which results in better analysis and interpretation.