I will try to describe briefly what I am doing:
my design matrix
designMatrix <- model.matrix( ~ 0 + Condition + Plant, data=grp2$annotation_)
 "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
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