Hi,

I got an experimental design made of 3 replicates. Each replicate consists in a time course comparing the response of a cell line to a treatment "A" vs a placebo "B". Replicates are not carried out on the same day, but within each replicate, treatment and placebo time course are carried out in parallel, usually starting from the cell line culture.

Experimental questions are:

• How do the two treatments differ in their response over time?
• Do the two treatments differ overall?

I read the section "9.7 Multi-level Experiments" of limma user guide and Aaron's answer in LIMMA time series spline model.matrix design setup. I think there are two blocks (one for subjects, one for replicates), but I am not sure how to specify them. In the code below, I give my view of the analysis with only blocking for subjects.

How to setup the design and the analysis?

Any help is appreciated.

    > cbind.data.frame(repl, trtm, timp, subj, g)[order(repl, trtm),]
       repl trtm timp subj  g
    1    R1    A    1  AR1 A1
    7    R1    A    2  AR1 A2
    13   R1    A    3  AR1 A3
    19   R1    A    4  AR1 A4
    25   R1    A    5  AR1 A5
    31   R1    A    6  AR1 A6
    2    R1    B    1  BR1 B1
    8    R1    B    2  BR1 B2
    14   R1    B    3  BR1 B3
    20   R1    B    4  BR1 B4
    26   R1    B    5  BR1 B5
    32   R1    B    6  BR1 B6
    3    R2    A    1  AR2 A1
    9    R2    A    2  AR2 A2
    15   R2    A    3  AR2 A3
    21   R2    A    4  AR2 A4
    27   R2    A    5  AR2 A5
    33   R2    A    6  AR2 A6
    4    R2    B    1  BR2 B1
    10   R2    B    2  BR2 B2
    16   R2    B    3  BR2 B3
    22   R2    B    4  BR2 B4
    28   R2    B    5  BR2 B5
    34   R2    B    6  BR2 B6
    5    R3    A    1  AR3 A1
    11   R3    A    2  AR3 A2
    17   R3    A    3  AR3 A3
    23   R3    A    4  AR3 A4
    29   R3    A    5  AR3 A5
    35   R3    A    6  AR3 A6
    6    R3    B    1  BR3 B1
    12   R3    B    2  BR3 B2
    18   R3    B    3  BR3 B3
    24   R3    B    4  BR3 B4
    30   R3    B    5  BR3 B5
    36   R3    B    6  BR3 B6

library(limma)

# design

trtm <- rep(c("A", "B"), 18)
timp <- rep(1:6, each=6)
repl <- rep(rep(c("R1", "R2", "R3"), each = 2), 6)
subj <- paste0(trtm, repl)

# g features points of interest to build contrasts
g <- factor(paste0(trtm, timp))

cbind.data.frame(trtm, timp, repl, subj, g)

# simulation + 1 experimental profile
med_sd = 0.04  # median sd observed
prfl = c(  # observed means
  A1=0.00, B1=0.20, A2=0.04, B2=0.60, A3=0.25, B3=0.95, 
  A4=0.10, B4=0.80, A5=0.15, B5=0.75, A6=0.15, B6=0.70)
set.seed(0)
meas = matrix(rnorm(250*length(g), sd = med_sd), ncol = length(g))
for (i in names(prfl)) {
  meas[1, g==i] = meas[1, g==i] + prfl[i]
}
set.seed(1)
meas[-1,] = meas[-1,] + runif(nrow(meas)-1, 0, .3)
rownames(meas) = sprintf("G%03d", 1:nrow(meas))
dim(meas)

# graphical check
library(ggplot2)

ggplot(data.frame(trtm, timp, repl, subj, g, value=meas[1,]), 
       aes(timp, value, group=subj, col=repl)) + geom_point() + geom_line()

# limma design
design <- model.matrix(~ 0 + g)
colnames(design) <- levels(g)

# with blocking on subject in duplicateCorrelation().
corfit <- duplicateCorrelation(meas, design, block=subj)
corfit$consensus

# Fit the model 
fit <- lmFit(meas, design, block=subj, correlation=corfit$consensus)


# How do the two treatments differ in their response over time?

con.1 <- makeContrasts(B1 - A1, levels=design)
fitcon.1 <- contrasts.fit(fit, con.1)
fitcon.1 <- eBayes(fitcon.1) 
topTable(fitcon.1, number=5, adjust.method="BH")

# and so on...
# con.2 <- makeContrasts(B2 - A2, levels=design)
# con.3 <- makeContrasts(B3 - A3, levels=design)
# con.4 <- makeContrasts(B4 - A4, levels=design)
# con.5 <- makeContrasts(B5 - A5, levels=design)
# con.6 <- makeContrasts(B6 - A6, levels=design)

# Do the two treatments differ overall?

anova.con <- makeContrasts(
  (B1 - A1), (B2 - A2), (B3 - A3), 
  (B4 - A4), (B5 - A5), (B6 - A6),
  levels=design)
fitcon.anova <- contrasts.fit(fit, anova.con)
fitcon.anova <- eBayes(fitcon.anova) 
topTable(fitcon.anova, number=5, adjust.method="BH")
# look at F test and p-value

# here we identify which tests to keep
rescon.anova = topTable(fitcon.anova, number=Inf, adjust.method="BH")
decideTests(fitcon.anova, method="nestedF")

The treatments are paired within replicate. If you want to adjust for both subject and replicate, you can include replicate in the design matrix and subject as a duplicateCorrelation block.