Dear Dr. Smyth,
As other limma users I got a 'no residual degrees of freedom' when
trying to calculate the Bayes statistics:
fit <- lmFit(MA, design)
contrast.matrix <- makeContrasts(R10, R40, N2L, R10-R40, R40-N2L,
R10-N2L, levels=design)
cfit <- contrasts.fit(fit, contrast.matrix)
eb <- eBayes(cfit)
As you indicated in an earliner message, the reason for this is that
all fit$df.residual are missing.
There is this code segment in the lm.series function which is called
by lmFit():
if(!is.null(weights)) {
weights <- as.matrix(weights)
if(any(dim(weights) != dim(M))) weights <-
array(weights,dim(M))
weights[weights <= 0] <- NA
M[!is.finite(weights)] <- NA
}
This replaces elements of the weight matrix which are 0 or less and
also the corresponding elements in MA$M with NA. My weight matrix
consists of zeros and negative numbers which were generated when
reading in the Stanford Microarray Database (SMD) files with the
command:
RG<-read.maimages(targets$FileName, source="smd", fill=T,
wt.fun=function(x) {return(x$FLAG)})
This copies the FLAG column values into the weight matrix and in the
raw smd file this column consists of zeros and negative numbers.
Modifying the wt.fun parameter of read.maimages and using
wt.fun=wtflags(w=0.1)
didn't work and generated an error message since there is no column
named 'Flags', the SMD file has a 'FLAG' column. A simple change in
the wtflags function clears all this error messages.
wtflags <- function(w=0.1)
# Quality weights based on Flags from GenePix output
# Gordon Smyth
# 9 March 2003. Last revised 11 June 2003.
function(gpr) {
# changed "Flags" to "FLAG"
flagged <- (gpr[,"FLAG"] < 0)
w*flagged + !flagged
}
Kind regards,
Juerg Straubhaar, PhD
Umass Med
Dear Juerg,
'Weights' refer to the weight that each observation gets in the
residual
analyses. Weight 1 means a full observation, 0.5 means half weights,
weight 2 is equivalent to having two observations at that value, zero
means that observation contributes nothing. Weights below zero have
no
meaning. Weights have the same meaning in limma as in any other
regression program in R such as lm, glm etc etc.
You are creating weights which are all zero or negative and you're
surprised that eBayes doesn't give you any results?
If you look at the bottom of your email you'll see that wtflags() is
only
for GenePix data. This means that you can't use it for SMD data.
You need to consider what values the SMD FLAGS column takes and you
need
to decide what actual weight you want to give to each FLAG value.
There
is no shortcut. Then you can create a meaningful weight function for
use
with read.maimages.
Gordon
> Dear Dr. Smyth,
>
> As other limma users I got a 'no residual degrees of freedom' when
> trying to calculate the Bayes statistics:
>
> fit <- lmFit(MA, design)
> contrast.matrix <- makeContrasts(R10, R40, N2L, R10-R40, R40-N2L,
> R10-N2L, levels=design) cfit <- contrasts.fit(fit, contrast.matrix)
> eb <- eBayes(cfit)
>
> As you indicated in an earliner message, the reason for this is that
all
> fit$df.residual are missing.
>
> There is this code segment in the lm.series function which is called
by
> lmFit():
>
> if(!is.null(weights)) {
> weights <- as.matrix(weights)
> if(any(dim(weights) != dim(M))) weights <-
array(weights,dim(M))
> weights[weights <= 0] <- NA
> M[!is.finite(weights)] <- NA
> }
>
> This replaces elements of the weight matrix which are 0 or less and
also
> the corresponding elements in MA$M with NA. My weight matrix
consists of
> zeros and negative numbers which were generated when reading in the
> Stanford Microarray Database (SMD) files with the command:
>
> RG<-read.maimages(targets$FileName, source="smd", fill=T,
> wt.fun=function(x) {return(x$FLAG)})
>
> This copies the FLAG column values into the weight matrix and in the
raw
> smd file this column consists of zeros and negative numbers.
>
> Modifying the wt.fun parameter of read.maimages and using
>
> wt.fun=wtflags(w=0.1)
>
> didn't work and generated an error message since there is no column
> named 'Flags', the SMD file has a 'FLAG' column. A simple change in
the
> wtflags function clears all this error messages.
>
> wtflags <- function(w=0.1)
> # Quality weights based on Flags from GenePix output
> # Gordon Smyth
> # 9 March 2003. Last revised 11 June 2003.
>
> function(gpr) {
> # changed "Flags" to "FLAG"
> flagged <- (gpr[,"FLAG"] < 0)
> w*flagged + !flagged
> }
>
> Kind regards,
> Juerg Straubhaar, PhD
> Umass Med
