# How to analyse multiple temporally correlated response variables?

I measured individuals with two different treatment in a task where they had to find 5 items. All individuals solved the task and found all items.

Here is a dummy dataset:

Individual  Treatment   Lat1    Lat2    Lat3    Lat4    Lat5
1           A           10      17      23      27      32
2           A           5       15      21      32      45
3           A           8       11      22      26      32
4           B           13      25      29      47      61
5           B           17      22      36      55      68
6           B           11      18      28      38      48


I would like to analyse if treatment had an effect on the latency to solve the task. One way to do this would be to just look at the latency to find the 5th item, but I was wondering how I could take into account all latencies at the same time. Of course the 5 latencies are temporally correlated within individual.

Any help is appreciated!

• What question would you hope to answer with all 5 timepoints at the same time? Commented Sep 11, 2016 at 3:52
• If treatment has an effect on the latency to find the items Commented Sep 12, 2016 at 6:08

It's quite easy. It's a question of repeated events. You have two options: that is to look at the effect on the gap time (i.e. for each of the 5 events) and another one is to look at the effect on the intensity of the process, i.e. on the "hazard" driving the sequence of 5 events.

First, I simulated some data in the same layout as yours.

individuals <- 1:20

treatment <- rep(c(1,2), 10)
gap_times <- matrix(rexp(100, rate = rep(treatment, each = 5) ), nrow = 20, ncol = 5, byrow = TRUE)
event_times <- t(apply(gap_times, 1, cumsum))

dat <- data.frame(
Individual = individuals,
Treatment = factor(treatment, labels = c("A", "B")),
event_times
)
names(dat)[3:7] <- c("Lat1", "Lat2", "Lat3", "Lat4", "Lat5")

Individual Treatment       Lat1      Lat2      Lat3     Lat4     Lat5
1          1         A 0.08615526 1.3719387 1.7466040 2.731752 3.102935
2          2         B 0.69924575 1.0118343 1.1675002 1.546387 1.810142
3          3         A 0.86302774 2.0503774 2.5264718 2.818729 4.053160
4          4         B 0.33189293 0.9338124 1.0977293 1.335918 1.351561
5          5         A 2.77120142 3.5721603 4.7291778 5.094909 6.729188
6          6         B 0.34863005 0.4164008 0.8473643 1.000529 2.687727


First thing first, you need your data set in a long format, in a survival time. I do this using the packages tidyr and dplyr:

dat_long1 <- tidyr::gather(dat, key = Lat, value = Time, Lat1:Lat5)
dat_long2 <- dplyr::arrange(dat_long1, Individual, Time)


This gives a data set like this:

> head(dat_long2)
Individual Treatment  Lat       Time
1          1         A Lat1 0.08615526
2          1         A Lat2 1.37193873
3          1         A Lat3 1.74660404
4          1         A Lat4 2.73175214
5          1         A Lat5 3.10293523
6          2         B Lat1 0.69924575


If you want to look at the intensity of the process which generates these you need to create a new variable Time0 (I'm sorry if you are not familiar with dplyr but I'll keep on using it):

dat_long3 <- dat_long2 %>%
group_by(Individual) %>%
mutate(Time0 = c(0, Time[-length(Time)])) %>%
ungroup() %>%
mutate(Status = 1)
# A tibble: 6 × 6
Individual Treatment   Lat       Time      Time0 Status
<int>    <fctr> <chr>      <dbl>      <dbl>  <dbl>
1          1         A  Lat1 0.08615526 0.00000000      1
2          1         A  Lat2 1.37193873 0.08615526      1
3          1         A  Lat3 1.74660404 1.37193873      1
4          1         A  Lat4 2.73175214 1.74660404      1
5          1         A  Lat5 3.10293523 2.73175214      1
6          2         B  Lat1 0.69924575 0.00000000      1


Now you can fit an Andersen-Gill model for the intensity of the recurrent event process. Since you do not have any censoring time, the status variable is always 1.

> mod_ag <- coxph(Surv(Time0, Time, Status) ~ Treatment + cluster(Individual), data = dat_long3)
> mod_ag
Call:
coxph(formula = Surv(Time0, Time, Status) ~ Treatment + cluster(Individual),
data = dat_long3)

coef exp(coef) se(coef) robust se    z       p
TreatmentB 0.642     1.901    0.220     0.156 4.13 3.6e-05

Likelihood ratio test=8.61  on 1 df, p=0.00334
n= 100, number of events= 100


The interpretation in this case is that TreatmentB increases the intensity of the process (the positive sign of the regression coefficient). This regression model basically says that each individual has a counting process N_i(t) which starts at time 0 (that is important) and the intensity of this process is $$\lambda_0(t) \exp(\beta x_i)$$

The second thing you can look at is the gaps between the events. For this I created a gap variable and fitted a model on the gaptimes:

dat_long4 <- dat_long3 %>%
mutate(Gap = Time - Time0)

mod_gap <- coxph(Surv(Gap, Status) ~ Treatment + cluster(Individual), data = dat_long4)
> mod_gap
Call:
coxph(formula = Surv(Gap, Status) ~ Treatment + cluster(Individual),
data = dat_long4)

coef exp(coef) se(coef) robust se    z       p
TreatmentB 0.765     2.148    0.209     0.204 3.74 0.00018

Likelihood ratio test=13.2  on 1 df, p=0.000277
n= 100, number of events= 100


The conclusion is quite similar; individuals with treatment B have shorter gap times. This regression model is $$\gamma_i(w) = \gamma_0(w) \exp(\beta x_i)$$ where $w$ is the time from the previous event.

Which one of the two approaches to use is debatable, but my impression is that since you do not have censoring you can use any of the two. There is a bunch of literature on the topic. In the book of Cook & Lawless Statistical Analysis of Recurrent Events, they say:

Poisson processes tend to be appropriate in settings where events for an individual or system are triggered or influenced by random external factors, whereas renewal processes tend to describe settings in which events flow from physical cycles that are internal to an individual or system. (...) Analyses based on waiting times are often useful when events are relatively infrequent, when some type of individual renewal occurs after an event, or when prediction of the time to the next event is of interest. Analyses based on waiting times are natural in studies of system failures, where repairs made at each failure return the system to a working state

So in the end it's up to you to decide which one is the most relevant one.

• Thanks for your elaborate answer. When I use your code I get the warning message "Warning message: In Surv(latency0, latency, Status) : Stop time must be > start time, NA created". My code is exactly like yours. Do you have any idea what to do? Commented Sep 26, 2016 at 12:48
• Also, in another part of the experiment not all individuals found all items. How should I code that? Thanks for your help! Commented Sep 26, 2016 at 13:29