At 09:15 PM 26/05/2007, Lev Soinov wrote:
>Dear List,
>It was a lot of discussion recently about paired designs.
>I suppose my question is very simple:
>I have 6 paired arrays (3 pairs, treatment vs control within each
>pair, 6 arrays in total).
>Could somebody explain briefly if this is correct that two scripts
>below (with different designs) should produce the same results for
>paired design in LIMMA(i am getting the same results)?
These are simply two different parametrizations for the same
experiment, hence they give the same result.
Is this so surprising? There are many examples in the limma User's
Guide which explain the equivalence of these sort of design matrices.
> And also, am I creating the design and targets objects correctly?
Is there any reason why you don't follow the section of the limma
User's Guide which deals with this design, Section 8.3 "Paired
Samples"?
You can experiment with other ways to define the design matrix if you
want. Surely you can check for yourself whether your approach is
correct by checking whether it gives the same answer as that from the
User's Guide.
Best wishes
Gordon
>1st script:
> >data<-ReadAffy()
> >temp<-rma(data)
> >targets <- readTargets("targets.txt")
> > targets
> FileName Pair Treatment
>1 1.CEL 1 C
>2 2.CEL 1 T
>3 3.CEL 2 C
>4 4.CEL 2 T
>5 5.CEL 3 C
>6 6.CEL 3 T
> >Pair <- factor(targets$Pair)
> >Treat <- factor(targets$Treatment, levels=c("C","T"))
> >design <- model.matrix(~Pair+Treat)
> > design
> (Intercept) Pair2 Pair3 TreatT
>1 1 0 0 0
>2 1 0 0 1
>3 1 1 0 0
>4 1 1 0 1
>5 1 0 1 0
>6 1 0 1 1
> >fit_pair <- lmFit(temp, design)
> >fit_pair <- eBayes(fit_pair)
> >topTable(fit_pair, coef="TreatT")
>
>2nd script:
> >data<-ReadAffy()
> >temp<-rma(data)
> > design2<-cbind(c(1,1,0,0,0,0), c(0,0,1,1,0,0), c(0,0,0,0,1,1),
> c(0,1,0,1,0,1))
> > colnames(design2)<-c("Pair1", "Pair2", "Pair3", "TreatT")
> > design2
> Pair1 Pair2 Pair3 TreatT
>[1,] 1 0 0 0
>[2,] 1 0 0 1
>[3,] 0 1 0 0
>[4,] 0 1 0 1
>[5,] 0 0 1 0
>[6,] 0 0 1 1
> > fit <- lmFit(temp, design2)
> > fit <- eBayes(fit)
> > topTable(fit, coef=4)
>
>With kind regards,
>Lev.
