bin/llp

   1 #!/usr/local/bin/perl
   2
   3 use FindBin qw($Bin);
   4 use lib "$Bin/../lib";
   5
   6 # Don't edit the line below, it must look exactly like this.
   7 # Everything above this line will be replaced #
   8
   9 use PsN;
  10 use model;
  11 use tool::llp;
  12 use debug;
  13 use strict;
  14 use Getopt::Long;
  15 use common_options;
  16 use ui;
  17 use Data::Dumper;
  18
  19 my $cmd_line = $0 . " " . join( " ", @ARGV );
  20
  21 my %options;
  22
  23 my %required_options = ( "thetas:s"=>'theta list',
  24                          "omegas:s"=>'omega list',
  25                          "sigmas:s"=>'sigma list' );
  26
  27 my %optional_options = ( "max_iterations:i" => '',
  28                          "mplots"=>'',
  29                          "rplots"=>'',
  30                          "normq:f" => '',
  31                          "outputfile:s" => '',
  32                          "ofv_increase:f" => '',
  33                          "significant_digits:i" => '',
  34                          "rse_thetas:s" => 'theta rel. SE list',
  35                          "rse_omegas:s" => 'omega rel. SE list',
  36                          "rse_sigmas:s" => 'sigma rel. SE list',
  37                          "theta_interval_ratio_check:f" => '',
  38                          "omega_interval_ratio_check:f" => '',
  39                          "sigma_interval_ratio_check:f" => '',
  40                          "within_interval_check:f" => '' );
  41
  42 my $res = GetOptions( \%options,
  43                       @common_options::get_opt_strings,
  44                       keys(%required_options),
  45                       keys(%optional_options) );
  46
  47 exit unless $res;
  48
  49 common_options::set_globals( \%options );
  50 common_options::get_defaults( \%options, 'llp' );
  51
  52 my %help_text;
  53
  54 $help_text{Pre_help_message} = <<'EOF';
  55     <h3 class="heading1">llp</h3>
  56
  57       Perl script for log-likelihood profiling of NONMEM runs.
  58
  59     <h3 class="heading1">Usage:</h3>
  60 EOF
  61
  62 $help_text{Description} = <<'EOF';
  63     <h3 class="heading1">Description:</h3>
  64
  65       The Log-likelihood Profiling tool can be used to assess
  66       confidence interval limits for parameter estimates. The
  67       -2*log-likelihood of hierarchical models are chi-square
  68       distributed. Fixing a parameter reduces the number of parameters
  69       of the model by one. To be able to say, for a given level of
  70       confidence, that there is a higher likelihood that the data has
  71       been produced by a system described by the full model than by one
  72       described by the reduced, the difference in the -2*log-likelihood
  73       should be at least X. For example, using a confidence level of
  74       95%, the difference (or X above) should be at least 3.84. The
  75       minimal number of arguments include a modelfile name and a
  76       listing of parameters, given that an output file with standard
  77       error estimates exist.
  78 EOF
  79
  80 $help_text{Examples} = <<'EOF';
  81     <h3 class="heading1">Example:</h3>
  82
  83       <p class="style2">llp -model=run89.mod -thetas='1,2'</p>
  84
  85       This will make the llp tool try to estimate the confidence
  86       intervals for thetas one and two of the model in run89.mod. It
  87       will base the first guesses on the standard error estimates from
  88       run89.lst.
  89
  90       <p class="style2">llp -model=run89.mod -thetas='1,2' -rse_thetas='20,30'</p>
  91
  92       In this example, we explicitly specify the relative standard
  93       errors which is necessary if we do not have an output file with
  94       standard error estimates.
  95 EOF
  96
  97 $help_text{Options} = <<'EOF';
  98     <h3 class="heading1">Options:</h3>
  99
 100       The options are given here in their long form. Any option may be
 101       abbreviated to any nonconflicting prefix. The <span class="style2">-threads</span> option
 102       may be abbreviated to <span class="style2">-t</span>(or even <span class="style2">-thr</span>) but <span class="style2">-debug</span> may not be
 103       abbreviated to <span class="style2">-d</span> because it conflicts with <span class="style2">-debug_packages</span> and
 104       <span class="style2">-debug_subroutines</span>.
 105     <br><br>
 106       The following options are valid:
 107 EOF
 108
 109 $help_text{-h} = <<'EOF';
 110       <p class="style2">-h | -?</p>
 111
 112       With -h or -? llp will print a list of options and exit.
 113 EOF
 114
 115 $help_text{-help} = <<'EOF';
 116       <p class="style2">-help</p>
 117
 118       With -help llp will print this, longer, help message.
 119 EOF
 120
 121 $help_text{-model} = <<'EOF';
 122       <p class="style2">-model='filename'</p>
 123
 124       The name of the model file. May be specified with or without the
 125       "'".
 126 EOF
 127
 128 $help_text{-outputfile} = <<'EOF';
 129       <p class="style2">-outputfile='filename'</p>
 130
 131       The name of the NONMEM output file. The default value is the
 132       name of the model file with a '.mod' substituted with
 133       '.lst'. Example: if the modelfile is run89.mod, the default name
 134       of the output file is run89.lst. If the name of the modelfile is
 135       cmd123 the default name of the output file is cmd123.lst. If the
 136       name of your output file does not follow this standard, you have
 137       to specify it with this option.
 138 EOF
 139
 140 $help_text{-thetas} = <<'EOF';
 141       <p class="style2">-thetas, -omegas, -sigmas='comma-separated list of parameter
 142       numbers'</p>
 143
 144       Specifies the parameters for which the llp should try to assess
 145       confidence intervals.
 146 EOF
 147
 148 $help_text{-rse_thetas} = <<'EOF';
 149       <p class="style2">-rse_thetas, -rse_omegas, -rse_sigmas='comma-separated list of
 150       relative standard errors'</p>
 151
 152       The relative standard errors should be specified in percent (%).
 153 EOF
 154
 155 $help_text{-max_iterations} = <<'EOF';
 156       <p class="style2">-max_iterations=integer</p>
 157
 158       This number limits the number of search iterations for each
 159       interval limit. If the llp has not found the upper limit for a
 160       parameter after max_iteration number of guesses it
 161       terminates. The default value is 10.
 162 EOF
 163
 164     $help_text{-significant_digits} = <<'EOF';
 165       <p class="style2">-significant_digits=integer</p>
 166
 167       Specifies the number of significant digits that is required for
 168       the test of the increase in objective function value.
 169
 170 EOF
 171     $help_text{-normq} = <<'EOF';
 172       <p class="style2">-normq=number</p>
 173
 174       This number is used for calculating the first guess of the
 175       confidence interval limits. If the standard errors exist, the
 176       first guess will be
 177
 178       <p class="style2">maximum-likelihood estimate ± normq * standard error</p>
 179
 180       otherwise it will be approximated with
 181
 182       <p class="style2">maximum-likelihood estimate ± normq * rse_parameter/100 * maximum-likelihood estimate</p>
 183
 184       where rse_parameter is rse_thetas, rse_omegas or rse_sigmas. The
 185       default value is 1.96 which translates a 95% confidence interval
 186       assuming normal distribution of the parameter estimates.
 187 EOF
 188
 189     $help_text{-ofv_increase} = <<'EOF';
 190       <p class="style2">-ofv_increase</p>
 191
 192       The increase in objective function value associated with the
 193       desired confidence interval. The default value is 3.84.
 194 EOF
 195
 196     $help_text{-mplots} = <<'EOF';
 197       <p class="style2">-mplots</p>
 198
 199       Generate matlab scripts for making various plots of the result.
 200 EOF
 201
 202     $help_text{-rplots} = <<'EOF';
 203       <p class="style2">-rplots</p>
 204
 205       Generate R scripts for making various plots of the result.
 206 EOF
 207
 208     $help_text{Post_help_message} = <<'EOF';
 209       Also see 'execute -help' for a description of common options.
 210 EOF
 211
 212     $help_text{-omegas} = <<'EOF';
 213       <p class="style2">-omegas</p>
 214
 215       Se corresponding help text for thetas.
 216 EOF
 217
 218     $help_text{-sigmas} = <<'EOF';
 219       <p class="style2">-sigmas</p>
 220
 221       Se corresponding help text for thetas.
 222 EOF
 223
 224     $help_text{-rse_omegas} = <<'EOF';
 225       <p class="style2">-rse_omegas</p>
 226
 227       Se corresponding help text for thetas.
 228 EOF
 229
 230     $help_text{-rse_sigmas} = <<'EOF';
 231       <p class="style2">-rse_sigmas</p>
 232
 233       Se corresponding help text for thetas.
 234 EOF
 235
 236 common_options::online_help( 'llp', \%options, \%help_text, \%required_options, \%optional_options);
 237
 238 ## Check that we do have a model file
 239 if ( scalar(@ARGV) < 1 ) {
 240   print "At least on model file must be specified. Use 'llp -h' for help.\n";
 241   exit;
 242 }
 243
 244 if( scalar(@ARGV) > 1 ){
 245   print "LLP can only handle one modelfile. Use 'llp -h' for help.\n";
 246   exit;
 247 }
 248
 249 unless( $options{'thetas'} or $options{'omegas'} or $options{'sigmas'}){
 250   print "You must specify one of the '--thetas', '--omegas' or '--sigmas' options\n";
 251   exit;
 252 }
 253
 254 ui -> category( 'llp' );
 255 ui -> silent(1) if( $options{'silent'} );
 256
 257 debug -> level( $options{'debug'} );
 258 debug -> package( $options{'debug_package'} );
 259 debug -> subroutine( $options{'debug_subroutine'} );
 260 debug -> warn_with_trace( $options{'warn_with_trace'} );
 261
 262 my @thetas     = split( ',',$options{'thetas'} );
 263 my @omegas     = split( ',',$options{'omegas'} );
 264 my @sigmas     = split( ',',$options{'sigmas'} );
 265 my @rse_thetas = split( ',',$options{'rse_thetas'} );
 266 my @rse_omegas = split( ',',$options{'rse_omegas'} );
 267 my @rse_sigmas = split( ',',$options{'rse_sigmas'} );
 268
 269 my ( %checked_rse_thetas, %checked_rse_omegas, %checked_rse_sigmas );
 270
 271 foreach my $param ( 'thetas', 'omegas', 'sigmas' ) {
 272   my $nse = eval('$#rse_'.$param)+1;
 273   if ( $nse > 0 ) {
 274     my $npa = eval('$#'.$param)+1;
 275     die "The number of $param ($npa) does not match the number of relative standard errors ($nse)\n"
 276       if ( not $npa == $nse );
 277     for ( my $i = 0; $i < $npa; $i++ ) {
 278       eval( '$checked_rse_'.$param.'{$'.$param.'['.$i.']} = $rse_'.$param."[$i]" ),"\n";
 279     }
 280   }
 281 }
 282
 283 my $eval_string = common_options::model_parameters(\%options);
 284
 285 my $model = model -> new ( eval( $eval_string ),
 286                            filename                    => $ARGV[0],
 287                            ignore_missing_output_files => 1);
 288
 289 if( $options{'nonparametric_etas'} or
 290     $options{'nonparametric_marginals'} ) {
 291   $model -> add_nonparametric_code;
 292 }
 293
 294 if( $options{'shrinkage'} ) {
 295   $model -> shrinkage_stats( enabled => 1 );
 296 }
 297
 298 ## Create new Llp object:
 299 my $llp = tool::llp ->
 300   new ( eval( $common_options::parameters ),
 301         max_iterations     => $options{'max_iterations'},
 302         models             => [ $model ],
 303         normq              => $options{'normq'},
 304         ofv_increase       => $options{'ofv_increase'},
 305         significant_digits => $options{'significant_digits'},
 306         run_thetas         => \@thetas,
 307         run_omegas         => \@omegas,
 308         run_sigmas         => \@sigmas,
 309         rse_thetas         => \%checked_rse_thetas,
 310         rse_omegas         => \%checked_rse_omegas,
 311         rse_sigmas         => \%checked_rse_sigmas,
 312         theta_interval_ratio_check => $options{'theta_interval_ratio_check'},
 313         omega_interval_ratio_check => $options{'omega_interval_ratio_check'},
 314         sigma_interval_ratio_check => $options{'sigma_interval_ratio_check'},
 315         within_interval_check => $options{'within_interval_check'} );
 316
 317 open(CMD, ">", $llp -> directory . "/command.txt");
 318 print CMD $cmd_line, "\n";
 319 close(CMD);
 320
 321 if ( $options{'summarize'} ) {
 322   $llp -> prepare_results();
 323   $llp -> print_summary;
 324 } else {
 325   $llp -> run;
 326   $llp -> prepare_results();
 327   $llp -> print_results;
 328 #  $llp -> print_summary;
 329 }
 330
 331 if( $options{'mplots'} ){
 332   $llp -> create_matlab_scripts();
 333 }
 334
 335 if( $options{'rplots'} ){
 336   $llp -> create_R_scripts();
 337 }
 338
 339
 340
 341
 342
 343