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