File Coverage

blib/lib/Statistics/ANOVA/KW.pm

Criterion	Covered	Total	%
statement	68	78	87.1
branch	15	28	53.5
condition	1	3	33.3
subroutine	14	16	87.5
pod	5	5	100.0
total	103	130	79.2

package Statistics::ANOVA::KW;

use 5.006;
use strict;
use warnings FATAL => 'all';
use base qw(Statistics::Data);
use Carp qw(croak);
use List::AllUtils qw(sum0);
use Math::Cephes qw(:dists);
use Statistics::Data::Rank;
use Statistics::Lite qw(mean);
$Statistics::ANOVA::KW::VERSION = '0.01';

=head1 NAME

Statistics::ANOVA::KW - Kruskall-Wallis statistics and test (nonparametric independent analysis of variance by ranks for nominally grouped data)

=head1 VERSION

This is documentation for B of Statistics::ANOVA::KW.

=head1 SYNOPSIS

 use Statistics::ANOVA::KW;
 my $kw = Statistics::ANOVA::KW->new();
 $kw->load({1 => [2, 4, 6], 2 => [3, 3, 12], 3 => [5, 7, 11, 16]});
 my $h_value = $kw->h_value(); # default used to correct for ties
 my $p_value = $kw->chiprob_test(); # H taken as chi^2 distributed
 my ($h_value, $df, $count, $p_value_by_chi, $phi) = $kw->chiprob_test(); # same as above, called in array context
 my ($f_value, $df_b, $df_w, $p_value_by_f, $omega_sq) = $kw->fprob_test(); # F-equivalent value tests

 # or without pre-loading, and specify correct_ties as well:
 $h_value = $kw->h_value(data => {1 => [2, 4, 6], 2 => [5, 3, 12]}, correct_ties => 1);
 # or test only a subset of the loaded data:
 $h_value = $kw->h_value(lab => [1, 3]);

=head1 DESCRIPTION

Performs calculations for the Kruskal-Wallis one-way nonparametric analysis of variance by ranks. This is for (at least) ordinal-level measurements of two or more samples of a nominal/categorical variable with equality of variances across the samples. The test is unreliable for small number of observations per sample (conventionally, all samples should have more than five observations). See L for more information, and discussions of the assumptions/interpretations, and pros/cons, of the test at L (pro) and L (con). Note that the Kruskall-Wallis test is often described as a test for I or more samples, in contrast to the Mann-Whitney test, which is restricted to two samples, but KW can also be used with only two samples: the absolute value of the I-value from a Mann-Whitney test equals the square-root of the KW statistic for two factors.

Data-loading and retrieval are as provided in L, on which this module's class object is Cd, so its other methods are available here.

Return values are tested on installation against published examples and output from other software (e.g., SPSS).

=head2 new

 $kw = Statistics::ANOVA::KW->new();

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

=head2 load, add, unload

 $kw->load('a' => [1, 4, 3.2], 'b' => [6.5, 6.5, 9], 'c' => [3, 7, 4.4]);

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. Once loaded, subsets of the loaded data can be tested by passing their names (or labels) in a referenced array to the argument B in the following methods (as supported by L). Once loaded, any non-numeric values in the samples are culled ahead of running the following methods.

Alternatively, without pre-loading the data, directly give the following methods the HOA of data as the value for the optional named argument B.

=cut

=head2 h_value

 $h_value = $kw->h_value(data => \%data, correct_ties => 1);
 $h_value = $kw->h_value(); # assuming data have already been loaded, & default of TRUE for correct_ties

Returns the Kruskall-Wallis I statistic.

=cut

sub h_value {
    my ( $self, %args ) = ( shift, @_ );
    my $data =
      $args{'data'}
      ? delete $args{'data'}
      : $self->get_hoa_by_lab_numonly_indep(%args);
    my $correct_ties = defined $args{'correct_ties'}
      and $args{'correct_ties'} == 0 ? 0 : 1;
    return ( _kw_stats( $data, $correct_ties ) )[0];
}

=head2 chiprob_test

 ($chi_value, $df, $count, $p_value, $phi) = $kw->chiprob_test(data => HOA, correct_ties => 1); # H as chi-square
 $p_value = $kw->chiprob_test(data => HOA, correct_ties => 1);
 $p_value = $kw->chiprob_test(); # assuming data have already been loaded, & default of TRUE for correct_ties

Performs the ANOVA and, assuming I-square distribution of the Kruskall-Wallis I value, returns its value, its degrees-of-freedom, the total number of observations (I), its I-square probability value, and I-coefficient as an estimate of effect-size ( = square-root of  (I-square divided by I) ). Returns only the I-value if called in scalar context. Default value of optional argument B is 1.

=cut

sub chiprob_test {
    my ( $self, %args ) = ( shift, @_ );
    my $data =
      $args{'data'}
      ? delete $args{'data'}
      : $self->get_hoa_by_lab_numonly_indep(%args);
    my $correct_ties = defined $args{'correct_ties'}
      and $args{'correct_ties'} == 0 ? 0 : 1;
    my ( $chi, $df_b, $count ) = _kw_stats( $data, $correct_ties );
    my $p_value = chdtrc( $df_b, $chi );    # Math::Cephes fn
    return
      wantarray
      ? ( $chi, $df_b, $count, $p_value, sqrt( $chi / $count ) )
      : $p_value;
}

=head2 chiprob_str

 $str = $kw->chiprob_str(data => HOA, correct_ties => 1);
 $str = $kw->chiprob_str(); # assuming data have already been loaded, & default of TRUE for correct_ties

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, $phi ) = $self->chiprob_test(%args);
    return "chi^2($df, N = $count) = $chi_value, p = $p_value, phi = $phi";
}

=head2 fprob_test

 ($f_value, $df_b, $df_w, $p_value, $es_omega) = $kw->fprob_test(data => HOA, correct_ties => BOOL);
 $p_value = $kw->fprob_test(data => HOA, correct_ties => BOOL);
 $p_value = $kw->fprob_test(); # assuming data have already been loaded, & default of TRUE for correct_ties

Performs the same test as above but transforms the I-square value into an I-distributed value, returning an array comprised of (1) this I-estimate value, its (2) between- and (3) within-groups degrees-of-freedom, (4) the associated probability of the value per the I-distribution, and (5) an estimate of the effect-size statistic, (partial) I-squared. The latter is returned only if L is installed and available. Called in scalar context, only the I-estimated I-value is returned. The default value of the optional argument B is 1. This method has not been tested against sample/published data (not being provided in the usual software packages).

=cut

sub fprob_test {
    my ( $self, %args ) = ( shift, @_ );
    my $data =
      $args{'data'}
      ? delete $args{'data'}
      : $self->get_hoa_by_lab_numonly_indep(%args);
    my $correct_ties = defined $args{'correct_ties'}
      and $args{'correct_ties'} == 0 ? 0 : 1;
    my ( $f_value, $df_b, $df_w ) = _f_stats( $data, $correct_ties );
    my $p_value = fdtrc( $df_b, $df_w, $f_value );    # Math::Cephes fn
    return $p_value if !wantarray;

    my $es_omega;
    eval { require Statistics::ANOVA::EffectSize; };
    if ( !$@ ) {
        $es_omega = Statistics::ANOVA::EffectSize->omega_sq_partial_by_f(
            f_value => $f_value,
            df_b    => $df_b,
            df_w    => $df_w
        );
    }
    return ( $f_value, $df_b, $df_w, $p_value, $es_omega );
}

=head2 fprob_str

 $str = $kw->chiprob_str(data => HOA, correct_ties => BOOL);
 $str = $kw->chiprob_str(); # assuming data have already been loaded, using default of TRUE for correct_ties

Performs the same test as for L but returns not an array but a string of the conventional reporting form, e.g., F(df_b, df_w) = f_value, p = p_value (and also, then, an estimate of partial I-squared, if available, see above).

=cut

sub fprob_str {
    my ( $self, %args ) = ( shift, @_ );
    my ( $f_value, $df_b, $df_w, $p_value, $es_omega ) =
      $self->fprob_test(%args);
    my $str = "F($df_b, $df_w) = $f_value, p = $p_value";
    if ( defined $es_omega ) {
        $str .= ', est_omega^2_p = ' . $es_omega;
    }
    return $str;
}

sub _kw_stats {
    my ( $data, $correct_ties ) = @_;
    my ( $ranks_href, $ties_aref, $gn, $ties_var ) =
      Statistics::Data::Rank->ranks_between( data => $data );
    my $num = sum0(
        map {
            scalar @{ $ranks_href->{$_} } *
              ( mean( @{ $ranks_href->{$_} } ) - ( ( $gn + 1 ) / 2 ) )**2
        } keys %{$ranks_href}
    );
    my $h = 12 / ( $gn * ( $gn + 1 ) ) * $num;

    # correction for ties:
    $h /= ( 1 - ( $ties_var / ( $gn**3 - $gn ) ) )
      unless defined $correct_ties and not $correct_ties;
    return ( $h, ( scalar keys %{$ranks_href} ) - 1, $gn ); # H, df, and grand N
}

sub _f_stats {
    my ( $data, $correct_ties ) = @_;
    my ( $h, $df_b, $n ) = _kw_stats( $data, $correct_ties );
    my $df_w  = sum0( map { scalar @{ $data->{$_} } - 1 } keys %{$data} );
    my $n_bt  = scalar keys( %{$data} );
    my $f_val = ( $h / ( $n_bt - 1 ) ) / ( ( $n - 1 - $h ) / ( $n - $n_bt ) );
    return ( $f_val, $df_b, $df_w );
}

=head1 DEPENDENCIES

L : used for summing.

L : used for probability functions.

L : used as base.

L : used to calculate between-group sum-of ranks.

=head1 REFERENCES

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

Rice, J. A. (1995). I. Belmont, CA, US: Duxbury.

Sarantakos, S. (1993). I. Melbourne, Australia: MacMillan.

Siegal, S. (1956). I. New York, NY, US: McGraw-Hill

=head1 SEE ALSO

L : Also a nonparametric ANOVA by ranks for independent samples, but where the ordinality of the numerical labels of the sample names (the order of the groups) is taken into account.

L : Returns the I-value and its I-square I-value (only), and implements the Newman-Keuls test for pairwise comparison.

=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::KW

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::KW