Statistics::Gtest - calculate G-statistic for tabular data

   use Statistics::Gtest;

   $gt = Statistics::Gtest->new($data);

    $degreesOfFreedom = $gt->getDF();
    $gstat = $gt->getG();

    $gt->setExpected($expectedvalues);
    $uncorrectedG = $gt->getRawG();

 Observed values:
   Spiny Hairy Smooth
     95    53    12

 Expected values (for a 9:3:3:1 ratio):
     90    60    10

 use Statistics::Gtest;
 use Statistics::Distributions;

 ...

 my $gt = Statistics::Gtest->new($data);
 my $df = $gt->getDF();
 my $g = $gt->getG();
 my $sig = '.05';
 my $chis=Statistics::Distributions::chisqrdistr ($df,$sig);
 if ($g > $chis) {
   print "$g: Sig. at the $sv level. ($chis cutoff)\n"
 } 

   $g = Statistics::Gtest->new($data);
   $g = new Statistics::Gtest($data);

 * whitespace-delimited string:          "95 53 12"  
 * reference to 1-dimensional array:     [ 95, 53, 12 ]   
 * reference to 2-dimensional array:     [ [ 10, 20 ], [ 20, 15 ] ]
 * external file or filehandle reference.

   $float = $g->getG();

   $float = $g->getRawG();

   $float = $g->getQ();

   $arrayref = $g->getObserved();

   $arrayref = $g->getExpected();

   $g->setExpected($string);
   $g->setExpected($arrayref);
   $g->setExpected($filename);
   $g->setExpected($filehandle);

   $integer = $g->getDF();

   $g->setDF($integer);

   $rowref = $g->getRow(rownum);

   $colref = $g->getCol(colnum);

   $integer = $g->rowSum($index);

   $integer = $g->colSum($index);

   $integer = $g->getRowNum();

   $integer = $g->getColNum();

   $integer = $g->getSumTotal();

#!/usr/bin/perl -w

package Statistics::Gtest;

use strict;
use vars qw($VERSION);
use Carp;
use IO::File;

$VERSION = '0.07';

my $self;

sub new {
   my ($type, $data) = @_;
   my $datahandle;
   if (! $data) {
      _error(1);
   }
   if (@_ eq 3) {
      _error(2);
   }      
   if ($datahandle = _getHandle($data)) {
      $self  = {};
      _initialize($datahandle);
      bless($self, $type);
      $self;
   }  
   else {
      _error(3);
   }
}

sub getG {
   my $self = shift;
   $self->{'G'} = $self->getRawG() / $self->getQ();
   $self->{'G'};
}

sub getRawG {
   my $self = shift;
   $self->{'logsum'} = _sumCellOp();
   $self->{'G'} = 2 * $self->{'logsum'};
   2 * $self->{'logsum'};
}

sub getQ {
   my $self = shift;
   $self->{'q'} = _williamsC();
   $self->{'q'};
}

sub setExpected {
   my ($self, $expData) = @_;
   my $datahandle;
   my @expect;
   my $expectTotal = 0;
   if (@_ eq 3) {
      _error(2);
   }      
   if ($datahandle = _getHandle($expData)) {
      foreach my $row_ref (@{$datahandle}) {
      if (!ref $row_ref) {
         _error(4);
      }
      my @row = @{$row_ref};
         foreach my $cell (@row) {
            _checkNumValidity($cell);
            push(@expect, $cell);
            $expectTotal += $cell;
         }    
         $self->{'expected'} = \@expect;
         # Expected values no longer intrinsic to observed data, so set
         # flag to 0.
         $self->{'intrinsic'} = 0;  
      }
   }
   else {
      _error(3);
   }
   # Data sanity checks
   if ($expectTotal != $self->{'sumtotal'}) {
      warn "Warning: Total of expected values ($expectTotal) does not ",
      "equal total of observed values (",$self->{'sumtotal'},").\nThis ",
      "will invalidate the test unless the discrepancy is very minor ",
      "\n(e.g., the result of rounding error).\n";
   }
}

sub getObserved {
   my $self = shift;
   my $exp = _formatData($self->{'observed'});
   $exp;
}

sub getExpected {
   my $self = shift;
   my $exp = _formatData($self->{'expected'});
   $exp;
}

sub getDF {
   my $self = shift;
   $self->{'df'};
}

sub setDF {
   my $self = shift;
   $self->{'df'} = shift;
   # Degrees of freedom set externally.
   $self->{'df_extern'} = 1;
}

sub getDFstate {
   my $self = shift;
   $self->{'df_extern'};
}

sub getHypothesis {
   my $self = shift;
   $self->{'intrinsic'};
}

sub getRow {
   my $self = shift;
   my $idx = shift;
   if ($idx > ($self->{'nrows'} - 1)) {
      return;
   }   
   my @row = ();
   my $skip = $self->{'ncols'};
   my $offset = $idx * $skip;
   my $final = $offset + $skip;
   for my $cell ($offset .. $final - 1) {
      push(@row, $self->{'observed'}->[$cell]);
   }    
   \@row;
}

sub getCol {
   my $self = shift;
   my $idx = shift;
   if ($idx > ($self->{'ncols'} - 1)) {
      return;
   }   
   my @col = ();
   my $skip = $self->{'ncols'};
   my $max = $skip * $self->{'nrows'};
   for (my $cell = $idx; $cell < $max; $cell += $skip) {
      push(@col,$self->{'observed'}->[$cell]);
   }    
   \@col;
}

sub rowSum {
   my $self = shift;
   my $idx  = shift;
   $self->{'rowsums'}->[$idx];
}

sub colSum {
   my $self = shift;
   my $idx  = shift;
   $self->{'colsums'}->[$idx];
}
   
sub getRowNum {
   my $self = shift;
   $self->{'nrows'};
}

sub getColNum {
   my $self = shift;
   $self->{'ncols'};
}

sub getSumTotal {
   my $self = shift;
   $self->{'sumtotal'};
}

# private methods

# Do all the initialization work: verify the data and write it into the 
# object; calculate the row and column sums and grand total; determine
# the format of the data table (rows x columns); calculate the likely
# degrees of freedom and generate default expected values.  Set state
# flags indicating that the expected values are 'intrinsic' (generated
# from the observed data), and that the degrees of freedom are also
# based on table characteristics.
sub _initialize {
   my $datahandle = shift;
   my @datastruct = ();
   my @rowsums = ();
   my @sumcols = ();
   my ($ncols, $nrows, $sumtotal);
   my $totalN = 0;
   $self->{'colsums'} = ();

   foreach my $row_ref (@{$datahandle}) {
      if (!ref $row_ref) {
         _error(4);
      }
      my @row = @{$row_ref};
      $ncols = scalar @row;
      my $i = 0;
      foreach my $cell (@row) {
         _checkNumValidity($cell);
         push(@datastruct, $cell);
         $sumtotal += $cell;
         $self->{'colsums'}->[$i] += $cell;
         $i++;
         $totalN++;
      }    
      $nrows++;
      push(@rowsums, _rowsum(\@row));
   }

   # Data sanity checks; make sure data exist and are internally 
   # consistent.
   if ($totalN == 0) {
      _error(5);
   }

   if ($totalN != ($nrows * $ncols)) {
      _error(6); 
   }
      
     if (($ncols == 0) || ($nrows == 0)) { 
      _error(7); 
   }
      
   my $tSum = 0;   
   foreach my $r (@rowsums) {
      $tSum += $r;
   }
   if ($sumtotal != $tSum) {
      _error(6); 
   }
      
    $tSum = 0;   
   foreach my $r (@{$self->{'colsums'}}) {
      $tSum += $r;
   }
   if ($sumtotal != $tSum) {
      _error(6); 
   }

   $self->{'observed'} = \@datastruct;
   $self->{'rowsums'} = \@rowsums;
   $self->{'nrows'} = $nrows;
   $self->{'ncols'} = $ncols;
   $self->{'sumtotal'} = $sumtotal;

   # tabletype: 0 = 2-way, 1 = single column, 2 = single row
   if (($ncols == 1) || ($nrows == 1)) {
      if ($ncols == 1) { 
         $self->{'df'} = $nrows - 1; 
         $self->{'tabletype'} = 1;
      }
      if ($nrows == 1) { 
         $self->{'df'} = $ncols - 1; 
         $self->{'tabletype'} = 2;
      }    
   }
   else {
      $self->{'df'} = ($nrows - 1) * ($ncols - 1);
      $self->{'tabletype'} = 0;
   }    
   $self->{'expected'}= _expected(\@datastruct);
   $self->{'intrinsic'} = 1;
   $self->{'df_extern'} = 0;
}

# Calculate and return William's correction, q. Calculation varies 
# depending on table format (1-way vs. 2-way).
sub _williamsC {
   my $sumtotal = $self->{'sumtotal'};
   my $rowsums = $self->{'rowsums'};
   my $colsums = $self->{'colsums'};
   if ($self->{'tabletype'} == 1) {
      $self->{'q'} = 
         1 + ($self->{'nrows'} + 1) / (6 * $sumtotal);
   }
   elsif ($self->{'tabletype'} == 2) {
      $self->{'q'} = 
         1 + ($self->{'ncols'} + 1) / (6 * $sumtotal);
   }
   else {
      my $recipRows = 0;
      my $recipCols = 0;
      foreach my $rval (@{$rowsums}) {
         $recipRows += (1 / $rval);
      }    
      foreach my $cval (@{$colsums}) {
         $recipCols += (1 / $cval);
      }    
      my $denom = 6 * $sumtotal * $self->{'df'};
      my $num = ($sumtotal * $recipRows - 1) * ($sumtotal * $recipCols - 1);
      $self->{'q'} = 1 + $num / $denom;
   }
}

# Generate the default expected values, based on a hypothesis of 
# independence (no association) between factors.
# Return ref to expected array.
sub _expected {
   my $data = shift;
   my @expected = ();
   my $rows = $self->{'nrows'};
   my $cols = $self->{'ncols'};
   my $sumtotal = $self->{'sumtotal'};
   for my $row (0 .. $rows-1) {
      my $marginalRow = $self->{'rowsums'}->[$row] / $sumtotal;    
      my @a;
      for my $col (0 .. $cols-1) {
         my $marginalCol = $self->{'colsums'}->[$col] / $sumtotal;    
         push(@expected, ($marginalRow * $marginalCol * $sumtotal));
      } 
   }
   \@expected;
}

# Calculate the log-likelihood statistic.
# Return sum of Obs * ln(Obs/Exp) for all cells.
sub _sumCellOp {
   my $logsum = 0;
   my $ncells = $self->{'nrows'} * $self->{'ncols'};
   my $o = $self->{'observed'};
   my $e = $self->{'expected'};
   for my $cell (0 .. $ncells-1) {
      $logsum += $o->[$cell] * log($o->[$cell] / $e->[$cell]);    
   }        
   $logsum;
}

# Return sum of the cells in this row.
sub _rowsum {
   my $row = shift;
   my $total = 0;
   foreach my $val (@{$row}) {
      $total += $val;
   }    
   $total;
}

# Order the data array back into the format originally input.
sub _formatData {
   my $ref = shift;
   my @data;
   my $len =  scalar @{$ref};
   my $rlen = $len / $self->{'nrows'};
   my $idx = $rlen - 1;
   for (my $r=0; $r<$len; $r+=$rlen) {
      my $row = [ @{$ref}[$r ... $idx] ];
      return $row if ($rlen eq $len);
      push(@data, $row);
      $idx += $rlen;
   }
   \@data;
}

# Check to make sure the values are numbers, greater than 0.
# Throw error if any problems are found.
sub _checkNumValidity {
   my $value = shift;
   my $eString = "Invalid data.";
   if ($value !~ /^\d+\.?\d*$/) {
      _error(8); 
   }

   # Very low cell counts (< 5) aren't great, either. Maybe a 
   # warning is in order.
   if ($value < 5) {
      warn "Warning: Very small cell frequency (freq. < 5) found.\n",
		"This will reduce the reliability of the test.\n";
   }
   if ($value == 0) {
      _error(9); 
   }
}

# Take a filehandle, return an array reference.
sub _fileread {
   my $data = shift;
   my @a;
   while (<$data>) {
      my @row = split(/\s+/,$_);
      push(@a, \@row);
   }
   \@a;
}

# Try to take whatever input format is thrown at us, and turn it into
# an array reference, or just bail out completely if we can't figure
# it out.
sub _getHandle {
   my $data = shift;
   if (ref $data) {
      if ($data =~ /\AARRAY/) {
         if (ref $data->[0]) {
            return $data;
         }
         else {
            return [ $data ];
         }   
      }
      if ($data =~ /\AGLOB/) {
         return _fileread($data);
      }
      if ($data =~ /\AIO::File=GLOB/) {
         return _fileread($data);
      }
      else {
         _error(10);
      }
   }
   else {
      my $fh = new IO::File;
      if ($fh->open("< $data")) {
         return _fileread($fh);
      }
      else {
         my @a = split(/\s+/,$data);
         return [ \@a ];

      }
      undef $fh;
   }
   0;
}

# Keep all the error strings in one place.
sub _error {
   my ($code) = shift;
	my $errorbar = "\n!-----------------!\n";
   my @e = (
   "", # Normal execution
   "Input error: Must pass in a filename for a file containing the data.\n\tEx: \$g = new Statistics::Gtest(\$datafile);",
   "Input error: Can't use array or hash as data; pass array ref.",
   "Input error: Invalid data format.",
   "Input error: Data must be structured as array of array refs, with each array ref
      pointing to one row of data.",
   "Data inconsistency: No data found in table.",
   "Data inconsistency: sum of table values does not match total.",
   "Data inconsistency: could not read row or column.",
   "Data invalid: All data must be positive integers.",
   "Data invalid: Found a zero frequency value. Can't continue.",
   "Input error: Reference to invalid data type; must be reference to array.",
   );

	print $errorbar;
	print "\n Error",Carp::shortmess();
   print " ", $e[$code], "\n";
	print Carp::longmess(" Stack trace:\n");
	print $errorbar;
   exit $code; 
}

1; 

__END__