"sagat" <- function(data,U,nEG=ncol(U),design=NULL,contrasts=1,ebayes=TRUE,orig.coef=TRUE) { if(nrow(data) != nrow(U)) stop("Dimensions of data and U aren't compatible!") N <- nrow(data) # Impute missing data if necessary if(any(is.na(data))) { require(EMV,quietly=TRUE) is.row.na <- apply(data,1,function(x) all(is.na(x))) data[is.row.na,] <- matrix(colMeans(data,na.rm=TRUE),nrow=sum(is.row.na),ncol=ncol(data),byrow=TRUE) data[!is.row.na,] <- knn(data[!is.row.na,],k=10)[[1]] } # Project data into eigenspace A <- t(U[,1:nEG]) %*% data # Calculate eigengene linear model EG.lm <- contrasts.fit(lmFit(A,design=design),contrasts=contrasts) # Calculate gene coefficients if(orig.coef) coef <- contrasts.fit(lmFit(data,design=design),contrasts=contrasts)$coefficients else coef <- U[,1:nEG] %*% EG.lm$coef # Calculate t-statistic if(ncol(data)==1) t <- coef / sqrt(U[,1:nEG]^2 %*% abs(EG.lm$coefficients)) else { if(ebayes) { EG.lm <- eBayes(EG.lm) s2 <- EG.lm$s2.post } else s2 <- EG.lm$sigma^2 t <- coef / sqrt(U[,1:nEG]^2 %*% (s2*EG.lm$stdev.unscaled^2)) } return(as.numeric(t)) }