I am interested in computing the zero-order cross correlation of two time series ($y$ and $z$) that represent neural activity (fMRI) in distinct regions. fMRI data tend to have a strong autoregressive structure and in my data, the errors are often fit by something in the ARMA(3,2) range. Given this temporal structure, I am interested in prewhitening the time series, then computing their correlation.

Having read similar posts and resources (e.g., https://onlinecourses.science.psu.edu/stat510/node/75/), I have learned that to prewhiten, one can fit an ARIMA model to one series, apply that model (with its fitted coefficients) to both series, then compute the cross-correlation of their residuals. This approach is implemented for example, in the prewhiten function of the TSA package in R. Likewise, I've been able to implement this using the Arima function from the forecast package, passing model as an argument when fitting the second series.

The dilemma is that I have a set of nuisance time series that I want to remove from the data at the same time. For example, the mean time series in cerebrospinal fluid can't be neural activity by definition, so I would like to remove it from both neural time series. I don't have any conceptual interest in the temporal structure of the nuisance time series, so a multivariate model like VAR seems like overkill for this purpose. But the nuisance series also have an autoregressive structure.

In short, how does one prewhiten two time series in the presence of nuisance covariates that could affect either one to varying degrees?

My sense is that one could include the nuisance series as exogenous covariates ($X$) when fitting the first time series of interest ($y$) in a regression with ARMA errors (as discussed here: https://robjhyndman.com/hyndsight/arimax/):

$$y_t = \beta_1 x_{1,t} + \beta_2 x_{2,t} + \eta_t$$ $$\eta_t = \phi_1 \eta_{t-1} + ... + \phi_p \eta_{t-p} - \theta_1 z_{t-1} - ... - \theta_q z_{t-q} + z_t$$

This seems tractable with the xreg argument of arima in R:

data(lh) #first time series of interest lh2 <- .5*lh + .5*lh^2 + 0.05*1:length(lh) #second time series of interest exog_x <- matrix(rnorm(96), ncol=2) #nuisance covariates m1 <- forecast::Arima(lh, order=c(2, 0, 1), xreg=exog_x) m2 <- forecast::Arima(lh2, model=m1, xreg=exog_x) cor(resid(m1), resid(m2)) #zero-order cross-correlation on whitened residuals

Is this a reasonable approach -- essentially passing the same xreg to both ARIMA models? My concern is that the nuisance variable may be associated with the time series of interest to different degrees (e.g., correlating with $y$ at $r = 0.5$, but $z$ at $r = 0.1$). For such reasons, people in my field often regress nuisance signals out in a separate step, allowing for unique fits for each region of interest. The residuals of these are then used for cross-correlations. It seems somewhat unprincipled to have two regression models, however, given that these are nuisance signals that should probably be fit at the same time as prewhitening.

One alternative I've considered is freeing the coefficients for the covariates when fitting the second series while fixing the ARIMA coefficients from the first model. Thanks in advance for any suggestions or resources for how to handle prewhitening in the presence of nuisance exogenous covariates.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.