P-value correction: Difference between revisions

Content added Content deleted
(added Ruby)
(→‎{{header|Perl}}: removed unused variable and unnecessary checks)
Line 3,607: Line 3,607:


use strict;
use strict;
use warnings;
use warnings FATAL => 'all';
use autodie ':all';
⚫
use List::Util 'min';
use feature 'say';


sub pmin {
sub pmin {
my $array_ref = shift;
my $array = shift;
my $x = 1;
my $x = 1;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @pmin_array;
my @pmin_array;
my $n = scalar @$array_ref;
my $n = scalar @$array;
for (my $index = 0; $index < $n; $index++) {
for (my $index = 0; $index < $n; $index++) {
if (@$array_ref[$index] < $x) {
$pmin_array[$index] = min(@$array[$index], $x);
$pmin_array[$index] = @$array_ref[$index];
} else {
$pmin_array[$index] = $x;
}
}
}
return @pmin_array;
@pmin_array
}
}


sub cummin {
sub cummin {
my $array_ref = shift;
my $array_ref = shift;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @cummin;
my @cummin;
my $cumulative_min = @$array_ref[0];
my $cumulative_min = @$array_ref[0];
Line 3,642: Line 3,633:
push @cummin, $cumulative_min;
push @cummin, $cumulative_min;
}
}
return @cummin;
@cummin
}
}


sub cummax {
sub cummax {
my $array_ref = shift;
my $array_ref = shift;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @cummax;
my @cummax;
my $cumulative_max = @$array_ref[0];
my $cumulative_max = @$array_ref[0];
Line 3,659: Line 3,646:
push @cummax, $cumulative_max;
push @cummax, $cumulative_max;
}
}
return @cummax;
@cummax
}
}


Line 3,675: Line 3,662:
die;
die;
}
}
}
unless ((ref $array_ref) =~ m/ARRAY/) {
print "You should have entered an array.\n";
die;
}
}
my @array;
my @array;
Line 3,687: Line 3,670:
@array = sort { @$array_ref[$b] <=> @$array_ref[$a] } 0..$max_index;
@array = sort { @$array_ref[$b] <=> @$array_ref[$a] } 0..$max_index;
}
}
@array
}
}


⚫
use List::Util 'min';


sub p_adjust {
sub p_adjust {
my $pvalues_ref = shift;
my $pvalues_ref = shift;
unless ((ref $pvalues_ref) =~ m/ARRAY/) {
print "p_adjust requires an array.\n";
die;
}
my $method;
my $method;
if (defined $_[0]) {
if (defined $_[0]) {
$method = shift;
$method = shift
} else {
} else {
$method = 'Holm';
$method = 'Holm'
}
}
my %methods = (
my %methods = (
Line 3,717: Line 3,696:
$method = $key;
$method = $key;
$method_found = 'yes';
$method_found = 'yes';
last;
last
}
}
}
}
Line 3,731: Line 3,710:
if ($method_found eq 'no') {
if ($method_found eq 'no') {
print "No method could be determined from $method.\n";
print "No method could be determined from $method.\n";
die;
die
}
}
my $lp = scalar @$pvalues_ref;
my $lp = scalar @$pvalues_ref;
Line 3,743: Line 3,722:
}
}
my @cummin = cummin(\@cummin_input);
my @cummin = cummin(\@cummin_input);
⚫
undef @cummin_input;
my @pmin = pmin(\@cummin);
my @pmin = pmin(\@cummin);
undef @cummin;
my @ro = order(\@o);
my @ro = order(\@o);
undef @o;
@qvalues = @pmin[@ro];
@qvalues = @pmin[@ro];
} elsif ($method eq 'bh') {
} elsif ($method eq 'bh') {
Line 3,756: Line 3,732:
}
}
my @ro = order(\@o);
my @ro = order(\@o);
undef @o;
my @cummin = cummin(\@cummin_input);
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
my @pmin = pmin(\@cummin);
my @pmin = pmin(\@cummin);
undef @cummin;
@qvalues = @pmin[@ro];
@qvalues = @pmin[@ro];
} elsif ($method eq 'by') {
} elsif ($method eq 'by') {
Line 3,773: Line 3,746:
$cummin_input[$index] = $q * ($n/($n-$index)) * @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
$cummin_input[$index] = $q * ($n/($n-$index)) * @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
}
}
⚫
# say join (',', @cummin_input);
# say '@cummin_input # of elements = ' . scalar @cummin_input;
my @cummin = cummin(\@cummin_input);
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
undef @cummin_input;
my @pmin = pmin(\@cummin);
my @pmin = pmin(\@cummin);
undef @cummin;
@qvalues = @pmin[@ro];
@qvalues = @pmin[@ro];
} elsif ($method eq 'bonferroni') {
} elsif ($method eq 'bonferroni') {
Line 3,786: Line 3,760:
$qvalues[$index] = 1.0;
$qvalues[$index] = 1.0;
} else {
} else {
print "Failed to get Bonferroni adjusted p.";
say 'Failed to get Bonferroni adjusted p.';
die;
die;
}
}
Line 3,804: Line 3,778:
@qvalues = @pmin[@ro];
@qvalues = @pmin[@ro];
} elsif ($method eq 'hommel') {
} elsif ($method eq 'hommel') {
my @i = 1..$n;
my @o = order($pvalues_ref);
my @o = order($pvalues_ref);
my @p = @$pvalues_ref[@o];
my @p = @$pvalues_ref[@o];
my @ro = order(\@o);
my @ro = order(\@o);
undef @o;
undef @o;
my @pa;
my (@q, @pa);
my @q;
my $min = $n*$p[0];
my $min = $n*$p[0];
for (my $index = 0; $index < $n; $index++) {
for (my $index = 0; $index < $n; $index++) {
my $temp = $n*$p[$index] / ($index + 1);
my $temp = $n*$p[$index] / ($index + 1);
if ($temp < $min) {
$min = min($min, $temp);
$min = $temp;
}
}
}
for (my $index = 0; $index < $n; $index++) {
for (my $index = 0; $index < $n; $index++) {
Line 3,823: Line 3,793:
}
}
for (my $j = ($n-1); $j >= 2; $j--) {
for (my $j = ($n-1); $j >= 2; $j--) {
⚫
my @ij = 0..($n - $j);#ij <- seq_len(n - j + 1)
# printf("j = %zu\n", j);
⚫
my @ij = 1..($n - $j + 1);#ij <- seq_len(n - j + 1)
for (my $i = 0; $i < $n - $j + 1; $i++) {
$ij[$i]--;#R's indices are 1-based, C's are 0-based
}
my $I2_LENGTH = $j - 1;
my $I2_LENGTH = $j - 1;
my @i2;
my @i2;
Line 3,838: Line 3,804:
for (my $i = 1; $i < $I2_LENGTH; $i++) {#loop through 2:j
for (my $i = 1; $i < $I2_LENGTH; $i++) {#loop through 2:j
my $TEMP_Q1 = $j * $p[$i2[$i]] / (2 + $i);
my $TEMP_Q1 = $j * $p[$i2[$i]] / (2 + $i);
if ($TEMP_Q1 < $q1) {
$q1 = min($TEMP_Q1, $q1);
$q1 = $TEMP_Q1;
}
}
}

for (my $i = 0; $i < ($n - $j + 1); $i++) {#q[ij] <- pmin(j * p[ij], q1)
for (my $i = 0; $i < ($n - $j + 1); $i++) {#q[ij] <- pmin(j * p[ij], q1)
$q[$ij[$i]] = min( $j*$p[$ij[$i]], $q1);
$q[$ij[$i]] = min( $j*$p[$ij[$i]], $q1);
Line 3,848: Line 3,811:


for (my $i = 0; $i < $I2_LENGTH; $i++) {#q[i2] <- q[n - j + 1]
for (my $i = 0; $i < $I2_LENGTH; $i++) {#q[i2] <- q[n - j + 1]
$q[$i2[$i]] = $q[$n - $j];#subtract 1 because of starting index difference
$q[$i2[$i]] = $q[$n - $j];
}
}


Line 3,862: Line 3,825:
} else {
} else {
print "$method doesn't fit my types.\n";
print "$method doesn't fit my types.\n";
die;
die
}
}
return @qvalues;
@qvalues
}
}
my @pvalues = (4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
my @pvalues = (4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
Retrieved from "https://rosettacode.org/wiki/P-value_correction"