MissScore: High-Order Score Estimation in the Presence of Missing Data

We thank the reviewer for their valuable insights and the time dedicated to reviewing our paper. We address each of the concerns and suggestions in detail below.

Changes Addressing Weaknesses

We have incorporated your feedback on the related work concerning missing data and the differences between MissScore and MissDiff into Appendix Sections A and B. Additionally, we have outlined a few limitations of MissScore in the Conclusion section. All changes have been highlighted in blue for your reference. Thank you once again for your valuable suggestions.

Various synthetic data generation methods used to generate synthetic records in the TSTR.

The baseline methods described in the Baselines paragraph of section 4 handle missing data by either imputing (filling in) or deleting it before training. These models are trained on a complete dataset (after imputation or deletion) and then generate synthetic data in the TSTR. In contrast, our approach trains directly on data with missing values and generates new synthetic data without needing to fill in or remove any missing entries beforehand.

The observed instability of MissDAG and MissDiffAN in higher dimensions (d = 50), as seen in Figure 12.

Both MissDAG and MissDiffAN are expected to face challenges in high-dimensional settings (e.g. where d=50). Generally, MissDAG struggles with high-dimensional data, while MissDiffAN performs poorly with higher missing rates, as shown in Figures 7 through 11. MissDAG's limitations are due to the complex, computationally intensive nature of its EM approach. This iterative EM process, combined with rejection sampling, leads to scalability issues as the number of nodes increases, resulting in degraded performance and greater computational demands. MissDiffAN, on the other hand, requires substantial data for effective diffusion model training. However, biased imputation of missing values negatively impacts training, leading to unreliable performance. Surprisingly, MissDAG performs well in this scenario, consistent with results from MissOTM (Vo et al., 2024), where MissDAG also showed strong performance at a missing ratio with 50 nodes (see Figure 3). This may indicate that the specific characteristics of this dataset are favorable for MissDAG's approach.

MissForest’s occasional competitiveness in certain conditions, in Figures 6-8.

MissForest’s occasional competitiveness in certain conditions, as seen in Figures 6-8, can be attributed to its unique strengths in handling complex data structures. In Figure 6, the task involves complete data, where DAGMA (Bello et al., 2022)—noted for its effectiveness in causal discovery within ANM—is applied, performs well. However, in Figures 7 and 8, where data imputation is required, MissForest emerges as a top performer due to its proven ability to deliver high-quality imputations. Studies such as MissOTM and MissDAG (Gao et al., 2022; Vo et al., 2024) consistently rank MissForest as one of the best imputers for causal discovery task. Its strengths lie in managing datasets with intricate interactions and non-linear relationships between variables of varying scales and types. Furthermore, MissForest performs effectively in high-dimensional datasets, even when the number of variables substantially exceeds the number of observations, consistently providing robust imputation results (Stekhoven & Buhlmann, 2012).

References:

Vy Vo, He Zhao, Trung Le, Edwin V Bonilla, and Dinh Phung. Optimal transport for structure learning under missing data. arXiv preprint arXiv:2402.15255, 2024.

Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization. Advances in Neural Information Processing Systems, 35:8226–8239, 2022.

Erdun Gao, Ignavier Ng, Mingming Gong, Li Shen, Wei Huang, Tongliang Liu, Kun Zhang, and Howard Bondell. Missdag: Causal discovery in the presence of missing data with continuous additive noise models. Advances in Neural Information Processing Systems, 35:5024–5038, 2022.

Daniel J Stekhoven and Peter Buhlmann. Missforest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1):112–118, 2012.

We appreciate the reviewer’s thoughtful comments and the time spent on reviewing our paper. Below, we address each of the concerns and suggestions provided.

Question 1: Direct comparison between the proposed method and the approach by Meng et al. (2021), using only complete-case data and a conservative view on the paper’s contribution.

In our work, we introduce a new approach to handling missing data with score matching. Our method addresses not only the MCAR case but also MAR case, ensuring that our estimator remains consistent with the oracle estimator. In MCAR, while using only complete-case data without imputation yields a consistent estimator comparable to the oracle, it performs worse when training a model (Ouyang et al., 2023; Kim et al., 2022; Gao et al., 2022; Vo et al., 2024), also as shown in the sampling experiments in Section 4. Additionally, at high missing rates, entire rows may be deleted, necessitating better methods for handling missing data in the MCAR scenario. Furthermore, our method supports k-th order estimation in the presence of missing data, which can be valuable for downstream tasks that require high-order score estimation.

To provide meaningful insights, we compare our approach to the method by Meng et al. (2021), which is based on complete-case data only. In Table 1, we use the complete data (i.e., $\alpha=0.0$) to calculate the first- and second-order scores for comparison. As the missing ratio increases, the performance of our first- and second-order scores remains comparable to that of the complete data case, demonstrating robustness. We have also included the MAR results in Tables 3 and 4 of the updated appendix for further reference. Additionally, our focus was not limited to the MCAR mechansim; instead, we provide substantial evidence of our method's effectiveness across both MCAR and MAR scenarios through experiments in sampling (Section 4) and causal discovery (Section 5).

Question 2: Theorem 6 when $\mathbb{P}[m=0|x=x]\rightarrow0$.

Thank you for pointing this out. Considering the case in Theorem 6 where $\mathbb{P}[\mathbf{m}=0|\mathbf{x}=\mathbf{x}] \to 0$, this has significant implications: the inverse probability weights, which are proportional to $1/\mathbb{P}[\mathbf{m}=0|\mathbf{x}=\mathbf{x}]$, become extremely large, leading to instability in the objective function. However, under the missing at random (MAR) assumption, the missing probability is generally treated as a fixed missing ratio in most studies, rather than approaching zero (Seaman & White, 2013; Vansteelandt et al., 2010; Zhang et al., 2023).

Question 3: Practical Applications with higher-order scores.

For instance, third- and fourth-order scores can capture more nuanced properties of the denoising distribution $p(\mathbf{x}|\tilde{\mathbf{x}})$, such as estimating its skewness and kurtosis, thus providing deeper insights into its shape and behavior.

Question 4: Tunable coefficient $w$.

The weight $\omega$ balances the contributions of the two loss components, $L_{\text{DSM}}(\theta)$ and $L_{\text{D2SM}}(\theta)$, in the overall objective function $L_{\text{joint}}(\theta)$. Selecting an appropriate $\omega$ is task-dependent and involves consideration of the relative importance of each score in the given context. For example, if the primary goal is to capture the broader structure of the data, as represented by the first-order scores, a smaller $\omega$ is preferred to place greater emphasis on the first-order term. However, if the task requires capturing finer-grained detail through higher-order relationships in the data, as is the case with second-order scores, a larger $\omega$ is advantageous and enhancing sensitivity to nuanced patterns. For consistency, $\omega$ is set to 10 across all experiments, as suggested by Meng et al. (2021).

Question 5: Impact of possible model misspecification with logistic regression model under MAR.

To illustrate the impact of potential model misspecification by the logistic regression model, we conduct an experiment comparing the ground truth $p$ with the estimated $\hat{p}$ from the logistic regression model on the synthetic dataset described in Section 3.2, with varying missing ratios. The results are provided in the additional discussions in Section B of the appendix, in Tables 3 and 4, in the updated manuscript. As you suggested, potential model misspecification in the logistic regression model does impact the estimation of missing probability. However, this effect remains within a reasonable range in the variance-reduced version.

References:

Stijn Vansteelandt, James Carpenter, and Michael G Kenward. Analysis of incomplete data using inverse probability weighting and doubly robust estimators. Methodology, 2010.

Shaun R Seaman and Ian R White. Review of inverse probability weighting for dealing with missing data. Statistical methods in medical research, 22(3):278–295, 2013.

Imke Mayer, Erik Sverdrup, Tobias Gauss, Jean-Denis Moyer, Stefan Wager, and Julie Josse. Doubly robust treatment effect estimation with missing attributes. The Annals of Applied Statistics, 14(3):1409–1431, 2020.

Yuqian Zhang, Abhishek Chakrabortty, and Jelena Bradic. The decaying missing-at-random framework: Doubly robust causal inference with partially labeled data. arXiv preprint arXiv:2305.12789, 2023.

Chenlin Meng, Yang Song, Wenzhe Li, and Stefano Ermon. Estimating high order gradients of the data distribution by denoising. Advances in Neural Information Processing Systems, 34:25359–25369, 2021.

We thank the reviewer for the comments, and we appreciate the time you spent on the paper. Below we address the concerns and comments that you have provided.

Provide additional details to ensure sufficient clarity

Thank you for the suggestion! We will add more details to clarify how these new theorems relate to the previous score-matching theorems, both for the first-order and higher-order scores, in the presence of missing data, making it more accessible to readers with less technical background.

Question 1: All the theorems (Theorem 1 - Theorem 7) are stated in terms of the joint distribution of $(x, \tilde{x}, m)$. However, we only have the distribution of $(x_{\text{obs}}, \tilde{x}_{\text{obs}}, m)$. Are all the theorems valid under this observed distribution?

No, the theorems do not hold if we rely solely on the observed data distribution. This limitation is primarily due to the score network, which is a neural network. However, element-wise multiplication with the mask naturally mitigates the impact of replacing missing values with zeros when computing the objective.

Question 2: Key difference between your method and the Missdiff of Ouyang et al.

MissScore differs from MissDiff by focusing on learning high-order score with missing data, specifically leveraging the Hessian to enhance downstream performance. In contrast, MissDiff emphasizes unbiased data synthesis with a first-order diffusion-based approach. Additionally, MissScore introduces an oracle estimator under the MAR mechanism and introduce a novel causal discovery method that operates effectively in the presence of missing data.

Question 3: It is not immediate from the experiment section how your method performs compared to MissDiff. Can you present some real data analysis (say using Census data) to compare your method with MissDiff?

MissDiff also conducted experiments on Census data, but since their code is not publicly available, we can only make indirect comparisons. Given that the diffusion model is, in some ways, equivalent to the DSM, MissDiff’s performance would be similar when applied to the Census data. However, our method incorporates second-order score information, which can provide better results. Therefore, our analysis primarily focuses on comparing first- and second-order models, where the first-order model aligns with the approach used in MissDiff. The results are presented in Section E.5 of the appendix in the updated manuscript.

Our analysis consists of two parts: numerical and categorical data in the census dataset. For the numerical data, which includes variables such as age, final weight (fnlwgt), years of education completed (education-num), and hours worked per week, we examine the distribution of both the original data and the sampled data. This comparison is done using MissScore with first-order information via Langevin dynamics (equivalent to MissDiff) and MissScore with second-order information through Ozaki sampling. We observe that MissScore with second-order information captures finer details in the distribution. For instance, in the age range of 20-40, the fnlwgt range of 0-0.25, education-num range of 10-11, and hours per week range of 40-45, the second-order MissScore better approximates the shape of the original distribution.

For the categorical data, we analyze the education variable and again find that the second-order MissScore aligns more closely with the original distribution. These results indicate that MissScore with second-order information captures finer-grained details in the sampled dataset, compared to MissDiff (MissScore with first-order information).

The paper claims that their results show that some baselines are introducing bias, since those baselines achieve worse results than MissScore. However, I don't think this is strong evidence of bias, since there could be other reasons for worse performance.

Yes, however, pinpointing a single reason for the poorer performance of these imputation methods is challenging. One contributing factor might be that some imputation methods don’t fully account for the uncertainty inherent in missing data, which can impact performance. Regarding bias, as Ouyang et al. (2023) shows, an “impute-then-generate” approach may introduce bias by potentially reducing generation quality, which can affect downstream tasks. Additionally, as Kyono et al. (2021) note, conditioning on missingness indicators during imputation may create spurious dependencies. In contrast, MissScore learns distributions directly from the data without relying on imputation, may help alleviate these issues. However, for some tasks, methods like MissDAG and MissOTM—which use directed acyclic graph (DAG) structures as inductive biases—could achieve better results, especially in causal discovery tasks where structured dependencies offer an advantage.

References:

Chenlin Meng, Yang Song, Wenzhe Li, and Stefano Ermon. Estimating high order gradients of the data distribution by denoising. Advances in Neural Information Processing Systems, 34:25359–25369, 2021.

Yidong Ouyang, Liyan Xie, Chongxuan Li, and Guang Cheng. Missdiff: Training diffusion models on tabular data with missing values. arXiv preprint arXiv:2307.00467, 2023.

Stef Van Buuren and Karin Groothuis-Oudshoorn. mice: Multivariate imputation by chained equations in r. Journal of statistical software, 45:1–67, 2011.

Lovedeep Gondara and Ke Wang. Mida: Multiple imputation using denoising autoencoders. In Advances in Knowledge Discovery and Data Mining: 22nd Pacific-Asia Conference, PAKDD 2018, Melbourne, VIC, Australia, June 3-6, 2018, Proceedings, Part III 22, pp. 260–272. Springer, 2018.

Alexander Gain and Ilya Shpitser. Structure learning under missing data. In International conference on probabilistic graphical models, pp. 121–132. Pmlr, 2018.

Zhenhua Wang, Olanrewaju Akande, Jason Poulos, and Fan Li. Are deep learning models superior for missing data imputation in large surveys? evidence from an empirical comparison. arXiv preprint arXiv:2103.09316, 2021.

Trent Kyono, Yao Zhang, Alexis Bellot, and Mihaela van der Schaar. Miracle: Causally-aware imputation via learning missing data mechanisms. Advances in Neural Information Processing Systems, 34:23806–23817, 2021.

Thank you for your valuable comments and suggestions. We appreciate the time and effort you dedicated to reviewing our paper. Below, we address each of the concerns and feedback you provided.

The statements of Theorems 3-7 exclude the case where the probability of a missing value is 0, so all values are observed. What happens to the theory in this case? Do the theorems fail?

Thank you for bringing up this important point. In Theorems 3-7, when all values are observed, the theorems remain valid. This corresponds to the original higher-order case, as discussed in Meng et al. (2021). We have incorporated the change in the revised version of our manuscript.

Theorem for k-th order scores with under MAR.

Under the MAR assumption, the result will align with Theorem 5, except that the term $(\mathbf{1}-\mathbf{m})$ is incorporating the inverse probability weights. The proof will also follow the same structure as the proof of the original Theorem 5, with the adjustment for the inverse probability weights. I have revised and updated this result in Theorem 7, along with the corresponding proof. Thank you for bringing this to our attention.

Fidelity

The fidelity metric measures how well synthetic data replicates the statistical properties of real data, specifically its distributions and inter-column relationships. In this work, we use SDMetrics to assess fidelity in two main areas: (1) Column Shapes: Assesses how closely the distribution in each synthetic column matches that of the real data. Using the Kolmogorov-Smirnov (KS) statistic, which compares CDFs, SDMetrics calculates the KSComplement (1 - KS statistic), where a higher score indicates a closer match; (2) Column Pair Trends: Evaluates if the relationships (e.g., correlation) between column pairs in synthetic data align with those in the real data. The overall fidelity score averages these two components into a single metric of similarity, which we use to compare data quality, as in Ouyang et al. (2023).

How do you aggregate the utility metric values from the different predictive models?

The utility metrics—Accuracy, AUROC, and weighted-F1—are reported separately for each predictive model without aggregation. In Table 2 of the main text, utility is represented by Accuracy, while additional evaluations, including AUROC and weighted-F1, are provided in Appendix Tables 6 and 7.

Is the number of iterations in Figure 2 the number of iterations for training the score models, or the number of iterations of Langevin dynamics?

Apologies for the confusion—the number of iterations in Figure 2 refers to the number of iterations of Langevin dynamics. We have made updates in the revised version of our manuscript.

Why do the utility metrics for MissScore in Table 5 only decrease significantly with increasing missingness under MCAR, but not MAR or MNAR?

Thank you for bringing up this important point. The observed results can be explained by the nature of the missingness mechanisms. Under MCAR, data is missing independently of both observed and unobserved values. As the missing rate increases, information is lost randomly across samples and features, leading to a consistent degradation in model performance. In contrast, MAR and MNAR involve dependencies that allow the model to better handle missing data. Under MAR, missingness depends on observed data, enabling the model to leverage correlations between observed variables and missingness patterns. With MNAR, missingness is related to unobserved data (e.g., missing values correlate with the missing data points themselves). Although MAR and MNAR generally present challenges for imputation, some methods can incorporate assumptions about these mechanisms to manage missing data more effectively (Van buuren & Groothuis-Oudshoorn, 2011; Gondara & Wang, 2018; Gain & Shpitser, 2018; Wang et al., 2021; Ouyang et al., 2023). Depending on the strength of these dependencies, the model may be able to sample new data more accurately from the learned model.

Thank you for the reply. You addressed my minor concerns and comments, but I did not see anything about my comment on the lack of realistic experiments. Since that is the main weakness of the paper in my opinion, I'm keeping my score for now.