File Coverage

blib/lib/Statistics/ChisqIndep.pm

Criterion	Covered	Total	%
statement	97	106	91.5
branch	13	26	50.0
condition	2	6	33.3
subroutine	8	9	88.8
pod	0	4	0.0
total	120	151	79.4

package Statistics::ChisqIndep;

use strict;
use vars qw ($VERSION $AUTOLOAD);
use Carp;
use Statistics::Distributions qw (chisqrprob);

my %fields = (obs => [],
              expected => [],
              rows => 0,
              cols => 0,
              df => 0,
              total => 0,
              rtotals => [],
              ctotals => [],
              chisq_statistic  => 0,
              p_value => 0,
              valid => 0,
              warning => 0);
$VERSION = '0.1';

sub new {
  my $proto = shift;
  my $class = ref($proto) || $proto;
  my $self= {%fields};
  bless($self,$class);
  return $self;
}

sub load_data {
  my $self = shift;
  $self->{valid} = 0;
  $self->{obs} = shift;
  my @obs = @{$self->{obs}};
  my $exp;
  my $rows = scalar(@obs);
  return if $rows < 2;
  my $cols = scalar(@{$obs[0]});
  return if $cols < 2;
  my @rtotals;
  my @ctotals;
  my $total = 0;
  my $df = 0;
  my $chisq = 0;
  my $chiprob = 0;
  for (my $i = 0; $i < $rows; $i++) {
     my $c = scalar(@{$obs[$i]});
     return if $c != $cols;
     for (my $j = 0; $j < $cols; $j++) {
        $rtotals[$i] = 0 if (!defined($rtotals[$i]));
        $ctotals[$j] = 0 if (!defined($ctotals[$j]));
        $rtotals[$i] += $obs[$i]->[$j];
        $ctotals[$j] += $obs[$i]->[$j];
        $total += $obs[$i]->[$j];
        $self->{warning} = 1 if $obs[$i]->[$j] < 25;
     }
  }
  return if ($total < 0);


  #compute the expected values for each cell
  for (my $i = 0; $i < $rows; $i++) {
    for (my $j = 0; $j < $cols; $j++) {
        if ($rtotals[$i] == 0 || $ctotals[$j] == 0) {
            $exp->[$i]->[$j] = 0
        } else {
            $exp->[$i]->[$j] = $rtotals[$i] * $ctotals[$j] / $total;
            $chisq += ($obs[$i]->[$j] - $exp->[$i]->[$j])**2 /$exp->[$i]->[$j];
        }
    }
  }

  $df = ($rows - 1) * ($cols - 1);
  $chiprob = chisqrprob($df, $chisq);

  #update the global fields
  $self->{expected} = $exp;
  $self->{rows} = $rows;
  $self->{cols} = $cols;
  $self->{df} = $df;
  $self->{total} = $total;
  $self->{rtotals} = \@rtotals;
  $self->{ctotals} = \@ctotals;
  $self->{chisq_statistic} = $chisq;
  $self->{p_value} = $chiprob;
  $self->{valid} = 1;

  return 1;
}

sub print_contingency_table {
  my $self = shift;
  my $rows = $self->{rows};
  my $cols = $self->{cols};
  my $obs = $self->{obs};
  my $exp = $self->{expected};
  my $rtotals = $self->{rtotals};
  my $ctotals = $self->{ctotals};
  my $total = $self->{total};

  for (my $j = 0; $j < $cols; $j++) {
    print "\t",$j + 1;
  }
  print "\trtotal\n";
  for (my $i = 0; $i < $rows; $i ++) {
    print $i + 1, "\t";
    for(my $j = 0 ; $j < $cols; $j ++) {
      print $obs->[$i]->[$j], "\t";
    }
    print $rtotals->[$i], "\n";
    print "\t";
    for(my $j = 0 ; $j < $cols; $j ++) {
      printf "(%.2f)\t", $exp->[$i]->[$j];
    }
    print "\n";
  }
  print "ctotal\t";
  for (my $j = 0; $j < $cols; $j++) {
    print $ctotals->[$j], "\t";
  }
  print $total, "\n";
}

sub print_summary {
  my $self = shift;
  if($self->{valid}) {
    print "Rows: ", $self->{rows}, "\n";
    print "Columns: ", $self->{cols}, "\n";
    print "Degree of Freedom: ", $self->{df}, "\n";
    print "Total Count: ", $self->{total}, "\n";
    print "Chi-square Statistic: ", $self->{chisq_statistic}, "\n";
    print "p-value: ", $self->{p_value}, "\n";
    print "Warning: some of the cell counts might be too low. \n" if ($self->{warning});
    $self->print_contingency_table();
  }
}

sub AUTOLOAD {
  my $self = shift;
  my $type = ref($self)
    or croak "$self is not an object";
  my $name = $AUTOLOAD;
  $name =~ s/.*://;     ##Strip fully qualified-package portion
  return if $name eq "DESTROY";
  unless (exists $self->{$name} ) {
    croak "Can't access `$name' field in class $type";
  }
  ##Read only method
  return $self->{$name};
}

1;
__END__