I was surprised that limma sometimes returns very small p-values in comparisons where one group is completely missing. And I am sorry if this is already a well-known behavior, I tried to search for descriptions of this, but couldn't find anything that matched my exact observation.

library(limma)
#> Warning: package 'limma' was built under R version 4.2.2

limma_test <- function(y, dat, formula, contrast){
  fit <- lmFit(y, design = model.matrix(formula, data = dat))
  fit <- contrasts.fit(fit, contrasts = contrast)
  fit <- eBayes(fit)
  topTable(fit)
}

dat <- data.frame(cond = rep(c("a", "b", "c"), each = 3))
y_compl <- matrix(c(2, 4, 8, 17, 19, 15, 2, 5, 3), nrow = 1)
limma_test(y_compl, dat, ~ cond, c(0, -1, 0))
#>       logFC  AveExpr         t      P.Value    adj.P.Val          B
#> 1 -12.33333 8.333333 -6.609954 0.0005770005 0.0005770005 -0.9740698
limma_test(y_compl, dat, ~ cond - 1, c(1, -1, 0))
#>       logFC  AveExpr         t      P.Value    adj.P.Val          B
#> 1 -12.33333 8.333333 -6.609954 0.0005770005 0.0005770005 -0.9740698

y_mis <- matrix(c(NA, NA, NA, 17, 19, 15, 2, 5, 3), nrow = 1)
limma_test(y_mis, dat, ~ cond, c(0, -1, 0))
#> Warning: Partial NA coefficients for 1 probe(s)
#>       logFC  AveExpr         t      P.Value    adj.P.Val         B
#> 1 -13.66667 10.16667 -9.406045 0.0007120111 0.0007120111 -1.005365
limma_test(y_mis, dat, ~ cond - 1, c(1, -1, 0))
#> Warning: Partial NA coefficients for 1 probe(s)
#> Warning in .ebayes(fit = fit, proportion = proportion, stdev.coef.lim =
#> stdev.coef.lim, : Estimation of var.prior failed - set to default value
#>   logFC  AveExpr  t P.Value adj.P.Val  B
#> 1    NA 10.16667 NA      NA        NA NA
limma_test(y_mis, dat, ~ cond - 1, c(0, -1, 1))
#> Warning: Partial NA coefficients for 1 probe(s)
#>       logFC  AveExpr         t      P.Value    adj.P.Val         B
#> 1 -13.66667 10.16667 -9.406045 0.0007120111 0.0007120111 -1.005365

(This is using R 4.2.1 and limma 3.54.1)

The first two calls to 'limma_test' show that it usually doesn't matter how the design is specified if the contrast is adapted accordingly. However, if the data contains missing values, in my specific example, limma returns a small p-value that actually corresponds to the comparison b vs. c.

Can someone reproduce this behavior, and is the behavior expected? And lastly, is there a workaround to ensure that, independent of the specification of the design matrix, a contrast between two groups where one is completely missing the p-value is NA?

With any linear model, NA predictors are dropped. As an example

> df <- data.frame(dat, y = t(y_compl))
> df2 <- df
> df2[1:3, 1] <- NA
> summary(lm(y~cond, df))

Call:
lm(formula = y ~ cond, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6667 -1.3333 -0.3333  1.6667  3.3333 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4.667      1.319   3.537 0.012263 *  
condb         12.333      1.866   6.610 0.000577 ***
condc         -1.333      1.866  -0.715 0.501705    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.285 on 6 degrees of freedom
Multiple R-squared:  0.9158,    Adjusted R-squared:  0.8877 
F-statistic: 32.62 on 2 and 6 DF,  p-value: 0.0005976

> summary(lm(y~cond, df2))

Call:
lm(formula = y ~ cond, data = df2)

Residuals:
         4          5          6          7          8          9 
-3.553e-15  2.000e+00 -2.000e+00 -1.333e+00  1.667e+00 -3.333e-01 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   17.000      1.027  16.547 7.81e-05 ***
condc        -13.667      1.453  -9.406 0.000712 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.78 on 4 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.9567,    Adjusted R-squared:  0.9459 
F-statistic: 88.47 on 1 and 4 DF,  p-value: 0.000712

The workaround is to not have NA predictor values, or to realize that they are going to be dropped and planning accordingly.

Yes, it is true that the meaning of the limma coefficients will change if you fit a model with an intercept and all the observations in the reference level of the treatment are missing. This is a consequence of subsetting the design matrix to remove missing y-values from the linear model.

Using your example, with the full design matrix:

> cond = rep(c("a", "b", "c"), each = 3)
> design <- model.matrix(~cond)
> design
  (Intercept) condb condc
1           1     0     0
2           1     0     0
3           1     0     0
4           1     1     0
5           1     1     0
6           1     1     0
7           1     0     1
8           1     0     1
9           1     0     1

the second coefficient corresponds to b vs a. However if the first three observations are removed:

> design[4:9,]
  (Intercept) condb condc
4           1     1     0
5           1     1     0
6           1     1     0
7           1     0     1
8           1     0     1
9           1     0     1

then the design matrix is no longer of full rank. The 3rd coefficient will be removed from the analysis and the second coefficient will compare b to c instead of b to a.

The work around is to use the group mean parametrization

design <- model.matrix(~0+cond)

then this problem will not occur.

I have discussed this issue previously on this forum. The basic rule is that the coefficients cannot be interpreted according to their original meanings if NAs cause the design matrix to become not of full rank. limma issues a warning whenever this occurs. The oneway layout with the group-mean parametrization is an exception for which the coefficients will retain their interpretations.