Entering edit mode
Dear Naomi,
duplicateCorrelation() calls mixedModel2Fit() (in the statmod
package),
which in turn calls glmgam.fit(). glmgam.fit() implements a
Levenberg-Marquardt modification to the Fisher scoring algorithm to
prevent divergence, so it is extremely reliable, much more so than a
call
to glm() would be. On the other hand, intraspotCorrelation() calls
remlscore() (also in the statmod package), which does ordinary Fisher
scoring for a related model without any modification to prevent
divergence. So it is entirely possible that duplicateCorrelation()
will
work for a dataset for which intraspotCorrelation() does not.
However, to get a NULL result from intraspotCorrelation(), the model
fit
would have to fail for every single gene in your dataset. Just
failing
for one or two would not cause a problem.
Best wishes
Gordon
> Date: Thu, 24 Feb 2011 20:20:12 -0500
> From: Naomi Altman <naomi at="" stat.psu.edu="">
> To: bioconductor at r-project.org
> Subject: [BioC] intraspotCorrelation vs duplicateCorrelation
> Message-ID: <6.2.3.4.1.20110224200727.01ffea98 at imap.stat.psu.edu>
> Content-Type: text/plain; charset="us-ascii"; format=flowed
>
> I was having problems with intraspotCorrelation for a single channel
> analysis, so I decided to reorganize the data to use
> duplicateCorrelation. I thought I would get the same answer, but I
> did not. I hope someone can tell me why. (As you will see below, I
> like the answer with duplicateCorrelation better!)
>
> MAp is the normalized 2-channel data.
> I have a targets list called targetp with the treatment
> information. The code is below.
>
> ##############################################################
> # separate channel approach
> #############################################################
> targetSingle <- targetsA2C(targetp)
> rep=c(0,0,1,1,0,0,1,1,0,0,1,1)
> designp=model.matrix(~-1+targetSingle$Target+rep,levels=unique(c(cy3
,cy5)))
> corfit =intraspotCorrelation(MAp, designp)
> corfit$cor
> NULL
> ##########################################################
> # the warning was
> #In remlscore(y, X, Z) : reml: Max iterations exceeded
> ##########################################################
>
> ############################################################
> # approach treating channels as separate arrays
> ############################################################
> RGp=RG.MA(MAp)
> Rp=log2(RGp$R)
> Gp=log2(RGp$G)
> exprP=cbind(Rp[,1],Gp[,1],Rp[,2],Gp[,2],Rp[,3],Gp[,3],Rp[,4],Gp[,4],
Rp[,5],Gp[,5],Rp[,6],Gp[,6])
> corfitP = duplicateCorrelation(exprP, designp, block =
> factor(c(1,1,2,2,3,3,4,4,5,5,6,6)))
> corfitP$cor
> [1] 0.4921254
>
> sessionInfo()
> R version 2.11.1 (2010-05-31)
> i386-pc-mingw32
>
> locale:
> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
> States.1252 LC_MONETARY=English_United States.1252
> [4] LC_NUMERIC=C LC_TIME=English_United
> States.1252
>
> attached base packages:
> [1] tools stats graphics grDevices
> utils datasets methods base
>
> other attached packages:
> [1] statmod_1.4.8 rat2302cdf_2.6.0
> limma_3.4.5 affy_1.26.1 Biobase_2.8.0
>
> loaded via a namespace (and not attached):
> [1] affyio_1.16.0 preprocessCore_1.10.0
>
>
> Thanks,
> Naomi
______________________________________________________________________
The information in this email is confidential and
intend...{{dropped:4}}
