Fixed $PRIOR support and bugfix in shrinkage results handling

[PsN.git] / bin / bootstrap

#!/usr/local/bin/perl

use FindBin qw($Bin);
use lib "$Bin/../lib";

# Don't edit the line below, it must look exactly like this.
# Everything above this line will be replaced #

use PsN;
use model;
use tool::bootstrap;
use strict;
use debug;
use ui;
use Getopt::Long;
use common_options;
use Data::Dumper;

my $cmd_line = $0 . " " . join( " ", @ARGV );

## Configure the command line parsing
Getopt::Long::config("auto_abbrev");

## Declare the options
my %options;

my %required_options = ();
my %optional_options = ("samples:i"=>'',
                        "sample_size:s"=>'',
                        "stratify_on:s"=>'',
                        "bca"=>'',
                        "mplots"=>'',
                        "rplots"=>'',
                        "skip_minimization_terminated!"=>'',
                        "skip_covariance_step_terminated!"=>'',
                        "skip_with_covstep_warnings!"=>'',
                        "skip_estimate_near_boundary!"=>'',
                        "minimization_successful_limit:f"=>'',
                        "covariance_step_successful_limit:f"=>'',
                        "covariance_step_warnings_limit:f"=>'',
                        "estimate_near_boundary_limit:f"=>'',
                        "estimate_near_boundary_limit:f"=>'',
                        "se_confidence_intervals_level:f"=>'',
                        "percentile_confidence_intervals_level:f"=>'',
                        "bca_confidence_intervals_level:f"=>'',
                        "se_confidence_intervals_check:f"=>'',
                        "percentile_confidence_intervals_check:f"=>'',
                        "bca_confidence_intervals_check:f"=>'',
                        "large_bias_limit:f"=>'' );

my $res = GetOptions( \%options, 
                      @common_options::get_opt_strings,
                      keys(%required_options),
                      keys(%optional_options) );

exit unless $res;

common_options::set_globals( \%options, 'bootstrap' );
common_options::get_defaults( \%options, 'bootstrap' );
common_options::sanity_checks( \%options, 'bootstrap' );

my %help_text;
$help_text{Pre_help_message} = <<'EOF';  
    <h3 class="heading1">bootstrap</h3>
      
      Perl script for non-parametric bootstrap of NONMEM runs.
      
    <h3 class="heading1">Usage:</h3>
EOF

    $help_text{Description} = <<'EOF';
   <h3 class="heading1">Description:</h3>

      The Bootstrap can be used to calculate bias, standard errors and
      confidence intervals. It does so by resampling with replacement
      from the data, see Efron B, An Introduction to the Bootstrap,
      Chap. & Hall, London UK, 1993. To compute standard errors for
      all parameters of a model using the non-parametric bootstrap
      implemented here, roughly 200 model fits are necessary. To assess
      95% confidence intervals approximatly 2000 runs will suffice.
EOF
    $help_text{Examples} = <<'EOF';
   <h3 class="heading1">Example:</h3>

      <p class="option">bootstrap -samples=200 run89.mod</p>

      This will run a non-parametric bootstrap of 200 samples and give
      you good estimates of the standard errors of the parameter
      estimates. You may get some estimates for the confidence
      intervals too, but they will generally not be of high quality.


      <p class="option">bootstrap -samples=2000 -bca run89.mod</p>

      This will run a non-parametric bootstrap using the BCa technique
      (See An introduction to the bootstrap, Efron, 1993). The BCa is
      intended for caclulation of second-order correct confidence
      intervals.


      <p class="option">bootstrap -samples=2000 -bca -stratify_on=5 run89.mod </p>

      This is the same BCa approach as above but with stratification
      on the factors of column five.
EOF

    $help_text{Options} = <<'EOF';
   <h3 class="heading1">Options:</h3>

      The options are given here in their long form. Any option may be
      abbreviated to any nonconflicting prefix. The -threads option
      may be abbreviated to <span class="option">-t</span> (or even <span class="option">-thr</span>) but -<span class="option">debug</span> may not be
      abbreviated to <span class="option">-d</span> because it conflicts with <span class="option">-debug_packages</span> and
      <span class="option">-debug_subroutines</span>.

      The following options are valid:
EOF
      
    $help_text{-samples} = <<'EOF';      
      <p class="option">-samples</p>

      The number of bootstrap samples.
EOF

    $help_text{-sample_size} = <<'EOF';
      <p class="option">-sample_size</p>

      The number of subjects in each bootstrap data set. The default
      value is set to the number of individuals in the original data
      set.

      When the resampling is stratified, the sample_size option can be
      used to specify the exact number of samples that should be drawn
      from each strata. Below follows an example of the syntax that
      should be used in such a case. Stratification is here done based
      on the study number, STUD, with the values 1001, 1002 and 1003.

      -sample_size='1001=>12,1002=>24,1003=>10'

      This example specifies that the bootstrap should use 12 samples
      from study 1001, 24 samples from 1002 and 10 from study 1003.

      If only one sample size is used together with stratified
      resampling (the default case; sample_size=number of individuals
      in the data set), the strata are assigned samples in proportion
      to their size in the data set. Please note that this usage of
      the sample_size option does not guarantee that the sum of the
      samples of the strata is equal to the given sample_size since
      PsN needs to round the figures to the closest integer. For a
      sample size equal to the number of individuals in the data set,
      the sum will however always be correct.

EOF

    $help_text{-bca} = <<'EOF';
      <p class="option">-bca</p>

      Using the <span class="option">-bca </span>option, the bootstrap
      utility will calculate the confidence intervals through the BCa
      method. The default approach however, is not to use the BCa (see
      Efron B, An introduction to the Bootstrap, 1993). The BCa is
      intended for calculation of second-order correct confidence
      intervals.  
EOF

    $help_text{-stratify_on} = <<'EOF';
      <p class="option">-stratify_on=integer|string</p>

      It may be necessary to use stratification in the resampling
      procedure. For example, if the original data consists of two
      groups of patients - say 10 patients with full pharmacokinetic
      profiles and 90 patients with sparse steady state concentration
      measurements - it may be wise to restrict the resampling
      procedure to resample within the two groups, producing bootstrap
      data sets that all contain 10 rich + 90 sparse data patients but
      with different compositions. The default is not to use
      stratification. Set <span class="option">-stratify_on</span> to
      the column that defines the two groups. If a string is used with
      stratify_on the header in the datafile is used to map the string
      to a column number.

      Note that the option sample_size has a different behavior when
      stratified resampling is used.
EOF

    $help_text{-skip_covariance_step_terminated} = <<'EOF';
      <p class="option">-skip_covariance_step_terminated</p>

      With this option enabled, the bootstrap will skip all samples
      where the NONMEM run terminated the covariance step.
EOF

    $help_text{-skip_with_covstep_warnings} = <<'EOF';
      <p class="option">-skip_with_covstep_warnings</p>

      With this option enabled, the bootstrap will skip all samples
      where the NONMEM run had warnings from the covariance step.
EOF
 
    $help_text{-skip_minimization_terminated} = <<'EOF';
      <p class="option">-skip_minimization_terminated</p>

      With this option enabled, the bootstrap will skip all samples
      where the NONMEM run terminated the minimization step.
EOF
  
    $help_text{-skip_estimate_near_boundary} = <<'EOF';
      <p class="option">-skip_estimate_near_boundary</p>

      With this option enabled, the bootstrap will skip all samples
      where the NONMEM run signal that some estimates are near its
      boundary.
EOF

    $help_text{-estimate_near_boundary_limit} = <<'EOF';
      <p class="option">-estimate_near_boundary_limit='number'</p>

      If the <span class="option">-summarize</span> or <span class="option">-summary</span> options are set, the bootstrap
      will do a set of diagnostics checks. Among other things it
      checks the ratio between runs with estimates near their
      boundaries and those without. If the ratio is to high ( by
      default 20% ) a warning is printed to the screen. You can change
      the ratio with <span class="option">-estimate_near_boundary_limit</span>.
EOF

    $help_text{-covariance_step_successful_limit} = <<'EOF';
      <p class="option">-covariance_step_successful_limit='number'</p>

      If the <span class="option">-summarize</span> or <span class="option">-summary</span> options are set, the bootstrap
      will do a set of diagnostics checks. Among other things it
      checks the ratio between runs with covariance step succesful and
      those without. If the ratio is to high ( by default 80% ) a
      warning is printed to the screen. You can change the ratio with
      <span class="option">-covariance_step_successful_limit</span>.
EOF

    $help_text{-covariance_step_warnings_limit} = <<'EOF';
      <p class="option">-covariance_step_warnings_limit='number'</p>

      If the <span class="option">-summarize</span> or <span class="option">-summary</span> options are set, the bootstrap
      will do a set of diagnostics checks. Among other things it
      checks the ratio between runs with warnings in the covariance
      step and those without. If the ratio is to high (by default 20%)
      a warning is printed to the screen. You can change the ratio
      with <span class="option">-covariance_step_warnings_limit</span>.
EOF

    $help_text{-minimization_successful_limit} = <<'EOF';
      <p class="option">-minimization_successful_limit='number'</p>

      If the <span class="option">-summarize</span> or <span class="option">-summary</span> options are set, the bootstrap
      will do a set of diagnostics checks. Among other things it
      checks the ratio between runs with successful minimizations and
      those without. If the ratio is to high (by default 80%) a
      warning is printed to the screen. You can change the ratio with
      <span class="option">-minimization_successful_limit</span>.
EOF

   $help_text{-mplots} = <<'EOF';
      <p class="option">-mplots</p>

      Generate matlab scripts for making various plots of the result.
EOF

    $help_text{-rplots} = <<'EOF';
      <p class="option">-rplots</p>

      Generate R scripts for making various plots of the result.
EOF
    $help_text{Post_help_message} = <<'EOF';
      Also see 'execute -help' for a description of common options.
EOF

common_options::online_help('bootstrap',\%options, \%help_text, \%required_options, \%optional_options);


## Check that we do have a model file
if ( scalar(@ARGV) < 1 ) {
  print "A model file must be specified. Use 'bootstrap -h' for help.\n";
  exit;
}

if( scalar(@ARGV) > 1 ){
  print "Bootstrap can only handle one modelfile. Use 'bootstrap -h' for help.\n";
  exit;
}

my $eval_string = common_options::model_parameters(\%options);

my $model = model -> new ( eval( $eval_string ),
                           filename                    => @ARGV[0],
                           ignore_missing_output_files => 1 );

if( $options{'shrinkage'} ) {
  $model -> shrinkage_stats( enabled => 1 );
}

my $type = defined $options{'bca'} ? 'bca' : undef;


if( defined $options{'sample_size'} ) {
  $options{'sample_size'} = 'default=>'.$options{'sample_size'} if( $options{'sample_size'} =~ /^\d+$/ );
} else {
  $options{'sample_size'} = 'default=>'.$model -> datas() -> [0] -> count_ind();
}
my %subj_hash = eval($options{'sample_size'});
if( $@ ) {
  die $@."\nThe sample_size option must be either a single digit or of the format 'strata1=>12,strata2=>34, ...' etc" ;
}

## Create new Bootstrap object:
my $bs = tool::bootstrap -> 
  new ( eval( $common_options::parameters ),
        top_tool => 1,
        prepend_model_file_name => 1,
        models                               => [ $model ],
        samples                              => $options{'samples'},
        subjects                             => \%subj_hash,
        type                                 => $type,
        stratify_on                          => $options{'stratify_on'},
        skip_minimization_terminated         => $options{'skip_minimization_terminated'},
        skip_covariance_step_terminated      => $options{'skip_covariance_step_terminated'},
        skip_with_covstep_warnings           => $options{'skip_with_covstep_warnings'},
        skip_estimate_near_boundary          => $options{'skip_estimate_near_boundary'},
        minimization_successful_limit        => $options{'minimization_successful_limit'},
        covariance_step_successful_limit     => $options{'covariance_step_successful_limit'},
        covariance_step_warnings_limit       => $options{'covariance_step_warnings_limit'},
        estimate_near_boundary_limit         => $options{'estimate_near_boundary_limit'},
        estimate_near_boundary_limit         => $options{'estimate_near_boundary_limit'},
        se_confidence_intervals_level         => $options{'se_confidence_intervals_level'},
        percentile_confidence_intervals_level => $options{'percentile_confidence_intervals_level'},
        bca_confidence_intervals_level        => $options{'bca_confidence_intervals_level'},
        se_confidence_intervals_check         => $options{'se_confidence_intervals_check'},
        percentile_confidence_intervals_check => $options{'percentile_confidence_intervals_check'},
        bca_confidence_intervals_check        => $options{'bca_confidence_intervals_check'},
        large_bias_limit                     => $options{'large_bias_limit'} );

open(CMD, ">", $bs -> directory . "/command.txt");
print CMD $cmd_line, "\n";
close(CMD);

if ( $options{'summarize'} ) {
  $bs -> prepare_results();
  $bs -> print_results();
  $bs -> print_summary();
} else {
  $bs -> run;
  $bs -> prepare_results();
  $bs -> print_results();
#  $bs -> print_summary();
}

if( $options{'mplots'} ){
  $bs -> create_matlab_scripts();
}

if( $options{'rplots'} ){
  $bs -> create_R_scripts();
}

ui -> print( category => 'bootstrap',
             message  => "Bootstrap done.");