6
$\begingroup$

I feel like I'm missing something obvious, but here we go. I have auto-correlated data measured in triplicate for two (or more) treatments. Something like this:

t <- 3:20 #times in my real dataset are possibly not always equidistant
a <- structure(c(0.652492388457625, 0.905172522010166, 1.23437705454616, 
                 1.48003667490842, 1.77876898946135, 1.99175317367897, 2.31666502140984, 
                 2.43520651415548, 2.67903421794922, 2.84115747823017, 2.89693734873647, 
                 2.91199679761145, 2.85645436179354, 2.99371033437697, 2.99965220711105, 
                 2.84984814715963, 2.64275376547326, 2.64060469520379, 0.481029734912324, 
                 0.8466803252367, 1.31126162780809, 1.56745630574946, 1.74865844658142, 
                 1.80367117155375, 2.06688393210808, 2.24500095501872, 2.52978288460243, 
                 2.69073206006205, 2.89657418056785, 2.93759772556246, 2.99305951550274, 
                 2.89146932307489, 2.88890777189028, 2.7974672802907, 2.70933381639295, 
                 2.66799551352975, 0.624178180970784, 0.867127935268765, 1.09752295578438, 
                 1.35037796202753, 1.60094288950107, 1.97949255710341, 2.15496378191076, 
                 2.42556913246041, 2.54331160179646, 2.67440414122285, 2.84249532365163, 
                 2.95278639560433, 3.06192227561515, 3.03297885461444, 3.04101341059534, 
                 3.01736966686846, 2.80061410999215, 2.69852643323913), 
               .Dim = c(18L, 3L), .Dimnames = list(NULL, c("a1", "a2", "a3")))
b <- structure(c(0.516527990622755, 0.84883434472028, 1.04202664437099, 
                 1.3100841689546, 1.48050413266838, 1.7824492800856, 1.96557179831706, 
                 2.17419105778186, 2.2453178060978, 2.35460428313729, 2.49308342865959, 
                 2.62343038370418, 2.70831189685371, 2.79459971623943, 2.94938536147398, 
                 3.04822554887815, 3.00287042052314, 2.91673487674283, 0.589490441973075, 
                 0.751768045201717, 0.917973959434798, 1.17617337222852, 1.39497560590896, 
                 1.65920945485901, 1.87749014780468, 2.11880355292648, 2.372755207219, 
                 2.46211141942227, 2.59688733749884, 2.72270421752644, 2.79848710425447, 
                 2.81134394947587, 2.75390203306788, 2.78499114431362, 2.86001341271914, 
                 2.95652300178809, 0.558662398944567, 0.834996005844121, 0.988238211915554, 
                 1.27569591423003, 1.38577342414377, 1.62664982549252, 1.83299700801392, 
                 2.04943560731628, 2.22950648854987, 2.38533269800646, 2.49845003387994, 
                 2.60036098089373, 2.61941602504858, 2.71298500309883, 2.78126388719353, 
                 3.04792375845498, 3.02691814463875, 3.06667590650438), 
               .Dim = c(18L, 3L), .Dimnames = list(NULL, c("b1", "b2", "b3")))

matplot(t,a,pch=1,xlab="",ylab="",col="blue")
matlines(t,a,col="blue", lty=2)

matpoints(t,b,pch=16,col="red")
matlines(t,b,col="red", lty=2)

simulated auto-correlated data

I would like to know in which time periods the treatments differ. I would like to avoid fitting any kind of model. (There are models for my kind of data from science, but they are known to be only an approximation for some ranges of my data and I'm afraid that model error might mask differences.) My idea is to calculate the mean and construct confidence intervals (using an assumption of normality) like this:

a_means <- apply(a,1,mean)
a_sds <- apply(a,1,sd)
a_lwr <- a_means-qt(0.975,3)*a_sds/sqrt(3)
a_upr <- a_means+qt(0.975,3)*a_sds/sqrt(3)

b_means <- apply(b,1,mean)
b_sds <- apply(b,1,sd)
b_lwr <- b_means-qt(0.975,3)*b_sds/sqrt(3)
b_upr <- b_means+qt(0.975,3)*b_sds/sqrt(3)

DF <- data.frame(treat=factor(rep(1:2, each=length(t))), 
                 time=rep(t, 2),
                 mean=c(a_means,b_means),
                 lwr=c(a_lwr,b_lwr),
                 upr=c(a_upr,b_upr))

library(ggplot2)
p <- ggplot(DF, aes(x=time, y=mean, ymin=lwr, ymax=upr)) +
  geom_ribbon(aes(fill=treat), alpha=0.3) +
  geom_line(aes(color=treat))
print(p)

means and 95 % confidence intervals

The way I'm constructing the confidence intervals obviously doesn't consider auto-correlation.

  • Is there a way to construct some kind of "auto-correlated confidence interval"?
  • Can I use the "un-correlated confidence interval"? Can I somehow estimate if it is too narrow or too wide in comparison to the auto-correlated confidence interval?
  • Is there a better approach to my problem?
$\endgroup$
1
  • $\begingroup$ Fascinating question, and it seems like this is a problem still being debated. Here's an interesting paper on the subject mpra.ub.uni-muenchen.de/31840 $\endgroup$
    – Drew75
    Oct 28, 2013 at 19:11

1 Answer 1

4
+100
$\begingroup$

Here are a couple thoughts that may be helpful:

  • Auto-correlation doesn't matter when you only look at a single t at a time. So, at a fixed time t, you could just run a t-test to check for a difference in means. If you run the t-test for each time separately, then you get a bunch of p-values. Because of auto-correlation these p-values are not independent, but each p-value considered alone is just fine.

  • So now you want to find the times for which there is a difference in means. I would try using false discovery rate (FDR) methods (see the "Benjamini-Hochberg procedure" at http://en.wikipedia.org/wiki/False_discovery_rate). Luckily, this procedure controls the FDR even when there is positive dependence among your p-values. (see "The Control of the False Discovery Rate in Multiple Testing under Dependency", free version here http://thom.jouve.free.fr/work/thesis/sitecopy_save/Biblio/ToCheck/fdr/Benjamini2001.pdf) This should give you a reasonable first answer to your original question.

  • Finally, I think the two plots you drew are very clear. They are probably more informative than any kind of statistical analysis you can run... Good luck!

Edit by Roland:

Here is an R implementation of the FDR method for the example in the question. The result looks reasonable.

dat <- setNames(cbind(stack(as.data.frame(t(a))), 
                      stack(as.data.frame(t(b)))), 
                c("a", "i", "b", "i"))
dat <- dat[,-4]
library(plyr)
p.raw <- ddply(dat, .(i), function(df) t.test(df$a, df$b)$p.value)
p.fdr <- cbind(p.adjust(p.raw[,2], method="fdr"),
               t[as.numeric(gsub("V","",p.raw[,1]))])
p.fdr[order(p.fdr[,2]),]

#             [,1] [,2]
#  [1,] 0.63001435    3
#  [2,] 0.19439226    4
#  [3,] 0.06200315    5
#  [4,] 0.07335654    6
#  [5,] 0.05336699    7
#  [6,] 0.06115999    8
#  [7,] 0.06115999    9
#  [8,] 0.06103370   10
#  [9,] 0.04324050   11
# [10,] 0.04324050   12
# [11,] 0.04324050   13
# [12,] 0.04324050   14
# [13,] 0.06103370   15
# [14,] 0.05533972   16
# [15,] 0.15489402   17
# [16,] 0.58234624   18
# [17,] 0.05533972   19
# [18,] 0.04324050   20
$\endgroup$
1
  • $\begingroup$ Thank you for pointing me to the FDR method. This looks promising. I took the liberty to add an implementation to your answer. $\endgroup$
    – Roland
    Nov 1, 2013 at 11:14

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.