cvdiag <- function(data,fileprefix="out")
 {
  #############################################################
  # Make CV diagnostic plots for replicate array data
  #
  # R 1.6.1 or later.
  #
  # License: GPL Version 2.0 or later.
  #
  # Author: M.A. Newton  newton@stat.wisc.edu
  #
  # Reference:  Newton, Noueiry, Sarkar, and Ahlquist (2003).
  #     Detecting differential gene expression with a semiparametric
  #     hierarchical mixture method.  Technical Report 1074, Department
  #     of Statistics, UW Madison.  January 2003.
  #     Figures 1 and 2.
  #
  # Web:  www.stat.wisc.edu/~newton/
  #
  #############################################################
  #
  # Input:
  #
  # data        a matrix, 1 column for each replicate, 1 row for each gene
  #             The code assumes there are lots of genes.
  #
  # fileprefix  a character string to name the output postscript files
  #
  # Output:
  #
  #  Two diagnostic plots.  One compares the ranked sample means
  #  to the ranked coefficients of variation.   The second looks
  #  at QQ plots against a common-shape Gamma distribution.
  ###############################################################

  means <- apply(data,1,mean)
  cvs <- (apply(data,1,cv))^2  

  # Rank plot 

  name1 <- paste( fileprefix,"-rank.ps", sep="")
  postscript(file=name1, height=6, width=6, hor=FALSE)

  plot(rank(means), rank(cvs), pch="." ,
     xlab=c("Rank of Mean"), ylab=c("Rank of CV") )

  # Make 9 bins
  mm <- round( length(means) * ppoints(9) )
  mmleft <- mm-50
  mmright <- mm+50
  rr <- c(1,length(means))
  for( i in 1:9 ) {
   polygon( c( rep( mmleft[i], 2 ), rep(mmright[i],2) ), c( rr, rev(rr) ),
            col="pink")
   }
  points(rank(means), rank(cvs), pch="." )
     
 dev.off()

 # Have data in matrix `data'
 lowers <- (sort(means))[mmleft]
 uppers <- (sort(means))[mmright]
 #
 shape.gg <- 1/( mean(cvs) )
 name1 <- paste( fileprefix,"-qq.ps", sep="")
 postscript(file=name1,height=6,width=6, hor=FALSE)
 par(mfrow=c(3,3),mar=c(1.2,1.2,1.2,1.2), mgp=c(3,.5,0))
 for(i in 1:9) {
   data.qq <- c(data[(means>=lowers[i]&means<=uppers[i]),])
   mu <- mean(data.qq)
   scale.gg <- shape.gg/mu
   xx <- qgamma( ppoints(data.qq), shape=shape.gg, rate=scale.gg )
   yy <- sort( data.qq )
   rr <- range( c(xx,yy) )
   plot(xx, yy, xlim=rr, ylim=rr )
   abline(0,1)
 }
 
 dev.off()
 
 }
 
cv <- function(x) { return(sqrt(var(x))/mean(x)) }