Interrupted Time Series Analysis - ARIMAX for High Frequency Biological Data? [closed]

Question

I have edited the below question to add more detail:

The Problem

I am currently working on doing an analysis on fluorescence data acquired from mice performing a behavioral task. As the data is acquired, specific time points are marked by external cues and specific actions performed by the mouse. Below is an example of the type of data we acquire:

My main goal is to describe the magnitude and the direction of the effects of each of the external events on signal amplitude. From my reading, it seems that the best way to go about doing this, when I consider the nature of the data as shown above, is to fit an ARIMA model to my data and to consider each of those external events an intervention in my model. Figure 1 above is a snippet (about 7 seconds) of a much longer recording (>75 min):

As you can see, over the course of a session, there is a downward shift in the mean signal - caused by bleaching of the signal. I account for that bleaching and readjust the values in Figure 3:

This is the data that I perform subsequent analyses on. Initially the data was sampled at 1000Hz, but I down-sampled the data to 100 Hz. We lose a lot of information, but for our experiment 100Hz is probably more appropriate than the previous and it makes analysis go much quicker.

Question 0: When I import my time series into R, what frequency do you think I should use? I have been leaving it as its default of 1, but when I try to change the frequency to match the data acquisiton speed ts(normDat, frequency = 100), all of my subsequent calculations take FOREVER. Like they don't finish after 2 hours. There is no inherent frequency in the data acuqision as there is with daily, weekly or yearly data.

As I mentioned, Figure 1 is merely a snapshot of the data - a represenation of a single "trial". In the course of a 75 minute session, our mice perform hundreds of these trials, so I can either segment the long trace into individual "trials" from single sessions, or consider the entire session as one unit. As I mentioned, my main question is how I can quantify the immediate and transient effects of each intervention on the signal profile. Let us take the "Consumption" Event as an example. When observing the experiment I noticed an unambiguous decrease in the signal every time the mouse was given a reward. Let us say that the mouse that produced the above traces did 300 trials and thus had 300 of those consumption events in a single session. If I align the consumption events and average the signals together to eliminate sporadic noise, the characteristic response pattern to the intervention that I observed as the animal was behaving reveals itself:

As you can see, once the event registers, there is a short delay followed by an unambiguous downward trend. I chose to fit an ARIMAX model to my data in order to account for the fact that the signal we see at any given moment is influenced by past events (that is, the signal isn't a perfect representation of internal state, but a representation of internal state contaminated by previous internal state). My thought was that this model would be most appropriate because it would let me take into account not just the lagged effects of previous signals, but also the intervention effects of external events given those lagged effects. You will note for instance that in figure 1, there are 5 external events that I am interested in profiling and each has a distinct profile as in figure 4. The time between the different external events varies from trial to trial, though they occur in the same order on every trial. My two main concerns in analyzing the data using a simpler method than the ARIMAX is that 1)the effects of adjecent events may overlap in either an opposing or additive fashion 2)because of trial to trial variability I can't simply look at the amplitude of the signal -for example, the consumption event in figure a starts when y = 7, but you can see that in figure 4 the amplitude of the averaged peak is < 3. Because of this variability I thought taking into account previous signal history and looking for impulse and transitory effects was the most ideal model.

Question 1: My first question is whether this approach seems reasonable or if there is a more straightforward way to get at the answer to the intervention effects. Basically, is my approach the question presented appropriate.

Originally, I had been using Matlab to try to identify an appropriate arima model to fit my data (for the rest of this post, I analyze the entire session as a whole, not individual trials - if anyone thinks doing this trial by trial is smarter, let me know!) but it seems most people seem to prefer R for this process. I thought the first step in developing my model would be to draw a correlogram for the entire session, but I quickly realized that because of the high rate of sampling, the correlogram basically showed a correlation of 1 out to about a large lag.

These data seemed to confirm my suspicion of an ARIMA process so I ran an auto.arima first (without exogenous regressors):

model1<-auto.arima(normDat, max.p = 6, max.q=6, seasonal = TRUE, stepwise = FALSE, parallel = TRUE, allowmean = TRUE, allowdrift = TRUE)

[Because I'm a noob, this first time I ran this I set max.p and q = 100. That was fun]

The output I got was: ARIMA(4,0,1) with non-zero mean

I was a little bit confused by the results for a number of reasons:

Question 2: Why would the auto arima return a model with no differencing when you consider the correlogram of the raw data? I thought a graph like that suggested non-stationarity? For what its worth I ran a KPSS test myself and found that the p value was >.05 so we couldn't reject the null of stationarity.

And a more general question:

Question 3: What is the relationship between lagged auto correlations and the order of the arima produced. Why, for instance is the AR order of the best model 4 when the data is autocorrelated out to like 100 lags.

When I graphed the autocorrelation of the residuals:

I believe that this graph shows seasonality that should be accounted for in the auto arima (I know that I haven't accounted for the external interventions yet; the sequence of external events produces a cyclic effect on short time scales). I set D=1 to try to force auto.arima to include a seasonal component, but the output is the same no matter how many times I tried. I then tried to manually add a seasonal effect (0,0,1) vs. (0,1,0) vs. (1,0,0). None seemed to do the trick. It didn't help that I'm unclear about how the seasonal errors I observed in the acf of the residuals should inform my seasonal period. It looks like the seasonal period is about 3. So I have been using that. Didn't help

Question 4: Does it seem like I need to account for seasonality before moving onto my intervention effects? Is there a systematic way of accounting for seasonality?

And finally I took the readout of the auto.arima (just trusting it until I find a better solution):

fit <- arimax(normDat, order = c(4, 0, 1), 
          xtransf = data.frame(pulse = consumption$X0, 
                               transit = consumption$X0),
          transfer = list(c(0, 0), c(1, 0)))

pulse (the immediate signal effect) and transit (effect over time) are identical. They are 0 at every timepoint except for the timestamp at which the consumption events occur. That means that if the mouse consumed 300 rewards, there are 300 1's buried in those vast data stuctures. As I mentioned, for now I am treating the data from all of the trials together so there are multiple pulses for this type of intervention.

I continually received this error:

Error in optim(init[mask], armafn, method = "BFGS", hessian = TRUE, control = optim.control, : non-finite finite-difference value 2

I tried a lot of solutions using optim.control, trying to readjust the method used or the hessian to FALSE, but no dice. I worry that this might be caused by the way my dummies are structured.

Question 5: What is the correct way to structure my dummy variable to answer the question as I have posed it above.

1.Should I only have one 1 corresponding to each event to model a pulse as I have done? Or since the time scale of my recordings are so small, should I have a group of 1's model my pulse?

2. Is it more appropriate to semgent the data into within session trials as I have previously alluded to and then have dummies that are 0 before intervention and 1 after? My concern was that this would model a step function and these interventions are not step functions.

Question 6: Can the arimax function deal with multiple types of interventions (each of which occurs multiple times [i.e. consumption occurs 300 times, choice occurs 300 times etc...]) at the same times. Or should I parallelize the analysis (model each type of intervention in its own arimax and possible do the same for each trial) and average everything together at the end

This is an obnoxiously long question but any help is welcome! I've beeen working on this exclusively for a week and my understanding has jumped greatly in that short time. I've only been working on time series a short time so please be gentle =]

Please let me know if there is anything that might be helpful in answering the question or if i've made an error in logic/the format of my question

Edit:

Intervention Analysis Coding in R TSA Package

The above question seems to answer question 5 in terms of structuring dummy variables for this type of analysis.

Answer 1

My first question is whether [an ARIMA] seems reasonable or if there is a more straightforward way to get at the answer to the intervention effects.

If there are seasonal effects or time-treatment interactions, you cannot naively employ an ARIMA. However, the basic approach of using a linear regression with an autoregressive term is a well accepted modeling strategy for time series data.

should I be concerned about abrogating the effect size of my external events by seasonally differencing the data?

Quite the opposite, in fact. Deseasoning the data will improve the efficiency of estimates of an intervention effect at a minimum. And, due to the imbalance of the intervention administration times, alleviate any possible bias. The correlogram should be calculated from residuals using a model that controls for intervention administration, otherwise the intervention effects are taken to be Gaussian noise, underestimating the actual autoregressive effect.

The periodicity of a time series should really not be evaluated under visual inspection. Again, if you are able to fit a model which subtracts off the intervention effects (as well as random effects due to repeated measures within each mouse), the residuals can be used to inspect periodicity using an FFT or (even better) you as an analyst can simply state the periodicity is X. I know that mice are relatively less predictable in their circadian and other rhythms. Modeling blood pressure in humans, however, was easily a matter of stating that we follow a 24 hours cycle (in vivo). An FFT can complement your researcher perspective in that matter.

Do you think it is more appropriate to treat each trial independently, segment the data and then find a way to average across trials?!

You should absolutely do this. What constitutes a trial may not be clear here. If mice were measured continuously, the best structure would be to consider each day (or appropriate rhythmic block) an observation. Then days are nested within mice, and replicates within each mouse can be handled using random intercepts, seasonal trends can be handled with fixed effects (like a spline or similar), and intervention effects are time-varying covariates.

Is it good practice to difference the entire data set anyway, even if I end up treating the trials independently?

What you are describing sounds like a learning effect (rarely actually "learning" in a strict sense, but a tendency for experimental time to be a lurking factor). "Differencing the entire dataset" doesn't make sense to me, but you can run more complex models adjusting for experimental time as a second time, and its interaction with intervention, to determine whether there is a learning or growth effect.