File Coverage

blib/lib/Statistics/Test/WilcoxonRankSum.pm

Criterion	Covered	Total	%
statement	306	352	86.9
branch	63	84	75.0
condition	14	21	66.6
subroutine	42	46	91.3
pod	15	15	100.0
total	440	518	84.9

package Statistics::Test::WilcoxonRankSum;

use warnings;
use strict;
use Carp;
use Carp::Assert;

use version; our $VERSION = qv('0.0.7');

use Contextual::Return;
use List::Util qw(sum);
use Set::Partition;
use Math::BigFloat;
use Math::Counting ':big';
use Statistics::Distributions;

use Class::Std;

{
  ############ Data ######################################################################

  my %EXACT_UPTO  : ATTR( :init_arg :default<20> );
  my %dataset1_of : ATTR( :get ); # array of numbers
  my %dataset2_of : ATTR( :get ); # array of numbers
  my %n1_of       : ATTR( :get );       # number of elements in dataset 1
  my %n2_of       : ATTR( :get );       # number of elements in dataset 2
  my %N_of        : ATTR( :get  );       # overall number of elements (ranks)
  my %MaxSum_of   : ATTR( :get ); # biggest possible ranksum
  my %ranks_of    : ATTR( :get    :set  ); # hash with ranked data
  my %rank_array_of : ATTR( :get ); # rank array from %ranks
  my %rankSum1_of : ATTR( :get ); # rank sum for dataset 1
  my %expected_rank_sum_1_of : ATTR( :get); # expected rank sum for dataset 1
  my %expected_rank_sum_2_of : ATTR( :get); # expected rank sum for dataset 2
  my %rankSum2_of : ATTR( :get ); # rank sum for dataset 2
  my %smaller_rank_sum_of : ATTR;
  my %smaller_ranks_count_of : ATTR;
  my %expected_rank_count_for_smaller_ranks_count_of : ATTR( :get);
  my %smaller_rank_sums_count_of : ATTR;            # number of possible arrangements with lesser rank sum
                                                    # than the smaller rank sum
  my %rank_sums_other_than_expected_count_of : ATTR; # number of possible arrangements with rank sum
                                                    # other than the smaller rank sum
  my %probability_of : ATTR;                   # probability for the ranking with smaller rank sum
  my %probability_normal_approx_of : ATTR;

  ############ Utility subroutines #######################################################

  sub _check_dataset {
    my ($dataset_ref) = @_;

    croak "Need array ref to dataset\n"
      unless ($dataset_ref);

    croak "Datasets must be passed as array references\n"
      unless (ref($dataset_ref) eq 'ARRAY');

    my @dataset = grep { $_ > 0 } @{ $dataset_ref };
    croak "dataset has no element greater 0\n" unless (@dataset);

    return \@dataset;

  }

  sub _compute_N_MaxSum {
    my ($id) = @_;

    my $N;
    unless ($N_of{$id}) {
      $N = $n1_of{$id} + $n2_of{$id};
      $N_of{$id} = $N;
    }

    unless ($MaxSum_of{$id}) {
      $MaxSum_of{$id} = $N*($N+1)/2;
    }

    unless ($expected_rank_sum_1_of{$id}) {
      $expected_rank_sum_1_of{$id} = $n1_of{$id}*$N/2;
    }

    unless ($expected_rank_sum_2_of{$id}) {
      $expected_rank_sum_2_of{$id} = $n2_of{$id}*$N/2;
    }

    return;
  }



  sub _reset_dependant_datastructures {
    my ($id) = @_;

    delete $ranks_of{$id};
    delete $rank_array_of{$id};
    delete $rankSum1_of{$id};
    delete $rankSum2_of{$id};
    delete $N_of{$id};
    delete $MaxSum_of{$id};
    delete $smaller_rank_sum_of{$id};
    delete $smaller_rank_sums_count_of{$id};
    delete $probability_of{$id};
    delete $probability_normal_approx_of{$id};
    delete $expected_rank_sum_1_of{$id};
    delete $expected_rank_sum_2_of{$id};

    return;
  }

  sub _rank_sum_for {
    my ($self, $dataset) = @_;

    my $id = ident $self;

    my @rank_array;

    if ($rank_array_of{$id} and @{ $rank_array_of{$id} }) {
      @rank_array = @{ $rank_array_of{$id} };
    } else {
      @rank_array = $self->compute_rank_array();
    }

    return sum map { $_->[0] } grep { $_->[1] eq $dataset } @rank_array;
  }


  sub _set_smaller_rank_for {
    my ($id, $rank_sum_1, $rank_sum_2) = @_;
    if ($rank_sum_1 <= $rank_sum_2) {
      $smaller_rank_sum_of{$id} = $rank_sum_1;
      $smaller_ranks_count_of{$id} = $n1_of{$id};
      $expected_rank_count_for_smaller_ranks_count_of{$id} = $expected_rank_sum_1_of{$id};
    } else {
      $smaller_rank_sum_of{$id} = $rank_sum_2;
      $smaller_ranks_count_of{$id} = $n2_of{$id};
      $expected_rank_count_for_smaller_ranks_count_of{$id} = $expected_rank_sum_2_of{$id};
    };
    return;
  };

  sub _NormalZ {  # ($Z) -> $p
    my ($x) = @_;
    #
    # P(x) = 1 - Z(x)(b1*t+b2*t**2+b3*t**3+b4*t**4+b5*t**5)
    # Z(x) = exp(-$x*$x/2.0)/(sqrt(2*3.14159265358979323846))
    # t = 1/(1+p*x)
    #
    # Parameters
    my @b = (0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429);
    my $p = 0.2316419;
    my $t = 1/(1+$p*$x);
    # Initialize variables
    my $fact = $t;
    my $Sum;
    # Sum polynomial
    foreach  my $bi (@b) {
      $Sum += $bi*$fact;
      $fact *= $t;
    };
    # Calculate probability
    $p = 2*$Sum*exp(-$x*$x/2.0)/(sqrt(2*3.14159265358979323846));
    #
    return $p;
  };

  ############ Methods ###################################################################

  sub set_dataset1 {
    my ($self, $dataset1_ref) = @_;

    $dataset1_ref = _check_dataset($dataset1_ref);

    my $id = ident $self;
    $dataset1_of{$id} = $dataset1_ref;
    $n1_of{$id} = scalar(@{ $dataset1_ref });

    _reset_dependant_datastructures($id);

    return;
  }

  sub set_dataset2 {
    my ($self, $dataset2_ref) = @_;

    $dataset2_ref = _check_dataset($dataset2_ref);

    my $id = ident $self;
    $dataset1_of{$id} = $dataset2_ref;
    $n2_of{$id} = scalar(@{ $dataset2_ref });

    _reset_dependant_datastructures($id);

    return;
  }

  sub load_data {
    my ($self, $dataset1_ref, $dataset2_ref) = @_;

    $dataset1_ref = _check_dataset($dataset1_ref);
    $dataset2_ref = _check_dataset($dataset2_ref);

    my $id = ident $self;

    $dataset1_of{$id} = $dataset1_ref;
    $dataset2_of{$id} = $dataset2_ref;
    $n1_of{$id} = scalar(@{ $dataset1_ref });
    $n2_of{$id} = scalar(@{ $dataset2_ref });

    _reset_dependant_datastructures($id);

    _compute_N_MaxSum($id);

    return;
  }

  sub compute_ranks {
    my ($self) = @_;
    my $id = ident $self;

    croak "Please set/load datasets before computing ranks\n" unless ($dataset1_of{$id} and $dataset2_of{$id});

    my @dataset1 = @{ $dataset1_of{$id} };
    my @dataset2 = @{ $dataset2_of{$id} };

    # at this point we are sure we have both data sets, so we may as well compute N and MaxSum - if not already computed
    _compute_N_MaxSum($id);

    my %ranks;

    foreach my $el (@dataset1) {
      $ranks{$el}->{in_dataset}->{ds1}++;
    }
    foreach my $el (@dataset2) {
      $ranks{$el}->{in_dataset}->{ds2}++;
    }

    my $rank=0;
    foreach my $value (sort { $a <=> $b } keys %ranks) {

      my $tied_ranks;

      foreach my $ds (keys %{ $ranks{$value}->{in_dataset} }) {
        $tied_ranks += $ranks{$value}->{in_dataset}->{$ds};
      }

      assert $tied_ranks if DEBUG;

      my $rs;
      for my $r ($rank+1 .. $rank+$tied_ranks) {
        $rs += $r;
      }
      $ranks{$value}->{rank} = $rs/$tied_ranks;
      $ranks{$value}->{tied} = $tied_ranks;

      $rank+=$tied_ranks;
    }

    $ranks_of{$id} = \%ranks;

    return $ranks_of{$id};
  }

  sub compute_rank_array {
    my ($self) = @_;
    my $id = ident $self;

    my @rank_array;
    if ($rank_array_of{$id} and @{ $rank_array_of{$id} } ) {
      @rank_array = @{ $rank_array_of{$id} };
    } else {

      my %ranks;

      if ($ranks_of{$id} and %{ $ranks_of{$id} } ) {
        %ranks = %{ $ranks_of{$id} };
      } else {
        %ranks = %{ $self->compute_ranks() };
      }

      foreach my $value (sort { $a <=> $b } keys %ranks) {
        foreach my $ds (keys %{ $ranks{$value}->{in_dataset} }) {
          for (1..$ranks{$value}->{in_dataset}->{$ds}) {
            push(@rank_array, [ $ranks{$value}->{rank}, $ds ]);
          }
        }
      }

      $rank_array_of{$id} = \@rank_array;

    }

    return (
            SCALAR { scalar @rank_array }      # How many?
            LIST   { @rank_array        }      # What are they?
           );
  }

  sub rank_sum_for {
    my ($self, $for_dataset) = @_;

    my $id = ident $self;

    my $rankSum;
    if ($for_dataset =~ m{1}) {
      if ($rankSum1_of{$id}) {
        return $rankSum1_of{$id};
      } else {
        $rankSum1_of{$id} = $self->_rank_sum_for('ds1');
        return $rankSum1_of{$id};
      }
    } elsif ($for_dataset =~ m{2}) {
      if ($rankSum2_of{$id}) {
        return $rankSum2_of{$id};
      } else {
        $rankSum2_of{$id} = $self->_rank_sum_for('ds2');
        return $rankSum2_of{$id};
      }
    } else {
      croak "Argument must match `1' or `2' (meaning dataset 1 or 2)\n";
    }

    return;

  }

  sub get_smaller_rank_sum {
    my ($self) = @_;

    my $id = ident $self;

    if ($smaller_rank_sum_of{$id} and $smaller_ranks_count_of{$id}) {

      return (
              SCALAR { $smaller_rank_sum_of{$id} }                                  # only the rank sum itselt
              LIST   { ($smaller_rank_sum_of{$id}, $smaller_ranks_count_of{$id} ) } # also the size of the corresponding ds
             );
    }

    my $rank_sum_1 = $rankSum1_of{$id};
    my $rank_sum_2 = $rankSum2_of{$id};

    if (not($rank_sum_1) and not($rank_sum_2)) {
      $rank_sum_1 = $self->rank_sum_for('ds1');
    }


    if ($rank_sum_1 and $rank_sum_2) {

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);

    } elsif ($rank_sum_1) {
      $rank_sum_2 = $MaxSum_of{$id} - $rank_sum_1;
      $rankSum2_of{$id} = $rank_sum_2;

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);

    } elsif ($rank_sum_2) {
      $rank_sum_1 = $MaxSum_of{$id} - $rank_sum_2;
      $rankSum1_of{$id} = $rank_sum_1;

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);

    }

    return (
            SCALAR { $smaller_rank_sum_of{$id} }                                  # only the rank sum itselt
            LIST   { ($smaller_rank_sum_of{$id}, $smaller_ranks_count_of{$id} ) } # also the size of the corresponding ds
           );


    return $smaller_rank_sum_of{$id};
  }

  sub smaller_rank_sums_count {
    my ($self) = @_;
    my $id = ident $self;

    if ($smaller_rank_sums_count_of{$id}) {
      return $smaller_rank_sums_count_of{$id};
    };

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $MaxSum = $MaxSum_of{$id};

    my @ranks = map { $_->[0] } $self->compute_rank_array();

    # let's do some checks before starting the big counting
    if ($W > $MaxSum) { croak "Rank sum bound $W is bigger than the maximum possible rank sum $MaxSum\n" };
    if ($N != scalar(@ranks))
      { croak "Sum of $nA and $nB must be equal to number of ranks: ".scalar(@ranks)."\n" };

    # compute all possible partitions
    my $s = Set::Partition->new(
                                list => \@ranks,
                                partition => [$nA, $nB],
                               );

    my $count_less_W = 0;

    while (my $p = $s->next()) {
      my @pA = @{ $p->[0] };
      my $sumA = sum @pA;
      if ($sumA <= $W) {
        $count_less_W++;
      }
    }

    return $count_less_W;

  };

  sub rank_sums_other_than_expected_counts {
    my ($self) = @_;
    my $id = ident $self;

    if ($rank_sums_other_than_expected_count_of{$id}) {
      return $rank_sums_other_than_expected_count_of{$id};
    };

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $W_exp = $self->get_expected_rank_count_for_smaller_ranks_count();

    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $MaxSum = $MaxSum_of{$id};

    my @ranks = map { $_->[0] } $self->compute_rank_array();

    # let's do some checks before starting the big counting
    if ($W > $MaxSum) { croak "Rank sum bound $W is bigger than the maximum possible rank sum $MaxSum\n" };
    if ($N != scalar(@ranks))
      { croak "Sum of $nA and $nB must be equal to number of ranks: ".scalar(@ranks)."\n" };

    # compute all possible partitions
    my $s = Set::Partition->new(
                                list => \@ranks,
                                partition => [$nA, $nB],
                               );

    my $count_other_W = 0;

    if ($W >= $W_exp) {

      while (my $p = $s->next()) {
	my @pA = @{ $p->[0] };
	my $sumA = sum @pA;
	if ($sumA >= $W) {
	  $count_other_W++;
	}
      }

    } else {

      while (my $p = $s->next()) {
	my @pA = @{ $p->[0] };
	my $sumA = sum @pA;
	if ($sumA <= $W) {
	  $count_other_W++;
	}
      }

    }

    return $count_other_W;

  };


  sub probability :NUMERIFY {
    my ($self) = @_;
    my $id = ident $self;

    if ($probability_of{$id}) {
      return $probability_of{$id};
    }

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};

    my $p;
    if ($N <= $EXACT_UPTO{$id}) {
      $p = $self->probability_exact();
    } else {
      $p = $self->probability_normal_approx();
    }

    $probability_of{$id} = $p;

    return $probability_of{$id};
  }

  sub probability_exact {
    my ($self) = @_;
    my $id = ident $self;

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};

    my $partition_count = bcomb($N, $nA);
    my $have_smaller_rank_sums = $self->rank_sums_other_than_expected_counts();
    my $p = Math::BigFloat->new($have_smaller_rank_sums) * 2.0 / Math::BigFloat->new($partition_count);

    if ($p > 1) { $p = 1 };

    return $p;
  }

  sub probability_normal_approx {
    my ($self) = @_;
    my $id = ident $self;

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $mean = $nA*($N+1)/2;
    my $deviation = sqrt($nA*$nB*($N+1)/12.0);
    my $continuity = (($W - $mean) >= 0) ? -0.5 : +0.5;
    my $z = ($W - $mean + $continuity)/$deviation;
    @{ $probability_normal_approx_of{$id} }{'mean', 'std deviation', 'z'} = ($mean, $deviation, $z);

    my $p = 2*Statistics::Distributions::uprob(abs($z));

    return $p;

  }


  sub probability_status {
    my ($self) = (@_);
    my $id = ident $self;

    my $return_string;
    if ($probability_of{$id}) {
      if ($probability_normal_approx_of{$id}) {
        $return_string = sprintf "Probability: %10f, normal approx w. mean: %10f, std deviation: %10f, z: %10f", $probability_of{$id}, map { $probability_normal_approx_of{$id}->{$_} } ('mean', 'std deviation', 'z');
      } else {
        $return_string = sprintf "Probability: %10f, exact", $probability_of{$id};
      }
    } else {
      $return_string = "Probability not yet computed";
    }

    return (
            STR    { "$return_string" }
            VOID   { print $return_string."\n" }
           );
  }

  sub as_hash :HASHIFY {
    my ($self) = @_;
    my $id = ident $self;

    return {
            dataset_1 => $dataset1_of{$id},
            dataset_2 => $dataset2_of{$id},
            n1 => $n1_of{$id},
            n2 => $n2_of{$id},
            N => $N_of{$id},
            rank_array => $rank_array_of{$id},
            rank_sum_1 => $rankSum1_of{$id},
            rank_sum_2 => $rankSum2_of{$id},
            rank_sum_1_expected => $expected_rank_sum_1_of{$id},
            rank_sum_2_expected => $expected_rank_sum_2_of{$id},
            probability => $probability_of{$id},
            probability_normal_approx => $probability_normal_approx_of{$id},
           };

  }


  sub summary :STRINGIFY {
    my ($self) = (@_);
    my $id = ident $self;

    my $hash = $self->as_hash();

    my $return_string;
    if (not($hash->{dataset_1})) {
      $return_string = "Dataset 1 is not yet initialised, no computations could be done\n";
    } elsif (not($hash->{dataset_2})) {
      $return_string = "Dataset 2 is not yet initialised, no computations could be done\n";
    } else {
      my $format = <
----------------------------------------------------------------
dataset |    n      | rank sum: observed / expected
----------------------------------------------------------------
   1    |%7d    |           %7d      /%7d
----------------------------------------------------------------
   2    |%7d    |           %7d      /%7d
----------------------------------------------------------------
N (size of both datasets): %7d
%s
END_FORMAT
      my $prob = $self->probability_status();
      $return_string = sprintf $format, @{ $hash }{'n1', 'rank_sum_1', 'rank_sum_1_expected', 'n2', 'rank_sum_2', 'rank_sum_2_expected', 'N'}, $prob;
      if ($hash->{probability} >= 0.05) {
        $return_string.="Not significant (at 0.05 level)\n";
      } else {
        $return_string.="Significant (at 0.05 level)\n";
        $return_string.= $hash->{rank_sum_1} > $hash->{rank_sum_1_expected} ?
          "Ranks of dataset 1 are higher than expected\n"
            : "Ranks of dataset 1 are lower than expected\n";

      }
      if ($hash->{N} < 5) {
        $return_string.="Warning: sample size ($hash->{N}) too small (<5)!\n";
      }
    }

    return (
            STR    { "$return_string" }
            VOID   { print $return_string }
           );
  }

}

1; # Magic true value required at end of module


__END__