library(stats)

friedmans_test = function( DF, block_key, treatment_key, response_key, alpha=0.05, debugging=FALSE ){
  #
  # Returns:
  #   Fr: the Friedman test statistic
  # 
  # See Page 620 from the book Probability and Statistics: For Engineering and the Sciences by Jay L. Devore
  #
  # Written by:
  # -- 
  # John L. Weatherwax                2009-04-21
  # 
  # email: wax@alum.mit.edu
  # 
  # Please send comments and especially bug reports to the
  # above email address.
  #
  #-----

  N = nrow(DF)

  # Get the unique labels of the blocks: 
  block_key_column_index = which( colnames(DF) == block_key )
  block_names = unique( DF[[block_key_column_index]] )
  J = length(block_names)
 
  # Get the unique labels of the treatments: 
  treatment_key_column_index = which( colnames(DF) == treatment_key )
  treatment_names = unique( DF[[treatment_key_column_index]] )
  I = length(treatment_names)

  response_key_column_index = which( colnames(DF) == response_key )

  if( debugging ){
    print( sprintf( "N= %10d", N ) )
    print( sprintf( "I= %10d", I ) )
    print( sprintf( "J= %10d", J ) )
  }

  # Calculate the r_i values:
  # 
  ri = matrix( 0, I, J )
  for( jj in 1:J ){
    inds = DF[[block_key_column_index]]==block_names[jj]
    block_jj_response = DF[inds,response_key]
    block_jj_ranks = rank( block_jj_response )
    ri[,jj] = block_jj_ranks 
  }
  if( debugging ){ 
    print( "ri= " ); print( ri ) 
  }

  Ri = rowSums(ri)
  if( debugging ){
    print( "Ri= " ); print( Ri )
  }

  Fr = ( 12 / ( I*J*(I+1) ) ) * sum( Ri^2 ) - 3 * J * (I+1) 

  f_crit = qchisq( 1 - alpha, I-1 ) 

  # Compute the P-value for this test:
  # 
  p_value = 1 - pchisq( Fr, I-1 )

  list( Ri=Ri, Fr=Fr, f_crit=f_crit, p_value=p_value )
}