8
$\begingroup$

Description: Temperature if six objects (obj_A, obj_B,...obj_F) was measured in 1 hour intervals (1,2...10) Objects were under the influence of two treatments (A and B). Treatment A = obj_A, obj_B, obj_C; treatment B = obj_D, obj_E, obj_F.

Problem is that measured values of each object are serially correlated, therefore I can not use classical ANOVA. How to take into account such fact?

# example data
my.data <- data.frame(object = rep(c("obj_A","obj_B","obj_C",
                                     "obj_D","obj_E","obj_F"),
                                      each = 10),
                      time = rep(c(1:10), times = 6),
                      treatment = rep(c("A","B"), each = 30),
                      value = c(4,4,7,8,8,10,8,12,14,12,
                                8,8,12,12,10,12,10,11,12,16,
                                12,12,11,13,12,16,16,14,16,20,
                                11,20,23,27,31,29,31,32,28,30,
                                12,16,17,23,22,24,33,31,31,32,
                                14,13,19,20,24,26,24,28,25,23))

# converting values to time series object
obj_A <- ts(my.data$value[my.data$object=="obj_A"],
            start = 1, end = 10, frequency = 1)
obj_B <- ts(my.data$value[my.data$object=="obj_B"],
            start = 1, end = 10, frequency = 1)
obj_C <- ts(my.data$value[my.data$object=="obj_C"],
            start = 1, end = 10, frequency = 1)
obj_D <- ts(my.data$value[my.data$object=="obj_D"],
            start = 1, end = 10, frequency = 1)
obj_E <- ts(my.data$value[my.data$object=="obj_E"],
            start = 1, end = 10, frequency = 1)
obj_F <- ts(my.data$value[my.data$object=="obj_F"],
            start = 1, end = 10, frequency = 1)

# plot -> blue = treatment A; red = treatment B
ts.plot(obj_A, obj_B, obj_C, obj_D, obj_E, obj_F,
        col=c("deepskyblue","deepskyblue1","deepskyblue2",
              "darkred","indianred","indianred1"),
        lwd = 2.5, lty = 2, xlab = "time", ylab = "temperature")

enter image description here

How to rigorously test whether temperature of objects differ according to used treatment, but without ignoring serial correlation?

Improve this question
$\endgroup$
3
  • $\begingroup$ You can't do it here because there is no one clear way to define if they "differ". Are you interested in the trajectory, or the final value after a set period of time, or something else? $\endgroup$
    – shadowtalker
    Commented Feb 23, 2016 at 10:58
  • $\begingroup$ In your case a t-test even on a single time-point will show significant difference between groups (for several different time points). $\endgroup$
    – amoeba
    Commented Feb 23, 2016 at 11:28
  • $\begingroup$ Hi @ssdecontrol and @amoeba! Let say that each curve is generated by some hidden process. Is it possible to test whether "red curves" are generated by different process than "blue curves"? $\endgroup$
    – Ladislav Naďo
    Commented Feb 23, 2016 at 12:45

2 Answers 2

Reset to default
8
$\begingroup$

I believe your problem is a textbook case of longitudinal data analysis for which there is an extensive statistical methodology developed. In particular, analysis of response profiles is an approach that I would employ in this case.

To motivate this, I will provide some background from Applied Longitudinal Analysis by G.M. Fitzmaurice:

1.2 Longitudinal and clustered data

The defining feature of longitudinal studies is that measurements of the same individuals are taken repeatedly through time, thereby allowing the direct study of change over the time. The primary goal of a longitudinal study is to characterize the change in response over the time and the factors that influence change. (...)

A distinctive feature of longitudinal data is that they are clustered. In longitudinal studies the clusters are composed of the repeated measurements obtained from a single individual at different occasions. Observations within a cluster will typically exhibit positive correlation, and this correlation must be accounted for in the analysis. Longitudinal data also have a temporal order; the first measurement within a cluster necessarily comes before the second measurement, and so on. The ordering of the repeated measures has important implications for analysis.

and regarding analyzing response profiles choice:

Methods for analyzing response profiles are appealing when there is a single categorical covariate (perhaps denoting different treatment or exposure groups) and when no specific a priori pattern for the differences in the response profiles between groups can be specified. When repeated measures are obtained at the same sequence of occasions, the data can be summarized by the estimated mean response at each occasion, stratified by levels of a group factor. At any given level of the group factor, the sequence of means over time is referred to as the mean response profile.

The nlme R library contains a glm function to evaluate such model, giving you a longitudinal "equivalent" of Anova:

library(nlme)

# Model assuming the same variance for each time point
gls.fit <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit)

# Model allowing for different variance structure for each time point 
gls.fit.diff.var <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      weights = varIdent(form = ~ 1 | factor(time)),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit.diff.var)

Unfortunately, in this case model estimation does not coverage (even though I changed control parameters for more convenient). I am afraid there is just too few cases in your data set to estimate parameters.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Why didn't you try corr = corSymm(from= ~1 | treatment)? I assume that among the groups with the same treatment, all experimental units are identically distributed, so have the same covariance matrix. Yet it still won't converge: You have more parameters in your 10x10 covariance matrices than experimental units. But you can use corr = corAR1(from=~1|treatment). Then it will converge, becaus you'll only have two parameters each. But is it reasonably? It assumes some kind of equal spacing of the time points wrt. correlation; this "ordering of the repeated measurements" from your reference. $\endgroup$
    – Horst Grünbusch
    Commented Feb 26, 2016 at 22:14
  • $\begingroup$ I assume that data is correlated within each of the subjects (values for the same subject are correlated thorough time), thus I define correlation structure in relation to a subject (here: object). $\endgroup$
    – Marta Karas
    Commented Feb 27, 2016 at 4:36
  • $\begingroup$ "I assume that among the groups with the same treatment, all experimental units are identically distributed" <- please note we are not talking about a cross-sectional study (data observed in one particular time point); our concern here is not assuming different settling for subjects from different groups at given time point, but assuming that thorough the study tim, observations within a particular subject are correlated (and there is nothing in this point about being from different / same treatment groups). $\endgroup$
    – Marta Karas
    Commented Feb 27, 2016 at 4:49
  • $\begingroup$ OK, it was misleading terminology. "Group" usually means the factor the independent replications are nested within (i.e. samples). The subjects within the same sample are iid. For the corSymm help, a "grouping factor" is the replication itself. $\endgroup$
    – Horst Grünbusch
    Commented Feb 28, 2016 at 13:11
  • $\begingroup$ Dear @MartaKaras I am very thankful for your answer, however solution given by HorstGrünbusch seem to be more suitable for my data (data provided in example are fabricated). Thank you Marta. $\endgroup$
    – Ladislav Naďo
    Commented Feb 29, 2016 at 14:15
9
+50
$\begingroup$

You notice correcly that usual ANOVA cannot handle this type of very heteroskedastic and highly dependent time point data. But also, other repeated measures or multivariate procedures would fail because you have only 6 replications but 10 time points, so called high dimensional data. In the end, it is proven that there cannot be an exact test for data like that.

However, there is a good approximative test, it's the Huynh-Feldt-procedure, see this technical report for references and for a generalization to possibly unequal covariance matrices that you may use. It is defined even for your very small data set. There are no assumptions on the variances and covariances. It is intended for normally distributed data, but I think you're fine with it.

You can test if the curves differ (set T = diag(10)) and if there is an interaction between treatment and time point (set T = diag(10) - 1/10). It's too simple for an extra R-package. Check this: 

Tmat <- diag(10) # or Tmat <- diag(10)-1/10
Y1 <- t(matrix(my.data[my.data$treatment=="B","value"],10,3)) %*% Tmat
    Y2 <- t(matrix(my.data[my.data$treatment=="A","value"],10,3)) %*% Tmat

Mean1 <- apply(Y1,2,mean)
Mean2 <- apply(Y2,2,mean)
S1 <- cov(Y1)
S2 <- cov(Y2)

F <- sum((Mean1-Mean2)**2) / (sum(diag(S1/3 + S2/3))) # see eq. 3.18
trS1sq <- 3/2 * ( (sum(diag(S1)))**2 - 2/3 * sum(S1**2)) # see eq. 3.26
trS2sq <- 3/2 * ( (sum(diag(S2)))**2 - 2/3 * sum(S2**2)) 
trS1S1 <- 3/2 * ( sum(S1**2) - (sum(diag(S1)))**2/2)  # 3.27
trS2S2 <- 3/2 * ( sum(S2**2) - (sum(diag(S2)))**2/2)  # 3.27

f <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1+trS2S2+2*sum(S1*S2))  # see p. 16
f0 <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1/2 + trS2S2/2)

p.value <- 1-pf(F,f,f0)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.