limma: NA coefficients due to missing values
1
0
Entering edit mode
Shi, Tao ▴ 720
@shi-tao-199
Last seen 7.1 years ago
Dear Gordon and list, Please see the example below and R output.? The coefficients for Probe #2 can't be estimated b/c patient 10 data were missing.? But in fact, they can be estimated just using data from other 9 patients.? Is there a way to tell limma to do that? Thanks! Tao library(limma) dat <- matrix(rnorm(400), ncol=40) ptID <- factor(rep(1:10, each=4)) day <- factor(rep(1:4,10)) dat[2, 37:40] <- NA design <- model.matrix(~0+ptID+day) contrast.matrix <- makeContrasts(day2=day2, day3=day3,day4=day4, levels=design) fit1 <- lmFit(dat, design=design) fit2 <- eBayes(contrasts.fit(fit1, contrast.matrix)) Warning message: Partial NA coefficients for 1 probe(s) > fit2$p.value ?????? Contrasts ????????????? day2????? day3????? day4 ?? [1,] 0.40027892 0.4917185 0.0387787 ?? [2,]???????? NA??????? NA??????? NA ?? [3,] 0.05281279 0.4914354 0.6088744 ?? [4,] 0.72648798 0.6816108 0.7551084 ?? [5,] 0.35528206 0.7962423 0.9827632 ?? [6,] 0.09552328 0.3458950 0.1327465 ?? [7,] 0.76918217 0.5945065 0.4511735 ?? [8,] 0.40125232 0.1066155 0.1269416 ?? [9,] 0.56923379 0.2981261 0.9285374 ? [10,] 0.04276991 0.4693320 0.3139434 > fit2$coef ?????? Contrasts ????????????? day2?????? day3???????? day4 ?? [1,]? 0.3737307? 0.3054398? 0.921284513 ?? [2,]???????? NA???????? NA?????????? NA ?? [3,] -0.8628711 -0.3056396 -0.227256710 ?? [4,] -0.1553395? 0.1821986? 0.138509886 ?? [5,] -0.4107858? 0.1146613 -0.009593154 ?? [6,]? 0.7421089 -0.4188821? 0.668951727 ?? [7,] -0.1303084 -0.2364268 -0.334736056 ?? [8,] -0.3729581 -0.7182350 -0.679192601 ?? [9,]? 0.2528092? 0.4624636? 0.039823091 ? [10,] -0.9030535 -0.3214420 -0.447554386 > sessionInfo() R version 3.0.0 (2013-04-03) Platform: x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252? LC_CTYPE=English_United States.1252??? LC_MONETARY=English_United States.1252 LC_NUMERIC=C????????????????????????? [5] LC_TIME=English_United States.1252??? attached base packages: [1] grDevices datasets? splines?? graphics? stats???? tcltk???? utils???? methods?? base???? other attached packages: [1] limma_3.16.3??? svSocket_0.9-55 TinnR_1.0-5???? R2HTML_2.2.1??? Hmisc_3.10-1.1? survival_2.37-4 loaded via a namespace (and not attached): [1] cluster_1.14.4? grid_3.0.0????? lattice_0.20-15 svMisc_0.9-69?? tools_3.0.0???
limma limma • 5.0k views
0
Entering edit mode
@gordon-smyth
Last seen just now
WEHI, Melbourne, Australia
Dear Tao, The coefficients for days 2-4 have already been estimated without NAs, see fit1$coef[,11:13]. Using contrasts will introduce NAs unless you remove the missing patient effect. This can be done by, > fit.day <- fit1[,11:13] > colnames(fit.day) [1] "day2" "day3" "day4" > head(fit.day$coef) day2 day3 day4 [1,] 0.1018951 0.14728872 0.01348279 [2,] 0.8222967 0.95183289 0.02500825 [3,] 0.6601298 -0.27436185 0.16624044 [4,] -0.7082047 0.04495111 0.25636112 [5,] -0.0628469 -0.04033322 -0.59622820 [6,] -0.0296513 -0.31283864 0.14225534 Then you can happily take as many contrasts as you like of the day effects, using fit3 instead of fit1. For example: cont.mat <- makeContrasts(day3-day2, levels=colnames(fit.day)) Best wishes Gordon --------------------------------------------- Professor Gordon K Smyth, Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, Vic 3052, Australia. http://www.statsci.org/smyth On Tue, 14 May 2013, Shi, Tao wrote: > Dear Gordon and list, Please see the example below and R output.? The coefficients for Probe #2 can't be estimated b/c patient 10 data were missing.? But in fact, they can be estimated just using data from other 9 patients.? Is there a way to tell limma to do that? Thanks! Tao library(limma) dat <- matrix(rnorm(400), ncol=40) ptID <- factor(rep(1:10, each=4)) day <- factor(rep(1:4,10)) dat[2, 37:40] <- NA design <- model.matrix(~0+ptID+day) contrast.matrix <- makeContrasts(day2=day2, day3=day3,day4=day4, levels=design) fit1 <- lmFit(dat, design=design) fit2 <- eBayes(contrasts.fit(fit1, contrast.matrix)) Warning message: Partial NA coefficients for 1 probe(s) > fit2$p.value ?????? Contrasts ????????????? day2????? day3????? day4 ?? [1,] 0.40027892 0.4917185 0.0387787 ?? [2,]???????? NA??????? NA??????? NA ?? [3,] 0.05281279 0.4914354 0.6088744 ?? [4,] 0.72648798 0.6816108 0.7551084 ?? [5,] 0.35528206 0.7962423 0.9827632 ?? [6,] 0.09552328 0.3458950 0.1327465 ?? [7,] 0.76918217 0.5945065 0.4511735 ?? [8,] 0.40125232 0.1066155 0.1269416 ?? [9,] 0.56923379 0.2981261 0.9285374 ? [10,] 0.04276991 0.4693320 0.3139434 > fit2$coef ?????? Contrasts ????????????? day2?????? day3???????? day4 ?? [1,]? 0.3737307? 0.3054398? 0.921284513 ?? [2,]???????? NA???????? NA?????????? NA ?? [3,] -0.8628711 -0.3056396 -0.227256710 ?? [4,] -0.1553395? 0.1821986? 0.138509886 ?? [5,] -0.4107858? 0.1146613 -0.009593154 ?? [6,]? 0.7421089 -0.4188821? 0.668951727 ?? [7,] -0.1303084 -0.2364268 -0.334736056 ?? [8,] -0.3729581 -0.7182350 -0.679192601 ?? [9,]? 0.2528092? 0.4624636? 0.039823091 ? [10,] -0.9030535 -0.3214420 -0.447554386 > sessionInfo() R version 3.0.0 (2013-04-03) Platform: x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252? LC_CTYPE=English_United States.1252??? LC_MONETARY=English_United States.1252 LC_NUMERIC=C????????????????????????? [5] LC_TIME=English_United States.1252??? attached base packages: [1] grDevices datasets? splines?? graphics? stats???? tcltk???? utils???? methods?? base???? other attached packages: [1] limma_3.16.3??? svSocket_0.9-55 TinnR_1.0-5???? R2HTML_2.2.1??? Hmisc_3.10-1.1? survival_2.37-4 loaded via a namespace (and not attached): [1] cluster_1.14.4? grid_3.0.0????? lattice_0.20-15 svMisc_0.9-69?? tools_3.0.0??? > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:5}}
0
Entering edit mode
Thank you very much, Gordon!? That will do.? I also just saw one of your older posts ( https://stat.ethz.ch/pipermail/bioconductor/2007-October/019889.html? ). ? It clears things up even better. Tao >________________________________ > From: Gordon K Smyth <smyth at="" wehi.edu.au=""> >To: "Shi, Tao" <shidaxia at="" yahoo.com=""> >Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> >Sent: Tuesday, May 14, 2013 4:34 PM >Subject: Re: limma: NA coefficients due to missing values > > >Dear Tao, > >The coefficients for days 2-4 have already been estimated without NAs, see >fit1$coef[,11:13]. > >Using contrasts will introduce NAs unless you remove the missing patient >effect.? This can be done by, > >? > fit.day <- fit1[,11:13] >? > colnames(fit.day) >? [1] "day2" "day3" "day4" >? > head(fit.day$coef) >? ? ? ? ? ? day2? ? ? ? day3? ? ? ? day4 >? [1,]? 0.1018951? 0.14728872? 0.01348279 >? [2,]? 0.8222967? 0.95183289? 0.02500825 >? [3,]? 0.6601298 -0.27436185? 0.16624044 >? [4,] -0.7082047? 0.04495111? 0.25636112 >? [5,] -0.0628469 -0.04033322 -0.59622820 >? [6,] -0.0296513 -0.31283864? 0.14225534 > >Then you can happily take as many contrasts as you like of the day >effects, using fit3 instead of fit1.? For example: > >???cont.mat <- makeContrasts(day3-day2, levels=colnames(fit.day)) > >Best wishes >Gordon > >--------------------------------------------- >Professor Gordon K Smyth, >Bioinformatics Division, >Walter and Eliza Hall Institute of Medical Research, >1G Royal Parade, Parkville, Vic 3052, Australia. >http://www.statsci.org/smyth > >On Tue, 14 May 2013, Shi, Tao wrote: > >> Dear Gordon and list, > >Please see the example below and R output.? The coefficients for Probe #2 >can't be estimated b/c patient 10 data were missing.? But in fact, they >can be estimated just using data from other 9 patients.? Is there a way to >tell limma to do that? > >Thanks! > >Tao > > > >library(limma) >dat <- matrix(rnorm(400), ncol=40) >ptID <- factor(rep(1:10, each=4)) >day <- factor(rep(1:4,10)) >dat[2, 37:40] <- NA >design <- model.matrix(~0+ptID+day) >contrast.matrix <- makeContrasts(day2=day2, day3=day3,day4=day4, levels=design) >fit1 <- lmFit(dat, design=design) >fit2 <- eBayes(contrasts.fit(fit1, contrast.matrix)) > > >Warning message: >Partial NA coefficients for 1 probe(s) >> fit2$p.value >?????? Contrasts >????????????? day2????? day3????? day4 >?? [1,] 0.40027892 0.4917185 0.0387787 >?? [2,]???????? NA??????? NA??????? NA >?? [3,] 0.05281279 0.4914354 0.6088744 >?? [4,] 0.72648798 0.6816108 0.7551084 >?? [5,] 0.35528206 0.7962423 0.9827632 >?? [6,] 0.09552328 0.3458950 0.1327465 >?? [7,] 0.76918217 0.5945065 0.4511735 >?? [8,] 0.40125232 0.1066155 0.1269416 >?? [9,] 0.56923379 0.2981261 0.9285374 >? [10,] 0.04276991 0.4693320 0.3139434 >> fit2$coef >?????? Contrasts >????????????? day2?????? day3???????? day4 >?? [1,]? 0.3737307? 0.3054398? 0.921284513 >?? [2,]???????? NA???????? NA?????????? NA >?? [3,] -0.8628711 -0.3056396 -0.227256710 >?? [4,] -0.1553395? 0.1821986? 0.138509886 >?? [5,] -0.4107858? 0.1146613 -0.009593154 >?? [6,]? 0.7421089 -0.4188821? 0.668951727 >?? [7,] -0.1303084 -0.2364268 -0.334736056 >?? [8,] -0.3729581 -0.7182350 -0.679192601 >?? [9,]? 0.2528092? 0.4624636? 0.039823091 >? [10,] -0.9030535 -0.3214420 -0.447554386 > > >> sessionInfo() >R version 3.0.0 (2013-04-03) >Platform: x86_64-w64-mingw32/x64 (64-bit) > >locale: >[1] LC_COLLATE=English_United States.1252? LC_CTYPE=English_United States.1252??? LC_MONETARY=English_United States.1252 LC_NUMERIC=C????????????????????????? >[5] LC_TIME=English_United States.1252??? > >attached base packages: >[1] grDevices datasets? splines?? graphics? stats???? tcltk???? utils???? methods?? base???? > >other attached packages: >[1] limma_3.16.3??? svSocket_0.9-55 TinnR_1.0-5???? R2HTML_2.2.1??? Hmisc_3.10-1.1? survival_2.37-4 > >loaded via a namespace (and not attached): >[1] cluster_1.14.4? grid_3.0.0????? lattice_0.20-15 svMisc_0.9-69?? tools_3.0.0??? >> > >_____________________________________________________________________ _ >The information in this email is confidential and intended solely for the addressee. >You must not disclose, forward, print or use it without the permission of the sender. >_____________________________________________________________________ _ > >