#' Compute quantile statistics
#'
#' @param exprs for counts use log2(raw counts + 1)), for MA use log2(raw intensities)
#' @param groups groups to which samples belong (character vector)
#' @param window window size for running median as a fraction on the number of rows of exprs
qstats = function (exprs, groups, window) {
# Compute sample quantiles
Q = apply(exprs, 2, sort)
# Compute quantile reference
Qref = rowMeans(Q)
# Compute SST
SST = rowSums((Q - Qref)^2)
# Compute SSB
f = factor(as.character(groups))
X = model.matrix(~ 0 + f)
QBETAS = t(solve(t(X) %*% X) %*% t(X) %*% t(Q))
Qhat = QBETAS %*% t(X)
SSB = rowSums((Qhat - Qref)^2)
# Compute weights
roughWeights = 1 - SSB / SST
roughWeights[SST < 1e-6] = 1
# Compute smooth weights
k = floor(window * nrow(Q))
if (k %% 2 == 0) k = k + 1
smoothWeights = runmed(roughWeights, k = k, endrule="constant")
list(Q=Q, Qref=Qref, Qhat=Qhat, QBETAS=QBETAS, SST=SST, SSB=SSB, SSE=SST-SSB,
roughWeights=roughWeights, smoothWeights=smoothWeights)
}