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
ADD COMMENT
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}}
ADD COMMENT
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. >_____________________________________________________________________ _ > >
ADD REPLY

Login before adding your answer.

Traffic: 168 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