File Coverage

blib/lib/Statistics/Multtest.pm

Criterion	Covered	Total	%
statement	186	197	94.4
branch	26	38	68.4
condition	12	24	50.0
subroutine	32	32	100.0
pod	0	6	0.0
total	256	297	86.2

package Statistics::Multtest;

use List::Vectorize;
use Carp;
use strict;
use constant {EPS => 1e-10};

require Exporter;
our @ISA = qw(Exporter);
our @EXPORT_OK = qw(bonferroni holm hommel hochberg BH BY qvalue);
our %EXPORT_TAGS = (all => [qw(bonferroni holm hommel hochberg BH BY qvalue)]);

our $VERSION = '0.13';

1;

BEGIN {
    no strict 'refs';
    for my $adjp (qw(bonferroni holm hommel hochberg BH BY qvalue)) {
        *{$adjp} = sub {
            my $p = shift;
            my $name;
            ($name, $p) = initial($p);
            *private = *{"_$adjp"};
            my $adjp = private($p, @_);
            return get_result($adjp, $name);
        }
    }
}

sub initial {
    my $p = shift;

    if(! List::Vectorize::is_array_ref($p) and ! List::Vectorize::is_hash_ref($p)) {
        croak "ERROR: P-values should be stored in an array reference or a hash reference.\n";
    }

    my $name = [];
    if(List::Vectorize::is_hash_ref($p)) {
        $name = [ keys %{$p} ];
        $p = [ values %{$p} ];
    }

    if(max($p) > 1 or min($p) < 0) {
        croak "ERROR: P-values should between 0 and 1.\n";
    }

    return ($name, $p);
}

sub get_result {
    my ($adjp, $name) = @_;

    if(is_empty($name)) {
        return $adjp;
    }
    else {
        my $result;
        for (0..$#$name) {
            $result->{$name->[$_]} = $adjp->[$_];
        }
        return $result;
    }
}

# R code: pmin(1, n * p)
sub _bonferroni {
    my $p = shift;
    my $n = len($p);

    my $adjp = multiply($n, $p);

    return pmin(1, $adjp);
}

# R code: i = 1:n
#         o = order(p)
#         ro = order(o)
#         pmin(1, cummax((n - i + 1) * p[o]))[ro]
sub _holm {
    my $p = shift;
    my $n = len($p);

    my $i = seq(1, $n);
    my $o = order($p);
    my $ro = order($o);

    my $adjp = multiply(minus($n + 1, $i), subset($p, $o));
    $adjp = cumf($adjp, \&max);
    $adjp = pmin(1, $adjp);

    return subset($adjp, $ro);
}

# R code: i = 1:n
#         o = order(p)
#         p = p[o]
#         ro = order[o]
#         q = pa = rep.int(min(n * p/i), n)
#         for (j in (n - 1):2) {
#             ij = 1:(n - j + 1)
#             i2 = (n - j + 2):n
#             q1 <- min(j * p[i2]/(2:j))
#             q[ij] <- pmin(j * p[ij], q1)
#             q[i2] <- q[n - j + 1]
#             pa <- pmax(pa, q)
#         }
#         pmax(pa, p)[ro]
sub _hommel {
    my $p = shift;
    my $n = len($p);

    my $i = seq(1, $n);
    my $o = order($p);
    $p = subset($p, $o);
    my $ro = order($o);

    my $pa = rep(min(divide(multiply($n, $p), $i)), $n);
    my $q = copy($pa);

    # set the first index as 1
    unshift(@$p, 0);
    unshift(@$q, 0);
    unshift(@$pa, 0);

    my $ij;
    my $i2;
    my $q1;
    for my $j (@{seq($n - 1, 2)}) {

        $ij = seq(1, $n - $j + 1);
        $i2 = seq($n - $j + 2, $n);
        $q1 = min(divide(multiply($j, subset($p, $i2)), seq(2, $j)));
        subset_value($q, $ij, pmin(multiply($j, subset($p, $ij)), $q1));
        subset_value($q, $i2, $q->[$n - $j + 1]);
        $pa = pmax($pa, $q);
    }

    shift(@$p);
    shift(@$q);
    shift(@$pa);

    my $adjp = pmax($pa, $p);
    return subset($adjp, $ro);
}

# R code: i = n:1
#         o <- order(p, decreasing = TRUE)
#         ro <- order(o)
#         pmin(1, cummin((n - i + 1) * p[o]))[ro]
sub _hochberg {

    my $p = shift;
    my $n = len($p);
    my $i = seq($n, 1);

    my $o = order($p, sub {$_[1] <=> $_[0]});
    my $ro = order($o);

    my $adjp = multiply(minus($n+1, $i), subset($p, $o));
    $adjp = cumf($adjp, \&min);
    $adjp = pmin(1, $adjp);
    return subset($adjp, $ro);
}

# R code: i <- n:1
#         o <- order(p, decreasing = TRUE)
#         ro <- order(o)
#         pmin(1, cummin(n/i * p[o]))[ro]
sub _BH {
    my $p = shift;
    my $n = len($p);
    my $i = seq($n, 1);

    my $o = order($p, sub {$_[1] <=> $_[0]});
    my $ro = order($o);

    my $adjp = multiply(divide($n, $i), subset($p, $o));
    $adjp = cumf($adjp, \&min);
    $adjp = pmin(1, $adjp);
    return subset($adjp, $ro);

}

# R code: i <- n:1
#         o <- order(p, decreasing = TRUE)
#         ro <- order(o)
#         q <- sum(1/(1L:n))
#         pmin(1, cummin(q * n/i * p[o]))[ro]
sub _BY {

    my $p = shift;
    my $n = len($p);
    my $i = seq($n, 1);

    my $o = order($p, sub {$_[1] <=> $_[0]});
    my $ro = order($o);

    my $q = sum(divide(1, seq(1, $n)));
    my $adjp = multiply(divide($q*$n, $i), subset($p, $o));
    $adjp = cumf($adjp, \&min);
    $adjp = pmin(1, $adjp);
    return subset($adjp, $ro);
}

sub pmin {
    my $array1 = shift;
    my $array2 = shift;

    return mapply($array1, $array2, sub {min(\@_)});
}

sub pmax {
    my $array1 = shift;
    my $array2 = shift;

    return mapply($array1, $array2, sub {max(\@_)});
}


# the default setting can not assure the success of calculation
# so this function should run under eval
# eval(q(qvalue $p));
# if ($@) change_other_settings
# @p = qw(0.76 0.30 0.40 0.43 0.27 0.79 0.66 0.36 0.12 0.16 0.52 0.04 0.67 0.44 0.40 0.09 0.51 0.38 0.49 0.68)
sub _qvalue {
    my $p = shift;
    my %setup = ('lambda'         => multiply(seq(0, 90, 5), 0.01),
                 'robust'         => 0,
                 @_ );

    my $lambda = $setup{'lambda'};
    my $robust = $setup{'robust'};

    if(! List::Vectorize::is_array_ref($lambda)) {
        croak "ERROR: lambda should be an array reference. If your lambda is a single number, you should set it as an array reference containing one element, such as [0.1].";
    }

    if(len($lambda) > 1 and len($lambda) < 4) {
        croak "ERROR: If length of lambda greater than 1, you need at least 4 values.";
    }
    if(len($lambda) > 1 and (min($lambda) < 0 or max($lambda) >= 1)) {
        croak "ERROR: Lambda must be within [0, 1).";
    }

    # m <- length(p)
    my $m = len($p);
    my $pi0;

    if(len($lambda) == 1) {
        $lambda = $lambda->[0];
        if($lambda < 0 or $lambda >= 1) {
            croak "ERROR: Lambda must be within [0, 1).";
        }
        # pi0 <- mean(p >= lambda)/(1-lambda)
        $pi0 = mean(test($p, sub {_approximately_compare($_[0], ">=", $lambda)})) / (1 - $lambda);
        # pi0 <- min(pi0,1)
        $pi0 = min(c($pi0, 1));
    }
    else {
        # pi0 <- rep(0,length(lambda))
        $pi0 = rep(0, len($lambda));
        # for(i in 1:length(lambda)) {
        #     pi0[i] <- mean(p >= lambda[i])/(1-lambda[i])
        # }
        $pi0 = sapply(seq(0, len($lambda) - 1),
                      sub {
                          my $current_lambda = $lambda->[$_[0]];
                          mean(test($p, sub {_approximately_compare($_[0], ">=", $current_lambda)})) / (1 - $current_lambda);
                      });

        # minpi0 <- min(pi0)
        my $minpi0 = min($pi0);

        # mse <- rep(0,length(lambda))
        my $mse = rep(0, len($lambda));
        # pi0.boot <- rep(0,length(lambda))
        my $pi0_boot = rep(0, len($lambda));
        for (1..100) {
            #  p.boot <- sample(p,size=m,replace=TRUE)
            my $p_boot = sample($p, $m, 'replace' => 1);
            # for(i in 1:length(lambda)) {
            #     pi0.boot[i] <- mean(p.boot>lambda[i])/(1-lambda[i])
            # }
            $pi0_boot = sapply(seq(0, len($lambda) - 1),
                               sub {
                                    my $current_lambda = $lambda->[$_[0]];
                                    mean(test($p_boot, sub {$_[0] > $current_lambda})) / (1 - $current_lambda);
                               });
            # mse <- mse + (pi0.boot-minpi0)^2
            $mse = plus($mse, power(minus($pi0_boot, $minpi0), 2));
        }

        # pi0 <- min(pi0[mse==min(mse)])
        my $min_mse = min($mse);

        $pi0 = min(subset($pi0, which(test($mse, sub {_approximately_compare($_[0], "==", $min_mse)}))));
        $pi0 = min(c($pi0, 1));
    }

    if($pi0 <= 0) {
        croak "ERROR: The estimated pi0 ($pi0) <= 0. Check that you have valid p-values or use another lambda method.";
    }

    # u <- order(p)
    my $u = order($p);

    # this function is not the same to the function in qvalue
    # but returns the same result
    # returns number of observations less than or equal to
    sub qvalue_rank {
        my $x = shift;
        return sapply($x, sub {
                          my $current_x = $_[0];
                          sum(test($x, sub {$_[0] <= $current_x}));
                      });
    }

    # v <- qvalue.rank(p)
    my $v = qvalue_rank($p);
    # qvalue <- pi0*m*p/v
    my $qvalue = divide(multiply($pi0, $m, $p), $v);

    if($robust) {
        # qvalue <- pi0*m*p/(v*(1-(1-p)^m))
        $qvalue = divide(multiply($pi0, $m, $p), multiply($v, minus(1, power(minus(1, $p), $m))));
    }

    # qvalue[u[m]] <- min(qvalue[u[m]],1)
    $qvalue->[$u->[$#$u]] = min(c($qvalue->[$u->[$#$u]], 1));
    # for(i in (m-1):1) {
    #     qvalue[u[i]] <- min(qvalue[u[i]],qvalue[u[i+1]],1)
    # }
    for my $i (@{seq($#$u-1, 0)}) {
        $qvalue->[$u->[$i]] = min([$qvalue->[$u->[$i]], $qvalue->[$u->[$i+1]], 1]);
    }

    return $qvalue;
}

# vectorized power calculation
sub power {
    my $array = shift;
    my $power = shift;

    my $res = sapply($array, sub {$_[0]**($power)});
    return $res;
}

# approximately
sub _approximately_compare {
    my $left = shift;
    my $sign = shift;
    my $right = shift;

    if($sign eq ">" or $sign eq ">=") {
        if($left > $right) {
            return 1;
        }
        elsif(abs($left - $right) < EPS) {
            return 1;
        }
        else {
            return 0;
        }
    }
    elsif($sign eq "<" or $sign eq "<=") {
        if($left < $right) {
            return 1;
        }
        elsif(abs($left - $right) < EPS) {
            return 1;
        }
        else {
            return 0;
        }
    }
    elsif($sign eq "==") {
        if(abs($left - $right) < EPS) {
            return 1;
        }
        else {
            return 0;
        }
    }
}

__END__