I will try to describe briefly what I am doing:

my design matrix

designMatrix <- model.matrix( ~ 0 + Condition + Plant, data=grp2$annotation_)

> unique(grp2$annotation_$Condition)

[1] "PC_16h_sys" "PC_16h_test" "PC_38h_sys" "PC_38h_test" "PC_96h_sys" "PC_96h_test"

Then I define the contrasts I want to compute:

contrasts = c( "PC_38h_sys - PC_16h_sys" , "PC_96h_sys - PC_16h_sys" , "PC_38h_test - PC_16h_test" , "PC_96h_test - PC_16h_test" )

Then I generate a fit with lmFit

fit <- lmFit(intmat , designMatrix)

since in the variable intmat in row 38 are missing values for condition PC_96h_test

> fit$coefficients[38,] PC_16h_sys PC_16h_test PC_38h_sys PC_38h_test PC_96h_sys PC_96h_test PlantB2 -4.0764436 -3.5428944 -3.4582324 -2.6018461 -3.4004178 NA -0.2996508

next, I do compute contrasts and all the contrasts are NAs. I would have expected that only the contrast involving PC_96h_test are NA but not all contrasts.

cont <- limma::makeContrasts(contrasts = contrasts, levels = levels) lmfit.cont <- contrasts.fit(fit, cont) > lmfit.cont[38,] An object of class "MArrayLM" $coefficients Contrasts PC_38h_sys - PC_16h_sys PC_96h_sys - PC_16h_sys PC_38h_test - PC_16h_test PC_96h_test - PC_16h_test CCnCP_000456 NA NA NA NA

I very much like the contrasts interface and I much prefer it to working with relevel. Is there a way to better handle NA's with contrasts?

Looking at the code for contrasts.fit I see that there is the matrix multiplication:

fit$coefficients[38,] %*% cont

Which gives NA for all contrasts if the coefficient row contains just a single NA.

Doing something on the lines :

fit$coefficients[38,][!is.na(fit$coefficients[38,])] %*% cont[!is.na(fit$coefficients[38,]),] Contrasts PC_38h_sys - PC_16h_sys PC_96h_sys - PC_16h_sys PC_38h_test - PC_16h_test PC_96h_test - PC_16h_test [1,] 0.6182112 0.6760258 0.9410483 3.542894

could give me the coefficients I am looking for. But of course, there is much more going on in contrasts.fit

cov.coefficients, stdev.unscaled,

Can anyone more informed, tell me if contrasts.fit could be rewritten to better handle NAs for instance by going through the coefficients ,stdev.unscaled row by row dropping NAs etc. And just for completeness

Thank you

Witold

Thanks, Gordon,

I have now a solution based on the post:

contrasts in limma and missing values

but it is a hack, workaround, and features usually do not need workarounds.

Still, if you think that contrast.fit not handling NA is a feature why not just add contrast.fit.NA to limma which is able to handle NA's if this is what is asked for?

Well, I don't consider one request in the past 9 years to be an overwhelming clamour of user requests. Believe me, one can't implement everything than everyone suggests and still write quality software.

For your part, if you want to comment on my answer, please post a comment (by clicking on "Add Comment") rather than answering your own question (by clicking on "Add your answer"). I'm a site moderator, so I've taken the liberty of moving your answer to be a comment.