Skip to content

Commit

Permalink
Add more tests for Chao2 calcs.
Browse files Browse the repository at this point in the history 
We check for cases (Q1==0 && Q2>0), (Q1 > 0 && Q2==0) and (Q1==0 && Q2==0).

Flush out an error in the confidence interval logic by doing so.

Updates issue #420
  • Loading branch information
@shawnlaffan
shawnlaffan committed Apr 25, 2016
1 parent 9a62216 commit 0e2fb57
Show file tree
Hide file tree
Showing 3 changed files with 374 additions and 223 deletions.
8 changes: 4 additions & 4 deletions lib/Biodiverse/Indices/RichnessEstimation.pm
View file
Expand Up @@ -258,7 +258,7 @@ sub calc_chao2 {
}
elsif ($Q1) {
$variance_uses_eq12 = 1;
$chao_formula = undef;
$chao_formula = 0;
}
# if only one singleton and no doubletons then chao stays zero
}
Expand All @@ -270,7 +270,7 @@ sub calc_chao2 {
if ($variance_uses_eq11) {
$variance = $c1 * ($Q1 * ($Q1 - 1)) / 2
+ $c2 * ($Q1 * (2 * $Q1 - 1) ** 2) / 4
- $c2 * $Q1 ** 4 / (4 * $chao);
- $c2 * $Q1 ** 4 / 4 / $chao;
}
elsif ($variance_uses_eq12) { # same structure as eq8 - could refactor
my %sums;
Expand All @@ -283,7 +283,7 @@ sub calc_chao2 {
$part2 += $i * exp (-$i) * $Q;
}
$variance = $part1 - $part2 ** 2 / $R;
$chao_formula = undef;
$chao_formula = 0;
}

$variance = max (0, $variance);
Expand Down Expand Up @@ -334,7 +334,7 @@ sub _calc_chao_confidence_intervals {

# and now the confidence interval
my ($lower, $upper);
if ($f1 && $f2) {
if (($f1 && $f2) || ($f1 > 1)) {
my $T = $chao - $richness;
my $K;
eval {
Expand Down
157 changes: 152 additions & 5 deletions t/24-Indices_richness_estimators.t
View file
Expand Up @@ -14,6 +14,10 @@ use Biodiverse::TestHelpers qw{
};
use Data::Section::Simple qw(get_data_section);

use Devel::Symdump;
my $obj = Devel::Symdump->rnew(__PACKAGE__);
my @subs = grep {$_ =~ 'main::test_'} $obj->functions();

exit main( @ARGV );

sub main {
Expand All @@ -31,16 +35,19 @@ sub main {
return 0;
}

test_indices($bd);
test_indices_1col($bd);


foreach my $sub (sort @subs) {
no strict 'refs';
$sub->($bd);
}

done_testing;
return 0;
}


sub test_indices {
my $bd = shift;
my $bd = shift->clone;

my $focal_gp = 'Broad_Meadow_Brook';
my @groups = grep {$_ ne $focal_gp} $bd->get_groups;
Expand All @@ -64,7 +71,7 @@ sub test_indices {
}

sub test_indices_1col {
my $bd = shift;
my $bd = shift->clone;

my $focal_gp = 'Broad_Meadow_Brook';
my @groups = grep {$_ ne $focal_gp} $bd->get_groups;
Expand Down Expand Up @@ -124,6 +131,146 @@ sub test_indices_1col {
}


# chao2 differs when q1 or q2 are zero
sub test_chao2_Q2_no_Q1 {
my $bd = shift->clone;

# need to ensure there are no uniques - make them all occur everywhere
foreach my $label ($bd->get_labels) {
if ($bd->get_range(element => $label) == 1) {
foreach my $group ($bd->get_groups) {
$bd->add_element (group => $group, label => $label);
}
}
}

my $focal_gp = 'Broad_Meadow_Brook';
my @nbr_groups = grep {$_ ne $focal_gp} $bd->get_groups;

my $results2 = {
CHAO2 => 26, CHAO2_SE => 0.629523,
CHAO2_Q1_COUNT => 0, CHAO2_Q2_COUNT => 8,
CHAO2_UNDETECTED => 0, CHAO2_VARIANCE => 0.396299,
CHAO2_CI_LOWER => 26.030028, CHAO2_CI_UPPER => 28.623692,
CHAO2_META => {
CHAO_FORMULA => 0,
CI_FORMULA => 14,
VARIANCE_FORMULA => 12,
},
};

my %expected_results = (2 => $results2);

run_indices_test1 (
calcs_to_test => [qw/
calc_chao2
/],
sort_array_lists => 1,
basedata_ref => $bd,
element_list1 => ['Broad_Meadow_Brook'],
element_list2 => \@nbr_groups,
expected_results => \%expected_results,
skip_nbr_counts => {1 => 1},
#generate_result_sets => 1,
descr_suffix => 'test_chao2_Q2_no_Q1',
);

return;
}

sub test_chao2_Q1_no_Q2 {
my $bd = shift->clone;

# need to ensure there are no uniques - make them all occur everywhere
foreach my $label ($bd->get_labels) {
if ($bd->get_range(element => $label) == 2) {
foreach my $group ($bd->get_groups) {
$bd->add_element (group => $group, label => $label);
}
}
}

my $focal_gp = 'Broad_Meadow_Brook';
my @nbr_groups = grep {$_ ne $focal_gp} $bd->get_groups;

my $results2 = {
CHAO2 => 39.636364, CHAO2_SE => 12.525204,
CHAO2_Q1_COUNT => 6, CHAO2_Q2_COUNT => 0,
CHAO2_UNDETECTED => 13.636364, CHAO2_VARIANCE => 156.880734,
CHAO2_CI_LOWER => 28.943875, CHAO2_CI_UPPER => 89.165183,
CHAO2_META => {
CHAO_FORMULA => 4,
CI_FORMULA => 13,
VARIANCE_FORMULA => 11,
},
};

my %expected_results = (2 => $results2);

run_indices_test1 (
calcs_to_test => [qw/
calc_chao2
/],
sort_array_lists => 1,
basedata_ref => $bd,
element_list1 => ['Broad_Meadow_Brook'],
element_list2 => \@nbr_groups,
expected_results => \%expected_results,
skip_nbr_counts => {1 => 1},
#generate_result_sets => 1,
descr_suffix => 'test_chao2_Q1_no_Q2',
);

return;
}

sub test_chao2_no_Q1_no_Q2 {
my $bd = shift->clone;

# need to ensure there are no uniques or doubles - make them all occur everywhere
foreach my $label ($bd->get_labels) {
my $range = $bd->get_range(element => $label);
if ($range == 1 || $range == 2) {
foreach my $group ($bd->get_groups) {
$bd->add_element (group => $group, label => $label);
}
}
}

my $focal_gp = 'Broad_Meadow_Brook';
my @nbr_groups = grep {$_ ne $focal_gp} $bd->get_groups;

my $results2 = {
CHAO2 => 26, CHAO2_SE => 0.3696727,
CHAO2_Q1_COUNT => 0, CHAO2_Q2_COUNT => 0,
CHAO2_UNDETECTED => 0, CHAO2_VARIANCE => 0.1366579,
CHAO2_CI_LOWER => 26, CHAO2_CI_UPPER => 26.910745,
CHAO2_META => {
CHAO_FORMULA => 0,
CI_FORMULA => 14,
VARIANCE_FORMULA => 12,
},
};

my %expected_results = (2 => $results2);

run_indices_test1 (
calcs_to_test => [qw/
calc_chao2
/],
sort_array_lists => 1,
basedata_ref => $bd,
element_list1 => ['Broad_Meadow_Brook'],
element_list2 => \@nbr_groups,
expected_results => \%expected_results,
skip_nbr_counts => {1 => 1},
#generate_result_sets => 1,
descr_suffix => 'test_chao2_no_Q1_no_Q2',
);

return;
}

sub get_basedata {

my $data = get_data_section ('SAMPLE_DATA');
Expand Down

0 comments on commit 0e2fb57

Please sign in to comment.