When 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]
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%.