File Coverage

blib/lib/Statistics/ANOVA/Page.pm

Criterion	Covered	Total	%
statement	69	81	85.1
branch	4	8	50.0
condition	1	2	50.0
subroutine	17	20	85.0
pod	7	7	100.0
total	98	118	83.0

package Statistics::ANOVA::Page;

use 5.006;
use strict;
use warnings FATAL => 'all';
use base qw(Statistics::Data);
use List::AllUtils qw(sum0);
use Math::Cephes qw(:dists);
use Statistics::Data::Rank;
use Statistics::Zed;

=head1 NAME

Statistics::ANOVA::Page - Nonparametric analysis of variance by ranks for trend across dependent variables (Page and sign tests).

=head1 VERSION

Version 0.01

=cut

our $VERSION = '0.01';

=head1 SYNOPSIS

 use Statistics::ANOVA::Page;
 my $page = Statistics::ANOVA::Page->new();
 $page->load({1 => [2, 4, 6], 2 => [3, 3, 12], 3 => [5, 7, 11, 16]}); # note ordinal datanames
 my $l_value = $page->observed(); # or expected(), variance()
 my ($z_value, $p_value) = $page->zprob_test(ccorr => 2, tails => 1);
 # or without pre-loading:
 $l_value = $page->observed(data => {1 => [2, 4, 6], 2 => [5, 3, 12]});
 # or for subset of loaded data:
 $l_value = $page->observed(lab => [1, 3]);

=head1 DESCRIPTION

Calculates Page statistics for nonparametric analysis of variance across given orders of dependent variables. Ranks are computed exactly as for the L test, but the ranks are weighted according to the ordinal position of the group/level to which they pertain. Also, the test of significance is based on a standardized value, with the I-value read off the normal distribution. Similarly to the relationship between the Kruskal-Wallis and L tests for non-dependent observations, the Friedman test returns the same value regardless of the ordinality of the variables as levels, but the Page test tests ranks in association with the ordinality of the variables (as levels rather than groups). These are weighted according to their Perl sort { $a <=> $b} order, so they should have sort-able names that reflect the ordering of the variables.

With only two groups, the test statistic is equivalent to that provided by a B.

Build tests include comparison of return values with published data, viz. from Hollander and Wolfe (1999, p. 286ff); passing these tests means the results agree.

=head1 SUBROUTINES/METHODS

=head2 new

 $page = Statistics::ANOVA::Page->new();

New object for accessing methods and storing results. This "isa" Statistics::Data object.

=head2 load, add, unload

 $page->load(1 => [1, 4], 2 => [3, 7]);

The given data can now be used by any of the following methods. This is inherited from L, and all its other methods are available here via the class object. Only passing of data as a hash of arrays (HOA) is supported for now. Alternatively, give each of the following methods the HOA for the optional named argument B.

=head2 observed

 $val = $page->observed(); # data pre-loaded
 $val = $page->observed(data => $hashref_of_arefs);

Returns the observed statistic I based on within-group rankings of the data weighted according to the ordinal position of the variable (by its numerical name) to which they pertain.

Optionally, if the data have not been pre-loaded, send as named argument B.

=cut

sub observed {
    my ( $self, %args ) = @_;
    return _calc_l_value( _get_data( $self, %args ) );
}

=head2 observed_r

 $val = $page->observed_r(); # data pre-loaded
 $val = $page->observed_r(data => $hashref_of_arefs);

This implements a "l2r" transformation: Hollander and Wolfe (1999) describe how Page's I-statistic is directly related to Spearman's rank-order correlation coefficient (see L), based on the observed and predicted order of each associated group/level per observation.

=cut

sub observed_r {
    my ( $self, %args ) = @_;
    return _calc_l2r_value( _get_data( $self, %args ) );
}

=head2 expected

 $val = $page->expected(); # data pre-loaded
 $val = $page->expected(data => $hashref_of_arefs);

Returns the expected value of the I statistic for the given data.

=cut

sub expected {
    my ( $self, %args ) = @_;
    return _calc_l_exp( _get_data( $self, %args ) );
}

=head2 variance

 $val = $page->variance(); # data pre-loaded
 $val = $page->variance(data => $hashref_of_arefs);

Return the variance expected to occur in the I values for the given data.

=cut

sub variance {
    my ( $self, %args ) = @_;
    return _calc_l_var( _get_data( $self, %args ) );
}

=head2 zprob_test

 $p_val = $page->zprob_test(); # data pre-loaded
 $p_val = $page->zprob_test(data => $hashref_of_arefs);
 ($z_val, $p_val) = $page->zprob_test(); # get z-score too

Calculates an expected I value and variance, to provide a normalized I for which the I-value is read off the normal distribution. This is appropriate for "large" samples. Optional arguments are B and B as in L.

=cut

sub zprob_test {
    my ( $self, %args ) = @_;
    my $href = _get_data( $self, %args );
    my $zed = Statistics::Zed->new();
    my ( $z_value, $p_value ) = $zed->z_value(
        observed => _calc_l_value($href),
        expected => _calc_l_exp($href),
        variance => _calc_l_var($href),
        %args
    );
    return wantarray ? ( $z_value, $p_value ) : $p_value;
}

=head2 chiprob_test

 $p_val = $page->chiprob_test(); # data pre-loaded
 $p_val = $page->chiprob_test(data => $hashref_of_arefs);
 ($chi_val, $p_val) = $page->chiprob_test(); # get z-score too

Calculates a chi-square statistic based on the observed value of , the number of ranked variables, and the number of replications; as per Page(1963, Eq. 4). This is a two-tailed test; if the optional argument B => 1, the returned probability, read off the chi-square distribution, is halved.

=cut

sub chiprob_test {
    my ( $self, %args ) = @_;
    my $data_href = _get_data( $self, %args );
    my $l         = _calc_l_value($data_href);
    my $n_bt      = scalar keys %{$data_href};
    my $n_wt      = __PACKAGE__->equal_n( data => $data_href );
    my $num = ( ( 12 * $l ) - ( 3 * $n_wt * $n_bt ) * ( $n_bt + 1 )**2 )**2;
    my $chi =
      $num / ( ( $n_wt * $n_bt**2 ) * ( $n_bt**2 - 1 ) * ( $n_bt + 1 ) );
    my $p_value = chdtrc( 1, $chi );    # Math::Cephes fn
    $args{'tails'} ||= 2;
    $p_value /= 2 if $args{'tails'} == 1;
    return wantarray ? ( $chi, 1, ( $n_bt * $n_wt ), $p_value ) : $p_value;
}

=head2 chiprob_str

 $str = $page->chiprob_str(data => HOA, correct_ties => 1);

Performs the same test as for L but returns not an array but a string of the conventional reporting form, e.g., chi^2(df, N = total observations) = chi_value, p = p_value.

=cut

sub chiprob_str {
    my ( $self, %args ) = ( shift, @_ );
    my ( $chi_value, $df, $count, $p_value ) = $self->chiprob_test(%args);
    return "chi^2($df, N = $count) = $chi_value, p = $p_value";
}

sub _calc_l_value {
    my $data  = shift;
    my $ranks = Statistics::Data::Rank->sum_of_ranks_within( data => $data );
    my $c     = 0;
    return sum0(
        map  { ++$c * $ranks->{$_} }
        sort { $a <=> $b } keys %{$ranks}
    );
}

sub _calc_l2r_value {
    my $data = shift;
    my $l    = _calc_l_value($data);
    my $n_bt = scalar keys %{$data};
    my $n_wt = __PACKAGE__->equal_n( data => $data );
    return ( ( ( 12 * $l ) / ( $n_wt * $n_bt * ( $n_bt**2 - 1 ) ) ) -
          ( ( 3 * ( $n_bt + 1 ) ) / ( $n_bt - 1 ) ) );
}

sub _calc_l_exp {
    my $data = shift;
    my $n_bt = scalar keys %{$data};
    my $n_wt = __PACKAGE__->equal_n( data => $data );
    return ( $n_wt * $n_bt * ( $n_bt + 1 )**2 ) / 4;
}

sub _calc_l_var {
    my $data = shift;
    my $n_bt = scalar keys %{$data};
    my $n_wt = __PACKAGE__->equal_n( data => $data );
    return ( $n_wt * $n_bt**2 * ( $n_bt + 1 ) * ( $n_bt**2 - 1 ) ) / 144;
}

sub _get_data {
    my ( $self, %args ) = @_;
    my $hoa;
    if ( ref $args{'data'} ) {
        $hoa = delete $args{'data'};
    }
    else {
        $hoa = $self->get_hoa_by_lab(%args);
    }
    return $hoa;
}

=head1 REFERENCES

Hollander, M., & Wolfe, D. A. (1999). I. New York, NY, US: Wiley.

Page, E. B. (1963). Ordered hypotheses for multiple treatments: A significance test for linear ranks. I, I<58>, 216-230. doi: L<10.1080/01621459.1963.10500843|http://dx.doi.org/10.1080/01621459.1963.10500843>. [L]

=head1 DEPENDENCIES

L : provides the handy sum0() function

L : used for probability functions.

L : used as a C for caching and retrieving data.

L : used to implement cross-case ranking.

L : for z-testing with optional continuity correction and tailing.

=head1 AUTHOR

Roderick Garton, C<<  >>

=head1 BUGS

Please report any bugs or feature requests to C, or through
the web interface at L.  I will be notified, and then you'll
automatically be notified of progress on your bug as I make changes.

=head1 SUPPORT

You can find documentation for this module with the perldoc command.

    perldoc Statistics::ANOVA::Page


You can also look for information at:

=over 4

=item * RT: CPAN's request tracker (report bugs here)

L

=item * AnnoCPAN: Annotated CPAN documentation

L

=item * CPAN Ratings

L

=item * Search CPAN

L

=back


=head1 ACKNOWLEDGEMENTS


=head1 LICENSE AND COPYRIGHT

Copyright 2015 Roderick Garton.

This program is free software; you can redistribute it and/or modify it
under the terms of either: the GNU General Public License as published
by the Free Software Foundation; or the Artistic License.

See L for more information.


=cut

1;    # End of Statistics::ANOVA::Page