Problem obtaining standard error from limma
1
0
Entering edit mode
@jose-luis-lavin-5531
Last seen 6.1 years ago
Spain

*Question posted as new question thread as requested by Aaron Lun*

Dear limma team,

First of all, I want to thank you for this wonderful package. This said, I need to ask you a question on the subject mentioned in this post, Standard error calculation from my fit object using the command in the post:

>std.error =  fit$stdev.unscaled * fit$sigma
> std.error
numeric(0)
> fit$stdev.unscaled
NULL
> fit$sigma
NULL

I'd ask if some experienced limma user may help me figuring out what the problem is.

Here is the structure of my fit object:

> str(fit)

List of 1
 $ :Formal class 'MArrayLM' [package "limma"] with 1 slot
  .. ..@ .Data:List of 12
  .. .. ..$ : num [1:19568, 1:2] 7.93 9.38 7.46 7.81 10.88 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:19568] "ILMN_2417611" "ILMN_2896528" "ILMN_2721178" "ILMN_3033922" ...
  .. .. .. .. ..$ : chr [1:2] "g1" "g2"
  .. .. ..$ : int 2
  .. .. ..$ : NULL
  .. .. ..$ :List of 5
  .. .. .. ..$ qr   : num [1:24, 1:2] -3.464 0.289 0.289 0.289 0.289 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. ..$ : chr [1:2] "g1" "g2"
  .. .. .. ..$ qraux: num [1:2] 1.29 1
  .. .. .. ..$ pivot: int [1:2] 1 2
  .. .. .. ..$ tol  : num 1e-07
  .. .. .. ..$ rank : int 2
  .. .. .. ..- attr(*, "class")= chr "qr"
  .. .. ..$ : int [1:19568] 22 22 22 22 22 22 22 22 22 22 ...
  .. .. ..$ : Named num [1:19568] 0.213 0.156 0.16 0.26 0.135 ...
  .. .. .. ..- attr(*, "names")= chr [1:19568] "ILMN_2417611" "ILMN_2896528" "ILMN_2721178" "ILMN_3033922" ...
  .. .. ..$ : num [1:2, 1:2] 0.0833 0 0 0.0833
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "g1" "g2"
  .. .. .. .. ..$ : chr [1:2] "g1" "g2"
  .. .. ..$ : num [1:19568, 1:2] 0.289 0.289 0.289 0.289 0.289 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:19568] "ILMN_2417611" "ILMN_2896528" "ILMN_2721178" "ILMN_3033922" ...
  .. .. .. .. ..$ : chr [1:2] "g1" "g2"
  .. .. ..$ : int [1:2] 1 2
  .. .. ..$ : Named num [1:19568] 7.85 9.39 7.46 7.89 10.86 ...
  .. .. .. ..- attr(*, "names")= chr [1:19568] "ILMN_2417611" "ILMN_2896528" "ILMN_2721178" "ILMN_3033922" ...
  .. .. ..$ : chr "ls"
  .. .. ..$ : num [1:24, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:2] "g1" "g2"

 

Please accept my apologies if this can be considered a very naive question...

Best wishes

JL

 

limma standard error • 2.2k views
ADD COMMENT
0
Entering edit mode

Here is my session info (I couldn't include it in the previous post...

And here is my session info:

Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252    LC_MONETARY=Spanish_Spain.1252
[4] LC_NUMERIC=C                   LC_TIME=Spanish_Spain.1252    

attached base packages:
 [1] grid      stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lumiMouseIDMapping_1.10.0 lumiMouseAll.db_1.22.0    org.Mm.eg.db_3.2.3       
 [4] fields_8.3-5              maps_3.0.0-2              spam_1.3-0               
 [7] igraph_1.0.1              gtools_3.5.0              gplots_2.17.0            
[10] KEGG.db_3.2.2             KEGGgraph_1.28.0          GOstats_2.36.0           
[13] Category_2.36.0           Matrix_1.2-2              lumiHumanIDMapping_1.10.0
[16] annotate_1.48.0           XML_3.98-1.3              lumiHumanAll.db_1.22.0   
[19] limma_3.26.3              lumi_2.22.0               Rgraphviz_2.14.0         
[22] graph_1.48.0              SparseM_1.7               clusterProfiler_2.99.1   
[25] DOSE_2.9.7                GO.db_3.2.2               org.Hs.eg.db_3.2.3       
[28] RSQLite_1.0.0             DBI_0.3.1                 AnnotationDbi_1.32.0     
[31] IRanges_2.4.4             S4Vectors_0.8.3           Biobase_2.30.0           
[34] AnnotationHub_2.2.2       BiocGenerics_0.16.1 

ADD REPLY
0
Entering edit mode

Dear Gordon,

Thank you for your swift answer.

I recreated the object and got the following fit:

> show(fit)
[[1]]
An object of class "MArrayLM"
$coefficients
                    g1        g2
ILMN_2417611  7.933408  7.760932
ILMN_2896528  9.382122  9.390301
ILMN_2721178  7.462175  7.456717
ILMN_3033922  7.812959  7.967204
ILMN_3092673 10.877547 10.836764
19563 more rows ...

$rank
[1] 2

$assign
NULL

$qr
$qr
             g1        g2
[1,] -3.4641016  0.000000
[2,]  0.2886751 -3.464102
[3,]  0.2886751  0.000000
[4,]  0.2886751  0.000000
[5,]  0.2886751  0.000000
19 more rows ...

$qraux
[1] 1.288675 1.000000

$pivot
[1] 1 2

$tol
[1] 1e-07

$rank
[1] 2

$df.residual
[1] 22 22 22 22 22
19563 more elements ...

$sigma
ILMN_2417611 ILMN_2896528 ILMN_2721178 ILMN_3033922 ILMN_3092673
   0.2125742    0.1558975    0.1598369    0.2595060    0.1346700
19563 more elements ...

$cov.coefficients
           g1         g2
g1 0.08333333 0.00000000
g2 0.00000000 0.08333333

$stdev.unscaled
                    g1        g2
ILMN_2417611 0.2886751 0.2886751
ILMN_2896528 0.2886751 0.2886751
ILMN_2721178 0.2886751 0.2886751
ILMN_3033922 0.2886751 0.2886751
ILMN_3092673 0.2886751 0.2886751
19563 more rows ...

$pivot
[1] 1 2

$Amean
ILMN_2417611 ILMN_2896528 ILMN_2721178 ILMN_3033922 ILMN_3092673
    7.847170     9.386212     7.459446     7.890081    10.857155
19563 more elements ...

$method
[1] "ls"

$design
     g1 g2
[1,]  1  0
[2,]  1  0
[3,]  1  0
[4,]  1  0
[5,]  1  0
19 more rows ...

I can see fit$stdev.unscaled and fit$sigma in the object, but I have the same problem when I try to calculate the standard error...

>std.error =  fit$stdev.unscaled * fit$sigma
> std.error
numeric(0)

I use a script created by a former member of my lab where the limma part of the code, where fit is created is the following:

fit = ebfit = list()
for (i in 1:nrow(contrasts)) {
  nelem = c(0,0)
  colvals = c()
  for (cc in 1:2) {
    elements = strsplit(contrasts[i,cc],";")[[1]]
    for (e in elements) {
      if (length(grep("!", e))>0) {
        subelements = strsplit(e,"!")[[1]]
        colval = exprs(lumi.N)[,samples.pass[,1]==subelements[1],drop=F] -
                 exprs(lumi.N)[,samples.pass[,1]==subelements[2],drop=F]
      } else {
        colval = exprs(lumi.N)[,samples.pass[,1]==e,drop=F]
      }
      colvals = cbind(colvals, colval)
      nelem[cc] = nelem[cc] + ncol(colval)
    }
  }
  designcol = c(rep(1,nelem[1]),rep(0,nelem[2]))
  design = cbind(g1=designcol,g2=as.integer(!designcol))
  fit[[i]] = lmFit(colvals, design)
  cont.matrix = makeContrasts(contrasts=paste(colnames(design),collapse="-"), levels=design)
  ebfit[[i]] = eBayes(contrasts.fit(fit[[i]], cont.matrix))
}

I also tried to calculate the standard error using the ebfit object, but it returns a NULL result as well.

I hope this helps to clear out this issue.

Best wishes and thanks in advance.

JL

ADD REPLY
2
Entering edit mode
@gordon-smyth
Last seen 3 hours ago
WEHI, Melbourne, Australia

Your fit object is not valid. I suggest that you recreate the object.

If the problem persists, then show us the code you used to create 'fit'.

By the way, to see the structure of a limma fit object, best to use show(fit) rather than str(fit).

14 hours later:

This is a question about programming in R rather than an issue with limma. The cause of the confusion is that ebfit is a list of fit objects rather than a fit object itself. Indeed you created it by ebfit <- list(). If you type

class(ebfit)

you will find it is a simple list. On the other hand, if you type

class(ebfit[[1]])
names(ebfit[[1]])

you will find that ebfit[[1]] is a fit object. Note that '[[1]]' selects the first component of the list. Double brackets are needed for list components. List objects are ubiquitous in R, so you need to become very familiar with them.

To compute standard errors you need

std.error <- ebfit[[1]]$stdev.unscaled * ebfit[[1]]$sigma

 

If I may say, the script written by your former colleague is unnecessarily complicated for the purpose you are using it. You are only creating a single fit object, so nesting it within another list object is unnecessarily obscure. Usually it would be better for you to write your own simpler code directly so that you understand all the steps. If you do have to use your former colleague's script, then it would be best to send questions directly to him or her, because the script has a lot of unusual aspects.

ADD COMMENT
0
Entering edit mode

Dear Gordon,

I totally agree with you.

This time, I asked in this list because I needed a clear explanation and a solution for the nth issue with this code...Since I do not have to use it frequently, I never seem to find the time to invest into reprograming this pipeline.

Please accept my apologies if the question was not appropriate, but I was completely overwhelmed by the code and its associated errors and missenses.

Best wishes and thank you again for your kind help.

JL

ADD REPLY

Login before adding your answer.

Traffic: 821 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6