I have edited the below question to add more detail:
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?
- 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
The above question seems to answer question 5 in terms of structuring dummy variables for this type of analysis.