Entering edit mode
*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
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
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...
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