[BioC] limma: NA coefficients due to missing values
Shi, Tao
shidaxia at yahoo.com
Thu May 16 00:05:31 CEST 2013
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.
>______________________________________________________________________
>
>
More information about the Bioconductor
mailing list