Geodesic Optimization for Predictive Shift Adaptation on EEG data

Dear Reviewer r2hf,

We thank you for your detailed review. We took into account your comments and modified the experiments accordingly:

Addressing Weaknesses:

1. Knowledge of the target label mean:
   We acknowledge your concern regarding potential information leakage due to the assumption of knowing the target label mean value. We now estimate target $\bar{y}$ on target splits (50% of the data) that do not overlap with the evaluation target splits (50% of the data), rather than assuming target $\bar{y}$ to be known, ensuring no leakage from the test set. The revised Figure 2, included in the rebuttal PDF, shows that the performance remains comparable to when target $\bar{y}$ is assumed to be known.
   We emphasize that target $\bar{y}$ is the only quantity computed from the target y. In many scenarios, it is acceptable to assume that target $\bar{y}$ is known, such as when a hospital has prior knowledge about its patient population.
2. Comparison with the state-of-the-art:
   The Riemannian method “No DA” we use is already a strong baseline from recent works on regression from SPD matrices [e, f]. In response to your review, we have added “GREEN” [b], a recently proposed deep-learning method for biomarker prediction, to our experiments. We utilized the “g2” variant, which is an SPD network [d], since the data are SPD matrices. The revised Figure 2 shows that “GREEN” performs better than “No DA,” but still falls short of “GOPSA”. Although “GOPSA” shows superior performance, it is worth noting that the learned parallel transport could be integrated into “GREEN,” which we plan to explore in future work.
   We also included the “Re-scale” baseline [a, f] in our experiments on real data. This method corrects second-order statistics on the SPD manifold. However, as shown in the revised Figure 2, this additional baseline suffers from similar issues as “Re-center” and does not effectively resolve the problem of joint shifts in (X, y). Additionally, [f] studied rotation corrections but found that a meaningful method could only be derived where data must be paired between domains (i.e., the same patients are present in several domains), which is not feasible in our experimental setup.
   Overall, the revised Figure 2 clearly indicates that GOPSA remains the best-performing method compared to all tested methods.

Addressing Questions:

1. The performance difference between “GOPSA” and “DO Intercept”:
   From the revised Figure 2, we highlight that “GOPSA” significantly outperforms the “DO Intercept” which is the best performing baseline. Indeed, “GOPSA” significantly improved Spearman’s rank in 4 out of 5 site combinations (t-test p-values below 0.001). The results for the R² score and MAE are more nuanced, but on the majority of splits, both scores tend to be improved with “GOPSA” compared to “DO Intercept.”
2. Statistical assumptions of the t-test:
   We thank you for your remark regarding the assumptions of the t-test. Firstly, we recall that we revised the computation of the target $\bar{y}$, leading to some changes in the p-values from the initial submission. To address your concerns, we examined the normality assumption of the score differences presented in the revised Figure 2B. Although space constraints prevented us from including Q-Q plots, our analysis confirms that the score differences generally follow normal distributions. Exceptions were noted for the MAE and R² scores in the site combination of the first row (Ba, Cho, G, S) and all scores in the site combination of the fourth row (Cu03, M, R, S). We also computed variances of scores of “GOPSA” and “DO Intercept” and obtained that they are equal except in the aforementioned cases. In these cases, the reported p-values are above 0.05, indicating that the score differences might not be statistically significant from zero and thus we do not mislead the reader. However, with more careful testing (without normal assumptions), some of the p-values could potentially be lower. We appreciate your attention to this detail, as it helped us ensure the robustness and reliability of our findings.

Addressing Limitation:

1. Dataset concerns:
   We use the multi-centric dataset HarMNqEEG (2022) [h], which effectively serves as a group of datasets by combining measurements from various hospitals, with data from over 1,500 subjects. This dataset encompasses varying conditions of measurements, including data from 9 countries, 14 studies, and 12 different EEG devices. This combination of heterogeneous datasets explains the presence of a (X, y) shift between sites. Consequently, testing across different combinations of sites allows us to simulate experiments across various measurement conditions effectively. This approach enables us to treat these site combinations as distinct datasets, thus providing a robust assessment of the generalizability of our method.
   For future work, we welcome any recommendations for additional datasets that could further validate the generalizability and robustness of our method.

[a] Rodrigues, P. L. C., et al. (2018). Riemannian Procrustes analysis: transfer learning for brain-computer interfaces. IEEE Transactions on Biomedical Engineering

[b] Paillard, J., et al. (2024). GREEN: a lightweight architecture using learnable wavelets and Riemannian geometry for biomarker exploration. bioRxiv

[d] Huang, Z., & Van Gool, L. (2017, February). A Riemannian network for spd matrix learning. In Proceedings of the AAAI conference on artificial intelligence (Vol. 31, No. 1).

[e] Sabbagh, D., et al. (2019). Manifold-regression to predict from MEG/EEG brain signals without source modeling. Advances in Neural Information Processing Systems, 32.

[f] Mellot, A., et al. (2023). Harmonizing and aligning M/EEG datasets with covariance-based techniques to enhance predictive regression modeling. Imaging Neuroscience, 1, 1-23.

[h] Li, M., et al (2022). Harmonized-multinational qEEG norms (HarMNqEEG). NeuroImage, 119190.

Dear Reviewer ucmS,

We thank you for your detailed review of our paper. We considered your feedback and addressed the concerns and suggestions you raised to improve the clarity and impact of our research.

Clarification of Motivation and Problem Setup:

1. Motivation for domain adaptation on the SPD manifold:
   The necessity for domain adaptation (DA) methods is illustrated by the “No DA” approach, as detailed in the paper. This baseline corresponds to a SOTA covariance-based regression pipeline [e, f] without any adaptation between domains. As shown in Figure 2 from the PDF rebuttal, the “No DA” approach performs worse than a dummy model returning the average y per domain in terms of R² and MAE, demonstrating the critical need for adaptation techniques.
   Classical domain adaptation methods on the SPD manifold are the “Re-center” and “Re-scale” methods [a, f] which correct shifts in X between domains. However, they exhibit very poor performance when both distributions in X and y vary between domains, as shown in Figure 2.
   This joint shift issue on EEG data recently emerged with the ‘HarMNqEEG’ dataset [h], introduced in 2022, which we use in our experiments. Thus, this submission aims at presenting this problem and providing a new method called “GOPSA” to solve it.
2. Riemannian mixed-effects models:
   In the existing literature (see e.g. [g]), Riemannian mixed-effects models typically assume that X belongs to a Euclidean space while y is in a Riemannian manifold. Our paper considers the opposite scenario, with X on a Riemannian manifold and y as a real-valued response. To the best of our knowledge, we are the first to address this problem, highlighting the novelty and significance of our approach. Thus, the proposed method is not made “on top of a Riemannian mixed-effects model” but the first to address mixed effects with X on a Riemannian manifold and y real-valued.
3. Advantages of GOPSA:
   We acknowledge the reviewer’s concern regarding the need for a clearer explanation of the main advantages of our method, “GOPSA”. To address this, we have added numerical experiments on simulated data that demonstrate how “GOPSA” effectively handles joint shifts in both X and y. The results are presented in Figure 1 of the rebuttal PDF. First, we leverage the classical instantaneous mixing model: $x_i(t) = A\eta_i(t)$ where $x_i$ is the observed time-series, $i$ is the underlying signal of the neural generators and $A$ is the mixing matrix whose columns are the observed patterns of the neural generators. Furthermore, we use the log-linear model proposed in [e, f] $y_i = \beta_0 + \sum_{\ell=1}^d \beta_{\ell} \log(p_{\ell i})$ where $p_{\ell i}$ is the variance of the $\ell$-th element of the underlying signal $\eta_i$. From this, we generate domains by applying two shifts. One on X that changes the mixing matrix per domain, $x_i \mapsto B_k^{\xi} x_i$ with $B_k$ a domain-specific SPD matrix, $k$ the domain number and $\xi$ the variable that controls the intensity of the shift. A second shift is applied on y by shifting the variances per domain, $p_{\ell i} \mapsto p_{\ell i}^{(1+k\xi)}.$ It should be noted that $\beta$ is kept the same across the domains.
   The results are presented in Figure 1 of the rebuttal PDF. First, it shows that, if there is no shift in X, then "No DA" perfectly estimates the y because the log-linear model is respected across domains even when the y distribution changes (Figure 1B). However, the clear advantage of “GOPSA” is to estimate this log-linear model with shifts in (X, y) per domain (Figure 1C) which other methods can not do.
   These experimental results on simulated results are consistent with the results on real data: “GOPSA” is the best method to estimate a shift on X while domains have different distributions in y.
4. Domain shifts in biomedical data:
   Domain shifts in X often occur in biomedical data due to variations in data collection methods, such as differences in equipment, electrode placement, and environmental conditions. These factors can lead to inconsistent signal characteristics across different study sites. Shifts in y, the outcome variable, can result from demographic and clinical variations in the populations being studied, leading to differences in response distributions. While these issues are prominent in EEG data, they are common across various biomedical fields, necessitating robust adaptation techniques like our “GOPSA” method to ensure accurate cross-domain predictions.

[a] Rodrigues, P. L. C., Jutten, C., & Congedo, M. (2018). Riemannian Procrustes analysis: transfer learning for brain-computer interfaces. IEEE Transactions on Biomedical Engineering, 66(8), 2390-2401.

[e] Sabbagh, D., Ablin, P., Varoquaux, G., Gramfort, A., & Engemann, D. A. (2019). Manifold-regression to predict from MEG/EEG brain signals without source modeling. Advances in Neural Information Processing Systems, 32.

[f] Mellot, A., Collas, A., Rodrigues, P. L., Engemann, D., & Gramfort, A. (2023). Harmonizing and aligning M/EEG datasets with covariance-based techniques to enhance predictive regression modeling. Imaging Neuroscience, 1, 1-23.

[g] Schiratti, J. B., Allassonniere, S., Colliot, O., & Durrleman, S. (2015). Learning spatiotemporal trajectories from manifold-valued longitudinal data. Advances in neural information processing systems, 28.

[h] Li, M., et al (2022). Harmonized-multinational qEEG norms (HarMNqEEG). NeuroImage, 256, 119190.

Dear Reviewer N9Jm,

Thank you for your positive assessment of our submission. In the following, we detailed answers to the weaknesses and questions you raised.

Addressing Weaknesses:

1. Evaluation of the methods on simulated scenarios:
   We agree with the reviewer and added numerical experiments on simulated data to illustrate the advantage of “GOPSA” compared to other methods. All methods are evaluated in 3 simulation scenarios: shift in X only, shift in y only, and joint shift in X and y. The results are presented in Figure 1 of the attached PDF file. We demonstrate that “GOPSA” effectively compensates for the joint shift in both X and y, and is robust to shifts in only X or only y.
2. Interpretability of Riemannian models:
   The interpretation of Riemannian-based models, such as the one we refer to as “No DA,” involves transforming the model’s parameters into interpretable patterns, as demonstrated in prior work [c]. Applying similar techniques to “GOPSA” could yield valuable insights into its decision-making process, revealing the patterns and features that contribute most significantly to its predictions. This analysis is left for future research, and could enhance the transparency and interpretability of our method and guide further refinements and applications.
3. Statistical significance of the results: The reviewer is right that the difference in performance between “GOPSA” and other techniques is not always significant. The previously mentioned analysis done on simulated data showed that “GOPSA” clearly outperforms other techniques in the presence of a joint shift in X and y. Thus, with real data, we expect “GOPSA” to perform better when there is a joint shift in X and y, and to perform similarly to other methods when it is not the case. In Figure 2, which we updated and included in the attached file, we report the p-values of a t-test on the difference “GOPSA” minus ”DO Intercept” for three metrics on all site combinations. We chose to perform this test with “DO Intercept” as it is the competing method with the best performances. “GOPSA” significantly improved Spearman’s rank in 4 out of 5 site combinations. The results for R2 score and MAE are more nuanced, but on the majority of splits both scores tend to be improved with “GOPSA”.

Addressing Questions:

1. Interpretation of the age prediction model:
   As explained previously, we agree with the reviewer that an interpretation analysis of the model would be interesting. In this work we mainly focused on the development of a novel approach to the specific setting of joint shift in both X and y. In our opinion, the specific work required to perform a clean and thorough interpretation represents a contribution in itself, see e.g. [c]. Indeed, the interpretation is not straightforward and would require taking into account many factors: age and frequency are linked, and confounding factors, like the sexe and the mental health of the participants.
2. Statistical test results:
   A statistical t-test was conducted to compare the performance of “GOPSA” with that of the best-performing baseline method, “DO Intercept”. The p-values resulting from this test are presented in Figure 2 and have been updated in the rebuttal PDF. The statistical analysis from Figure 2 ensures a rigorous evaluation of the significance of the performance differences between “GOPSA” and the baselines. “GOPSA” significantly improved Spearman’s rank in 4 out of 5 site combinations. The results for R2 score and MAE are more nuanced, but on the majority of splits both scores tend to be improved with “GOPSA”. Regarding Table 1, we agree that including p-values would be beneficial. However, due to space constraints, we will add them in a subsequent version of the manuscript.
3. Comparison with DL model:
   We added in the benchmark a deep learning (DL) approach, called “GREEN”, that has recently been developed for EEG applications like age prediction [b]. We utilized the “g2” variant, which is an SPD network [d], because the data are SPD matrices. Even though the architecture does not include any adaptation layer, the results on simulated (Figure 1 of the rebuttal PDF) and real data (Figure 2) showed this model to be relatively robust to shifts. On real data “GREEN” performed better than “No DA”, the classical Riemannian-based regression pipeline but it still falls short of the proposed “GOPSA” method which is tailored for the studied problem. However, as mentioned in the conclusion of the paper, since “GOPSA” is implemented using PyTorch, a future perspective would be to integrate “GOPSA” into Riemannian DL models like “GREEN”.
4. Application of the transport to the EEG time series:
   In this work, the dataset we used did not provide the raw EEG time series, so we focused on covariance matrices as features. However, the reviewer’s remark is correct, and the transportation operators learned from the covariance representation could be applied to the EEG time series. We think that it would be interesting to investigate this approach with datasets that provide the raw EEG signals.

[b] Paillard, J., Hipp, J. F., & Engemann, D. A. (2024). GREEN: a lightweight architecture using learnable wavelets and Riemannian geometry for biomarker exploration. bioRxiv, 2024-05.

[c] Kobler, R. J., Hirayama, J. I., Hehenberger, L., Lopes-Dias, C., Müller-Putz, G. R., & Kawanabe, M. (2021, November). On the interpretation of linear Riemannian tangent space model parameters in M/EEG. In 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC) (pp. 5909-5913). IEEE.

[d] Huang, Z., & Van Gool, L. (2017, February). A Riemannian network for SPD matrix learning. In Proceedings of the AAAI conference on artificial intelligence (Vol. 31, No. 1).

Dear Reviewer 85Rg,

Thank you for your positive evaluation of our submission. We took into account your concerns, and answered your questions point by point in the following:

Addressing Weaknesses and Questions:

1. Other similar problem settings and Combat harmonization algorithm:
   The studied problem setting is indeed present in other biomedical applications. For example, you mentioned ComBat which was first developed for gene expression analysis. Other linear mixed-effect models are also usually applied to biomedical data in this context. However, they are often used to harmonize data to perform statistical analysis and thus are not adapted for predictive applications of machine learning. Indeed, to the best of our knowledge, mixed-effects models, like ComBat, assume to have access to the biomedical outcome (y) of source and target domains to compute corrections. In the studied setting, we assume to only have access to the y mean of each target site which is a much harder problem. Furthermore, the proposed benchmark includes the DO Intercept method, which is well suited for machine learning applications, as an alternative to linear mixed-effect models.
2. Overlap of y distributions between source and target domains:
   The second question raises an important point. Theoretically, based on the generative model of the simulated data (as described in the caption of Figure 1 in the rebuttal PDF), the data and outcome (y) are linked by a log-linear relationship. This implies that, knowing the shift in X for the target domain, predictions can be made even when y distributions do not overlap between the source and target. Since GOPSA estimates the target shift in X by minimizing $(\bar{y} - \hat{y}_i)^2$, it is capable of handling such scenarios.
3. Noisy estimate of the target $\bar{y}$ :
   First, we want to point out to the reviewer that we changed how the target $\bar{y}$ is estimated in our empirical benchmark on real data. Indeed, we estimate target $\bar{y}$ on splits (50% of the data) that do not overlap with the evaluation target splits (50% of the data), rather than assuming target $\bar{y}$ to be known. This leads to noisier target $\bar{y}$ estimates. We noticed that the results are similar to the submission’s one. We still note that the performances of 3 splits (out of 100) in site combination 4 (Cu03, M, R, S) have deteriorated. Figure 2 of the rebuttal PDF presents the updated results.
   Our intuition is that the noisier the target $\bar{y}$ estimation is, the worse the performance will get to the point of suppressing the interest in applying a domain-specific adaptation. A future perspective of this work is to study the point at which the adaptation does more harm than good.