Hi all. 
 
I have tried to read previous posts on this matter but I found the
answers kind of vague (perhaps I missed something due to my limited
knowledge in statistics, in that case please guide me to an explaining
post and I apologize for asking the same question again).
 
I have affydata on 4 patients (paired) before and after treatment
(unbalance due to different number of biopsies from each patient). I
want to fit the patients as a fixed effect. I am interested in the
effect due to treatment.
 
In previous posts I have seen Gordon give directions about two ways to
fit patients as a fixed effect with LIMMA.
 
http://files.protsuggest.org/biocond/html/8297.html
 
and
 
http://files.protsuggest.org/biocond/html/8175.html
 
In the first post the t-statistic is according to Gordon suppose to be
the corresponding paired t-statistic.
But in the second mail he states “All forms of blocking are much the
same from a mathematical point of view, whether they are technical
replicates or blocking on patient. 
The direct extension of a paired t-test would actually arise from
fitting patient as a fixed effect. If you did that, limma would compute
for you paired t-tests with moderated denominators. If you use the
'block' argument instead of using a fixed effect, I would generally
recommend that you estimate the common correction rather than take a
preset value. However I suspect that, in this case, if you use a large
correlation like 0.75 or higher the results might be very similar to
using a fixed patient effect.”
 
My situation look like this: 
#########################################################
#Fitting the patients as a fixed effect, I am not totally shure of my
code, please correct me if I am #wrong
eset<-rma(Data)
tmp<-pData(eset) 
tmp<-cbind(tmp,Before=c(1,0,1,1,1,0,1,0,1,0,1,0,1,0))
tmp<-cbind(tmp,After=c(0,1,0,0,0,1,0,1,0,1,0,1,0,1))
tmp<-cbind(tmp,Patient=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
pData(eset)<-tmp 
design <- model.matrix(~(After - Before) + Patient,
data=as.data.frame(tmp)) 
fit<-lmFit(eset,design) 
fitModel <- eBayes(fit) 
 
#########################################################
#Using the block argument to “block” patient.
eset<-rma(Data)
tmp<-pData(eset) 
tmp<-cbind(tmp,Before=c(1,0,1,1,1,0,1,0,1,0,1,0,1,0))
tmp<-cbind(tmp,After=c(0,1,0,0,0,1,0,1,0,1,0,1,0,1))
pData(eset) <- tmp
design <- model.matrix(~(After - Before), data=pData(eset)) 
cor<-duplicateCorrelation(eset,
design,block=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
cor$cor
# [1] 0.2739284
fit<-lmFit(eset, design, block=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
fitcor <- eBayes(fit) 
 
# I also try this model using the default value, cor=0.75, when using
lmFit.
eset<-rma(Data)
tmp<-pData(eset) 
tmp<-cbind(tmp,Before=c(1,0,1,1,1,0,1,0,1,0,1,0,1,0))
tmp<-cbind(tmp,After=c(0,1,0,0,0,1,0,1,0,1,0,1,0,1))
pData(eset) <- tmp
design <- model.matrix(~(After - Before), data=pData(eset)) 
fit<-lmFit(eset, design, block=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
fitCorDefault <- eBayes(fit) 
fitCorDefault$cor
[1] 0.75
#########################################################
 
First of all, is my code right? I am kind of sceptic since Gordon states
that “if you use a large correlation like 0.75 or higher the results
might be very similar to using a fixed patient effect”. When I plot the
B-scores of both models they differ somewhat in the ranking and in the
number of genes with B-score >0. They look kind of the same but the
B-scores is higher when I use the block argument to “block” patient
using the default correlation. 
 
colnames(fitModel$lods)
#[1] "(Intercept)" "After"       "Patient" 
length(which(fitModel$lods[,2] > 0))
#[1] 0
 
fitcor$cor
#[1] 0.2739284
colnames(fitcor$lods)
#[1] "(Intercept)" "After"
length(which(fitcor$lods[,2] > 0))
#[1] 4
 
fitCorDefault$cor
#[1] 0.75
colnames(fitCorDefault$lods)
#[1] "(Intercept)" "After"
length(which(fitCorDefault$lods[,2] > 0))
#[1] 109
 
#################
#To get a feeling of similarity between the ranking of the genes
cor(fitcor$lods[,2], fitCorDefault$lods[,2],
use="pairwise.complete.obs")
[1] 0.9781474
cor(fitcor$lods[,2], fitModel$lods[,2], use="pairwise.complete.obs")
[1] 0.9753385
cor(fitCorDefault$lods[,2], fitModel$lods[,2],
use="pairwise.complete.obs")
[1] 0.9610673 
################
 
So my questions are, if my code is right, which way is the right way to
fit patients as a fixed effect since the result differ?
 
If my code is wrong, please guide me of how to get it right. 
 
Thank you for your time
 
Best regards
 
Johan Lindberg
 
 
************************************************************************
*******************
Johan Lindberg
Royal Institute of Technology
AlbaNova University Center
Stockholm Center for Physics, Astronomy and Biotechnology
Department of Molecular Biotechnology
106 91 Stockholm, Sweden
 
Phone (office) +46 8 553 783 44
Fax + 46 8 553 784 81
Visiting adress Roslagstullsbacken 21, Floor 3
Delivery adress Roslagsvägen 30B
http://www.biotech.kth.se/molbio/microarray/index.html
************************************************************************
******************* 
 
 

	[[alternative HTML version deleted]]

