Nonparametric Covariance Regression for Massive Neural Data on Restricted Covariates via Graph

Authors: Ganchao Wei
Primary area: applications to neuroscience & cognitive science
Keywords: covariance regression, latent variable modeling, Gaussian process, graph Laplacian, time series, spatiotemporal data, neural data

Abstract:
Modern recording techniques enable neuroscientists to simultaneously study neural activity across large populations of neurons, with capturing predictor-dependent correlations being a fundamental challenge in neuroscience. Moreover, the fact that input covariates often lie in restricted subdomains, according to experimental settings, makes inference even more challenging. To address these challenges, we propose a set of nonparametric mean-covariance regression models for high-dimensional neural activity with restricted inputs. These models reduce the dimensionality of neural responses by employing a lower-dimensional latent factor model, where both factor loadings and latent factors are predictor-dependent, to jointly model mean and covariance across covariates. The smoothness of neural activity across experimental conditions is modeled nonparametrically using two Gaussian processes (GPs), applied to both loading basis and latent factors. Additionally, to account for the covariates lying in restricted subspace, we incorporate graph information into the covariance structure. To flexibly infer the model, we use an MCMC algorithm to sample from posterior distributions. After validating and studying the properties of proposed methods by simulations, we apply them to two neural datasets (local field potential and neural spiking data) to demonstrate the usage of models for continuous and counting observations. Overall, the proposed methods provide a framework to jointly model covariate-dependent mean and covariance in high dimensional neural data, especially when the covariates lie in restricted domains. The framework is general and can be easily adapted to various applications beyond neuroscience.

Decision: reject
Reviewer scores: 5, 6, 5, 5

Reviews and discussion:
<thread>
creation date:
1730377839561

thread_type:
official_review

<thread_object>
creator: reviewer_dlgm

<content>
summary:
The authors present a non-parametric covariance regression method that allows modeling non euclidean via the  Graph Laplacian Gaussian process proposed by Dunson et al. 2022. They perform numerical experiments on simulated and neuron data.

strengths:
The paper is sound, correctly written, and well-motivated. Overall, the work reads solid, and the method proposed seems to work relatively well.

weaknesses:
The paper essentially relies on the work of Dunson et al. and the presented work seems to be closer to an adaptation of this work than a novel prior /method in the context of covariance regression. While this is interesting in principle, the method contribution seems a bit limited to a slight variation of the original work of Dunson and colleagues. Similarly, the authors do not tackle any particularly novel complex sampling scheme for their MCMC part and rely on existing work that they blend with the GL-GP prior to Dunson and colleagues.

 Further, the simulations are somehow limited (perhaps due to the heaviness of the computation for inference).  

Figure 6 seems to advocate the most for the method, but Figures 1, 4, and 5 do not strike me as proper improvements.
Lastly, it would have been nice to compare with non-GL-GP methods (perhaps a standard vanilla GP approach).

questions:
Have run similar experiment using some of the competing methods?
Could you clarify the data size (covariate and neurons data).
</content>
</thread_object>

</thread>
<thread>
creation date:
1730592638996

thread_type:
official_review

<thread_object>
creator: reviewer_1hvp

<content>
summary:
The authors present a framework for modeling neural data and covariates through including graph properties of covariates to account for restricted subspaces. They demonstrate their method in simulation and in two experimental datasets.

strengths:
- The inclusion of accounting for restricted space / inputs is smart and very relevant for neuroscience applications in particular. 
- The model components and math are well explained individually, though the thread of motivation throughout the methods section could be strengthened to connect back to the overall picture.  
- Simulations, training, selection of parameters, and results on experimental datasets are all very clear with a high level of detail (needed for any future work or replication, which is excellent).

weaknesses:
- The results overall are underwhelming. In the simulated data, all model fits look extremely similar. How much of the observed variation would be needed in a neuroscience context? In the fits to experimental data, the significance of the results is not explained. The PC plots are results, but what do they tell us? What new findings did we learn from this method of modeling the neural data? 
- Computational costs are discussed but not shown empirically.
- Minor notes: typo (extra ')') in line 66 
- References to e.g. Figure 4 (which is really in the supplement) are confusing. 
- Grammar throughout could use some editing for correctness and clarity (e.g., the sentence spanning lines 82-84 doesn't make sense to me). 
- The provided code is complex (many files/functions) and lacks documentation, and is in Matlab, which lessens its overall impact in the ML field (and likely the efficiency and speed of computation for MCMC).

questions:
- The code implementation appeared custom. Why not use established packages where appropriate (e.g. MCMC)?
- "we observe that the inference can be sensitive to hyper-parameters ({, K, t}) for GL-GP" (line 222) Is there evidence to show here?
- How does this scale to higher numbers of latent factors and increased dimensionality?  The simpler experiments showed "L = 10 is large enough" (line 262), so perhaps it is not a concern for low-dimensional tasks, but I would be curious to know the authors' thoughts on scalability (to tasks involving more freely moving behaviors, for example). 
- The observed data for the simulations in Figure 1 were fairly widely spread out throughout the space. How well would this method work for the cases where test or new data extended outside the range of previous observations?
</content>
</thread_object>

</thread>
<thread>
creation date:
1730679453938

thread_type:
official_review

<thread_object>
creator: reviewer_wnz7

<content>
summary:
This work proposes a nonparametric mean-covariance regression model that leverages latent factor analysis and Gaussian Processes (GPs) to jointly model the mean and covariance of the high-dimensional neural data. Despite the a bit hard to follow writing flow, a notable claimed contribution is the integration of graph-based Gaussian processes to manage covariates in resricted subspaces with the MCMC algorithm. The model is validated through simulation dataset and applied to real-world LFP and HC neural datasets.

strengths:
1. The introduction of graph-based Gaussian processes to handle restricted covariates adds a novel dimension to mean-covariance modeling, which I believe is helpful for neural data. Because this approach can be effective while maintaining the interpretability of the method.
2. The authors employ an MCMC algorithm with data augmentation techniques to handle intractable count models effectively, which can be more computationally effficient.
3. The authors argued and claimed that this proposed framework can be extended beyond neuroscience resaerch to other domains requiring mean-covariance modeling under similar constraints.

weaknesses:
1. The paper is a bit hard to follow in writing, the scientific question or mathematical task you hope to solve is not very clear.
2. There has the assumption of independent GPs for each factor dimension, although computationally efficient, might overlook some inherenet underlying dependencies between factors that could be relevant in most neuroscience applications.

questions:
1. What's the potential broader applications of this framework beyond neuroscience?
2. My other concerns please relate to the Weakness section.
</content>
</thread_object>

</thread>
<thread>
creation date:
1731029305742

thread_type:
official_review

<thread_object>
creator: reviewer_4jkz

<content>
summary:
The authors consider modeling neuroscience data (such as LFP and neural spike count data) in settings where experimental conditions can vary in complicated ways (that may not be well-modeled by assuming conditions lie in a Euclidean space, for instance).  They propose a non-parametric mean-covariance regression model.   They impose smoothness along both experimentation conditions and latent factors. To handle restricted domains they incorporate Graph Laplacian (GL) kernels.  For count data they handle non-conjugate count likelihood using a data augmentation technique.  In the inference step, they propose an MCMC algorithm with a data augmentation technique and derive full conditional to get posterior samples. They consider both continuous (LFP) and count observation real world neuroscience datasets.

strengths:
## Strengths
- Analyzing neuroscience data (esp. count also LFP) is challenging but important.  
- The problem of handling restricted spaces arising from experimental conditions is well-motivated.
- The use of a graph-Laplacian GP to allow for more flexibility than assuming an unknown manifold is interesting.
- The presentation overall was good.
- Experiments are included from two real-world neuroscience experiments, one on LFP data and one for count data.

weaknesses:
## Weaknesses
### Major
- My major concern regards technical novelty.  From my reading, the results seem to largely follow from integrating two key prior works, Fox and Dunson 2015 which proposed  predictor-dependent factor loadings modeled using GPs and Dunson et al 2021 which studied Graph Laplacian based GP regression for restricted domains.  Perhaps the authors could elaborate on any significant technical challenges that were overcome.
- The title includes “MASSIVE NEURAL DATA” and “massive” is mentioned several times but not described in particular with respect to computational or sample complexities
    – can you elaborate what ranges of which data dimensions qualify as “massive”?  
    - The computational complexity of the proposed method is not discussed, especially important in settings with massive data; how does the complexity vary with different data dimensions?
    -  Relatedly, empirical run-times for your experiments and the GPWP baseline are not reported
    - sample complexity is not analyzed theoretically or experimentally
    - The HC experiment involved 36 neurons in one recording session  with 200 ms bins.  The LFP data set I am not sure of the dimensions (eg LFPs from 13 areas – how many LFPs recorded per area?), but from the description did not appear to be ``massive.’’  
    - (minor) in the discussion, a point is made regarding computational complexity in lines 467-470 “in applications to massive data” which sounds like saying in settings with even larger data dimensions than is considered here


### Minor
- For the GPWP baseline (Nejatbakhsh et al., 2023), the authors mention the used "single trial", but did not include further discussion on the choice.  
- The GPWP reported value in the simulation experiment for the Gaussian case differs dramatically – any thoughts why that is?
- The L-GP baseline’s scores in most experiments (to me) seemed pretty close to the proposed method. 


### Very minor (Typos/etc.)
- line 71 ‘massive neurons.’
- 081 "(p >> N)" although standard notations, p and N are have not been introduced before this line
- a few commas start a new line after equation blocks (eg line 148)
- line 196 ‘identifiablility.’

questions:
(I have some doubts related to technical novelty, computational complexity, sample complexity, and baselines mentioned above -- clarification on those points would be appreciated)
</content>
</thread_object>

</thread>
<thread>
creation date:
1731939880672

thread_type:
official_comment

<thread_object>
creator: authors

<content>
title:
General Comments

comment:
Thanks all the reviewers for their constructive comments, and we have updated the manuscript accordingly. Here, we try to clarify the following three points raised by most reviewers.

## 1. Novelty

Overall, the primary purpose of this paper is to introduce a new method and perspective for addressing neuroscience problems for large number of neurons/ subjects. Previously, researchers mainly focused on dimension reduction (e.g., LDS and GPFA), which can at most parse to mean regression problem. Meanwhile, many mean-covariance regression methods do not account for dimension reduction. In addition to merging these two perspectives, we incorporate the geometric information of restricted inputs to enhance performance.

Besides, another significant contribution is providing a tractable inference method for count data (e.g., neural spikes) via P\'olya-Gamma (PG) augmentation. Furthermore, here are some additional points we'd like to highlight. 1) Combining building blocks is non-trivial for neural data analysis, as brute-force inference becomes infeasible for lengthy recordings. Here, we sequentially sample $\xi(x)$. 2) We also observed that while all models show relative robustness for Gaussian responses, the fitting results for count responses are sensitive to the latent dimension $k$ in L-GP and L-GLGP-fixed, while L-GLGP-adaptive model is not.This empirically suggests the need to consider all effects in count data analysis.

## 2. Computational Complexity

The computational xomplexity was added to the updated manuscript. Basically, all the following results are from a laptop with an Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz 2.11 GHz.

In simulations, the training response for both continuous and counting cases has $N = 50$ and $p = 100$, and use $q=2$ covariates. The following two tables show the time consumed for each iteration, under different response types and latent dimensions. The "-fixed" means the hyperparameters for kernel function is fixed, while "-adaptive" means they are sampled in MCMC.

For continuous response:

|  | L-GP-fixed | L-GP-adaptive | L-GLGP-fixed| L-GLGP-adaptive|
|----------|----------|----------|----------|----------|
| $k=2, L = 10$ | 0.12s | 0.20s | 0.13s | 0.21s |
| $k=5, L = 10$ | 0.18s | 0.28s | 0.22s | 0.25s |


For counting response:

|  | L-GP-fixed | L-GP-adaptive | L-GLGP-fixed| L-GLGP-adaptive|
|----------|----------|----------|----------|----------|
| $k=2, L = 10$ | 0.31s | 0.38s | 0.32s | 0.75s |
| $k=5, L = 10$ | 0.37s | 0.45s | 0.37s | 0.75s |

In terms of application to neural data: 1) For LFP with a continuous response, the training response matrix is $14 \times 700$, and the covariate dimension is 3. Each iteration takes ~3.5 seconds. 2) For HC with a counting response, the training response matrix is $36 \times 48$, and the covariate dimension is 2. Each iteration takes ~3.3 seconds. 

The number of neurons does not significantly impact computation (as in simulations), because we use a latent factor model for dimensionality reduction. At present, the main bottleneck for computation time is the number of experimental conditions (recording length for time series), as we require large matrix inversions for the GP. For application to data with long recording/ large number of experiment conditions we suggest to do approximated inference (e.g. variational inference) or use some trick to reduce computation cost for GP matrix inversion (e.g. sparse GP or [1]), as mentioned in the discussion. 

## 3. Performance among different methods

All latent models with smoothness,by (GL)GP, over different conditions (L-models) are somewhat similar in terms of overall metrics such as held-out log-likelihood (although in LFP application, L-GLGP-adaptive improves more significantly). 

However, 1) these L-models are significantly better than non latent models, such as GPWP and dCMP; 2) a factor model without smoothness (e.g. by GP) is consistent only when $p/N\to 0$ ([2]). These suggest that considering dimension reudction and smoothness (to some degree, reduce effective $p$) over experiement conditions/ time can handle most cases, and the overall performance of GP and GLGP are similar (this also applies to the original GLGP paper, Dunson et al 2022). 

**However**, some local performance (estimation at specific conditions/ time points/..., e.g. especially in figure A3) can be quite different for different models, and for problems requiring better local inference, considering geoemtric information is necessary. To make things more visually significant, we can design a more complicated geometric constraint.

[1] Radial neighbours for provably accurate scalable approximations of Gaussian processes, Zhu et al., Biometrika, 2024

[2] On Consistency and Sparsity for Principal Components Analysis in High Dimensions, Johnstone I.M. & Lu, A.Y., Journal of the American Statistical Association, 2009
</content>
</thread_object>

</thread>


Metareview:
metareview:
This paper proposes a non-parametric mean-covariance regression model for neural data analysis, integrating Graph Laplacian Gaussian Processes (GL-GPs) to address restricted experimental subspaces. Validated on simulated datasets and real-world neural data (LFP and spike counts), the method employs MCMC-based inference, showcasing versatility across data types.

Despite this, the paper suffers from poor clarity in notation, insufficient justification for key modeling choices, and overstated claims regarding scalability for "massive neural data." The further factorization of latent components is not well-motivated, and the kernel definitions for Gaussian processes lack transparency, leaving computational complexity underexplored. Additionally, the scalability of the inference process is questionable, given the challenges of GP-based methods when handling large datasets or time bins. These factors, coupled with unclear writing and underwhelming results, diminish the paper’s overall impact. I recommend rejection due to the limited technical novelty and these critical issues.