lmFit is used with a weights matrix,
contrasts.fit produces a wrong
stdev.unscaled, which is far off from the true value. Here's an example:
Y <- matrix(0, nrow=2, ncol=6) X <- matrix(c(1,1,1,1,0,0,1,1,0,0,1,0,0,0,1,1,0,1), ncol=3) weights <- matrix(c(rep(1,4),.2,.4,rep(1,6)), byrow=T, nrow=2) contr <- c(0,-1,1) fit <- lmFit(Y, X, weights=weights) fit <- contrasts.fit(fit, contrasts=contr) # wrong value fit$stdev.unscaled 1.384713 # correct value Xw <- diag(sqrt(weights[1,])) %*% X sqrt(contr %*% chol2inv(qr(Xw)$qr) %*% contr) 0.9393364
I know the note in
contrasts.fit mentions that in case of precision weights or missing values,
stdev.unscaled is approximate because the design matrix is not re-factorized for every gene, but for the end user it is difficult to judge if a test statistic is well approximated or not. It would be good to display a warning if
contrasts.fit is run in a weights-matrix setting. In my case, the approximation caused a drop in power from 70% to 30%.