#####################################################################
# A Library of R routines to compute Cui's Shrinkage-Based F-statistic 
#
# Description: 
# This file provides a set of R routines to compute the 
# shrinkage-based one-way ANOVA statistic proposed by Cui et al (2005).
# Instead of relying on simulation-based estimates of the constants
# V and B described by Cui et al (2005), this libary uses the values 
# of those constants that Pounds (2007) derived analytically.
# 
# Terms for Use:  
# This software may be freely used provided that Pounds (2007)
# and Cui et al. (2005) are cited in all resulting presentations 
# and publications.  See specific citation details below.
#
# References:
# 
# Pounds (2007) Computational Enhancement of a Shrinkage-Based
# ANOVA F-test Proposed for Differential Gene Expression Analysis.
# Biostatistics, to appear.
#
# Cui et al. (2005) Improved statistical tests for differential gene
# expression by shrinking variance components estimates.
# Biostatistics 6: 59-75.
#           
# Copyright 2007 Stanley B. Pounds, All Rights Reserved. 
#####################################################################

###################################################################
# Function: array.cui.Fstat
# Purpose: Compute Cui's shrinkage-based one-way ANOVA F-statistic
# Arguments: grps - vector of group labels
#            Y - expression matrix (rows=genes, cols=subjects)
# Returns: a vector of F-statistics
# Reference: Cui et al (2005) Biostatistics 6: 59-75.
###################################################################

array.cui.Fstat<-function(grps,Y)

{
	# Some preliminary values
	ugrp<-unique(grps)  # Determine unique group labels
	k<-length(ugrp)     # Determine number of groups
	g<-dim(Y)[1]        # Determine number of genes
	n<-dim(Y)[2]        # Determine overall sample size

	# Initialize matrix for group means and vector for SS within
	grpmns<-matrix(NA,g,k)
	SSW<-rep(0,g)

	# Compute group means and SS within
	for (i in 1:k)
	{
		ths.grp<-(grps==ugrp[i])
		grpmns[,i]<-rowMeans(Y[,ths.grp])
		SSW<-SSW+rowSums((Y[,ths.grp]-grpmns[,i])^2)
	}

	# Compute usual MSE and then shrinkage-based MSE
	MSE<-SSW/(n-k)
	shr.sig<-cui.mse(MSE,n-k)

	# Compute SS total and SS between
	ovlmn<-rowMeans(Y)
	SST<-rowSums((Y-ovlmn)^2)
	SSB<-SST-SSW

	# Compute and return the F-statistic
	cui.F<-(SSB/(k-1))/shr.sig
	return(cui.F)
}

###################################################################
# Function: cui.mse
# Purpose: Compute the denominator of Cui's shrinkage-based MSE estimator
# Arguments: mse.vec - vector of gene-specific MSEs
#            df - the degrees of freedom associated with the MSE
# Returns: a vector of shrinkage-based MSE's
# Reference: Cui et al (2005) Biostatistics 6: 59-75.
###################################################

cui.mse<-function(mse.vec,df)

{
	G<-length(mse.vec)
	log.mse<-log(mse.vec)
	
	log.mse.bar<-mean(log.mse)
	gmn<-exp(log.mse.bar/df)
		
	ss.log.mse<-sum((log.mse-log.mse.bar)^2)
	B<-cui.B(df)
	V<-cui.V(df)
	
	pos.term<-(1-(G-3)*V/ss.log.mse)
	pos.term<-(pos.term>0)*pos.term

	res<-exp(log.mse.bar)*B*exp(pos.term*(log.mse-log.mse.bar))
	return(res)
}

###################################################################
# Function: cui.B
# Purpose: Compute the bias correction factor for Cui's shrinkage-based MSE estimator
# Arguments: df - the degrees of freedom associated with the MSE
# Returns: a vector of shrinkage-based MSE's
# Reference: Cui et al (2005) Biostatistics 6: 59-75.
###################################################

cui.B<-function(df)

{
	return(exp(-log(2)-digamma(df/2)+log(df)))
}

###################################################################
# Function: cui.V
# Purpose: Compute the V-term for Cui's shrinkage-based MSE estimator
# Arguments: df - the degrees of freedom associated with the MSE
# Returns: a vector of shrinkage-based MSE's
# Reference: Cui et al (2005) Biostatistics 6: 59-75.
###################################################

cui.V<-function(df)

{
	E2<-log(2)^2+2*log(2)*digamma(df/2)+trigamma(df/2)+digamma(df/2)^2
	m<-log(2)+digamma(df/2)-log(df)
	res<-E2-2*log(df)*(m+log(df))+log(df)^2-m^2
	return(res)
}