bin/se_of_eta

   1 #!/usr/local/bin/perl
   2
   3 # Only for Development
   4 use FindBin qw($Bin);
   5 use lib "$Bin/../PsN_2_2_0_patches_serial/lib";
   6
   7 # Don't edit the line below, it must look exactly like this.
   8 # Everything above this line will be replaced #
   9
  10 # Perl includes #
  11 use strict;
  12 use Getopt::Long;
  13 # External modules #
  14 use Math::Random;
  15 # PsN includes #
  16 use PsN;
  17 use tool::modelfit;
  18 use model;
  19 use debug;
  20 use common_options;
  21 use ui;
  22 use File::Copy 'cp';
  23 use Data::Dumper;
  24
  25 my $cmd_line = $0 . " " . join( " ", @ARGV );
  26
  27 my %options;
  28
  29 my %optional_options = ( 'dv_column:s'
  30                         );
  31
  32 my $res = GetOptions( \%options,
  33                       @common_options::get_opt_strings,
  34                       keys(%optional_options) );
  35
  36 exit unless $res;
  37
  38 common_options::set_globals( \%options, 'se_of_eta' );
  39 common_options::get_defaults( \%options, 'se_of_eta' );
  40 common_options::sanity_checks( \%options, 'se_of_eta' );
  41
  42 my %help_text;
  43 $help_text{Pre_help_message} = <<'EOF';
  44     <h3 class="heading1">se_of_eta</h3>
  45
  46       Perl script for getting standard errors of etas.
  47
  48     <h3 class="heading1">Usage:</h3>
  49 EOF
  50
  51     $help_text{Description} = <<'EOF';
  52    <h3 class="heading1">Description:</h3>
  53
  54       se_of_eta requires "MAXEVALS" to be abreviated to "MAX" if it
  55       exists in the model file.
  56
  57 EOF
  58     $help_text{Examples} = <<'EOF';
  59    <h3 class="heading1">Example:</h3>
  60
  61       <p class="style2">se_of_eta pheno.mod</p>
  62
  63        This will create a new directory se_of_eta_dirX where X is a
  64        number increased every time you run the program. Inside that
  65        directory it will create a new model file for each OMEGA in the
  66        file. You will be able to retrieve the results from the
  67        corresponding table file.
  68 EOF
  69
  70 common_options::online_help('se_of_eta',\%options, \%help_text, {}, {});
  71
  72 if ( scalar( @ARGV ) < 1 ) {
  73   print "A model file must be specified. Use '$0 -h' for help.\n";
  74   exit;
  75 }
  76
  77 unless( defined $options{'dv_column'} ) {
  78   $options{'dv_column'} = 'DV';
  79 }
  80
  81 my( $base_dir, $dummy) = OSspecific::absolute_path( '', '' );
  82
  83 my $directory;
  84
  85 if ( defined $options{'directory'} ) {
  86   my $dummy;
  87   ( $directory, $dummy ) = OSspecific::absolute_path( $options{'directory'} );
  88 } else {
  89   $directory = OSspecific::unique_path( 'se_of_eta_dir' ,
  90                                         $base_dir );
  91 }
  92
  93 unless( -d $directory ){
  94   mkdir( $directory );
  95
  96   open(CMD, ">", $directory . "/command.txt");
  97   print CMD $cmd_line, "\n";
  98   close(CMD);
  99 }
 100
 101 my $eval_string = common_options::model_parameters(\%options);
 102
 103 my $model;
 104
 105 $model = model -> new ( eval( $eval_string ),
 106                         filename                    => $ARGV[0],
 107                         ignore_missing_output_files => 1 );
 108
 109 my $datafile_name = $model -> datas -> [0] -> full_name;
 110
 111 chdir( $directory );
 112
 113 my $modelfit;
 114
 115 $modelfit = tool::modelfit ->
 116     new ( eval( $common_options::parameters ),
 117           directory => 'orig_modelfit_dir',
 118           models => [$model] );
 119
 120 $modelfit -> run;
 121
 122 my $nr_of_old_etas = $model -> nomegas -> [0];
 123 my $nr_of_old_sigmas = $model -> nsigmas -> [0];
 124 my $nr_of_old_thetas = $model -> nthetas;
 125
 126 my %seen_eta;
 127 my %seen_eps;
 128 my %seen_err;
 129
 130 my @models;
 131
 132 my $data = $model -> datas -> [0];
 133
 134 $data -> synchronize;
 135
 136 my $ids = $data -> individuals;
 137
 138 my $nrec = scalar @{$ids -> [0] -> subject_data};
 139
 140 foreach my $id (0..$#{$ids}) {
 141
 142   my $model_copy = $model -> copy( filename => 'se_model_'.($id+1).'.mod',
 143                                    copy_data => ($id+1), # Trick to copy data once
 144                                    data_file_names => ['se_model.dta']);
 145
 146   push(@models,$model_copy);
 147
 148   $model_copy -> datas -> [0] -> synchronize;
 149   $model_copy -> datas -> [0] -> filename('se_model.dta');
 150   $model_copy -> datas -> [0] -> directory( $directory );
 151   $model_copy -> datafiles(new_names => ['se_model.dta']);
 152
 153   $model_copy -> update_inits( from_output => $model -> outputs -> [0] );
 154
 155   # 1 Change all ETA(n) to THETA(n + $nr_of_old_thetas)
 156
 157   my $code_block;
 158   my $pk = 0;
 159   # Find PK or PRED block
 160   if( defined $model_copy -> pk ){
 161     $code_block = $model_copy -> pk;
 162     $pk = 1;
 163   } elsif( defined $model_copy -> pred ){
 164     $code_block = $model_copy -> pred;
 165   } else {
 166     die "Error: No \$PK or \$PRED found... \n";
 167   }
 168
 169   for( my $row = 0; $row <= $#{$code_block}; $row++ ){
 170
 171     # High-tech RegExp to find ETA(n) and put n in $dummy. Then update
 172     # $dummy and replace ETA(n) with THETA($dummy). Also register all
 173     # n in ETA(n) by storing them in %seen_eta.
 174
 175     $code_block -> [$row] =~ s/
 176                                 ([^H])ETA\((\d+)
 177                                 (?{$seen_eta{$2}=1;$dummy=$2+$nr_of_old_thetas;})\)
 178                               /$1THETA($dummy)/gx;
 179   }
 180
 181   my $nr_of_new_thetas = $nr_of_old_thetas + scalar keys %seen_eta;
 182
 183   if( $pk ){
 184     $model_copy -> pk( new_pk => $code_block );
 185   } else {
 186     $model_copy -> pred( new_pred => $code_block );
 187   }
 188
 189   # Fixate all old thetas and omegas.
 190
 191   $model_copy -> fixed( parameter_type => 'theta',
 192                         new_values => [[(1)x $nr_of_old_thetas]] );
 193
 194   # Set low(0.01) initial values for the new thetas.
 195
 196   $model_copy -> initial_values( parameter_type => 'theta',
 197                                  add_if_absent => 1,
 198                                  parameter_numbers => [[($nr_of_old_thetas+1)..($nr_of_new_thetas)]],
 199                                  new_values => [[(0.01) x (($nr_of_new_thetas) - $nr_of_old_thetas)]]);
 200
 201   # Change all sigmas to omegas
 202
 203   my $sigmas = $model_copy -> initial_values( parameter_type => 'sigma' );
 204
 205   $model_copy -> remove_records( type => 'sigma' );
 206
 207   my $nr_of_old_etas_with_corr = $model_copy -> nomegas( with_correlations => 1 ) -> [0];
 208
 209   $model_copy -> initial_values( parameter_type => 'omega',
 210                                  add_if_absent => 1,
 211                                  parameter_numbers => [[($nr_of_old_etas_with_corr+1)..($nr_of_old_etas_with_corr+$nr_of_old_sigmas)]],
 212                                  new_values => $sigmas );
 213
 214   $model_copy -> fixed( parameter_type => 'omega',
 215                         new_values => [[(1) x ($nr_of_old_etas_with_corr+$nr_of_old_sigmas)]] );
 216
 217   # Add EBAY column to dataset and change ID to IID
 218
 219   my $input = $model_copy -> record( record_name => 'input' );
 220
 221 #  $input->[0]->[0] =~ s/ID/IID/;
 222   $input->[0]->[0] =~ s/\$INPUT//;
 223   chomp( $input->[0]->[ $#{$input->[0]} ] );
 224   $input->[0]->[ $#{$input->[0]} ] .= ' EBAY L2';
 225
 226   $model_copy -> set_records( type => 'input',
 227                               record_strings =>  $input->[0]);
 228
 229   # Get error block to replace EPS(n) with ETA(n + $nr_of_old_eta) and
 230   # replace ERR(n) with ERR(n+$nr_of_old_eta)
 231
 232   my $error = $model_copy -> record( record_name => 'error' );
 233
 234   for( my $row = 0 ; $row <= $#{$error->[0]}; $row ++ ){
 235
 236     $error -> [0] -> [$row] =~ s/
 237                                 (.)EPS\((\d+)
 238                                 (?{$seen_eps{$2}=1;$dummy=$2+$nr_of_old_etas;})\)
 239                               /$1ETA($dummy)/gx;
 240
 241     $error -> [0] -> [$row] =~ s/
 242                                 (.)ERR\((\d+)
 243                                 (?{$seen_err{$2}=1;$dummy=$2+$nr_of_old_etas;})\)
 244                               /$1ETA($dummy)/gx;
 245
 246
 247     if( $error -> [0] -> [$row] =~ /(.*)Y(\s*)=(.*)/ ){
 248
 249       @{$error -> [0]} = ( @{$error -> [0]}[0..($row-1)],
 250                            'IF( EBAY.LE.0 ) ' . $error -> [0] -> [$row],
 251                            @{$error -> [0]}[$row+1..$#{$error->[0]}] );
 252     }
 253
 254   }
 255
 256   # my $nr_of_new_etas = $nr_of_old_etas + scalar keys %seen_eps;
 257
 258   foreach( 1..$nr_of_old_etas ){
 259     push(@{$error -> [0]}, '       IF( EBAY.EQ.'.($_).' ) Y = THETA(' .
 260                            ($nr_of_old_thetas+$_) .
 261                            ") + ETA($_)\n");
 262   }
 263
 264   $model_copy -> set_records( type => 'error',
 265                               record_strings => $error -> [0] );
 266
 267   # Set options in $EST and $COV
 268
 269   $model_copy -> maxeval( new_values  => [[9999]] );
 270
 271   if( scalar @{$model_copy -> record( record_name => 'covariance' )} > 0 ){
 272     $model_copy -> remove_records( type => 'covariance' );
 273   }
 274
 275   $model_copy -> add_records( type => 'covariance',
 276                               record_strings => ['MAT=R'] );
 277
 278   $model_copy -> set_records( type => 'estimation' );
 279
 280   $model_copy -> remove_records( type => 'table' );
 281
 282 #  $model_copy -> set_option( record_name => 'data',
 283 #                            option_name => 'NREC',
 284 #                            option_value => $nrec );
 285
 286   $model_copy -> set_option( record_name => 'data',
 287                              option_name => 'ACCEPT',
 288                              option_value => '(ID.EQ.'.$ids -> [$id] -> idnumber().')' );
 289
 290   $model_copy -> _write;
 291
 292 }
 293
 294 my $data_copy = $models[0] -> datas -> [0];
 295
 296 my $ids = $data_copy -> individuals;
 297
 298 my $omegas = $models[0] -> problems -> [0] -> omegas();
 299 my @L2_column;
 300 my $L2_value = 0;
 301
 302 foreach my $omega( @{$omegas} ){
 303   $L2_value++;
 304   for( 1..$omega -> size() ){
 305     push( @L2_column, $L2_value );
 306   }
 307 }
 308
 309 foreach my $id( @{$ids} ){
 310   my @last_row_of_id = split(/,/,$id -> subject_data -> [$#{$id -> subject_data}]);
 311   my $DV_column = $data_copy -> column_head_indices -> {$options{'dv_column'}};
 312
 313   my $L2_first_value=0;
 314   foreach my $row( @{$id -> subject_data} ){
 315     $L2_first_value++;
 316     $row .= ",0,$L2_first_value";
 317   }
 318
 319   for( 1..$nr_of_old_etas ) {
 320     $last_row_of_id[$DV_column-1] = '0';
 321     push( @{$id -> subject_data}, join(',', @last_row_of_id, $_, $L2_first_value + @L2_column[$_ - 1]  ) );
 322   }
 323 }
 324
 325 push( @{$data_copy -> {'header'}}, 'EBAY', 'L2' );
 326 $data_copy -> synchronize;
 327
 328 $data_copy -> _write;
 329
 330 $modelfit = tool::modelfit ->
 331     new ( eval( $common_options::parameters ),
 332           directory => 'se_of_eta',
 333           models => \@models );
 334
 335 $modelfit -> run;
 336
 337 # New Results section 9:th of June.
 338
 339 my $num = scalar @{$ids};
 340
 341 my (@theta,@setheta);
 342 my $ofv;
 343
 344 open( FILE1,'>termination.csv' );
 345 open( FILE2,'>parameters.csv' );
 346
 347 print FILE1 "ID, OFV\n";
 348
 349 print FILE2 "ID,";
 350
 351 foreach my $i( 1..$nr_of_old_thetas ){
 352   print FILE2 "TH$i,";
 353 }
 354
 355 foreach my $i( 1..$nr_of_old_etas ){
 356   print FILE2 "ETA$i,";
 357 }
 358
 359 foreach my $i( 1..$nr_of_old_thetas ){
 360   print FILE2 "SEth$i,";
 361 }
 362
 363 foreach my $i( 1..$nr_of_old_etas ){
 364   if( $i == $nr_of_old_etas ){
 365     print FILE2 "SEeta$i\n";
 366   } else {
 367     print FILE2 "SEeta$i,";
 368   }
 369 }
 370
 371 foreach my $i(1..$num) {
 372   $ofv     = get_ofv($i);
 373   @theta   = get_theta($i);
 374   @setheta = get_setheta($i);
 375
 376   print FILE1 "ID$i,$ofv\n";
 377   print FILE2 "ID$i,", join( ',', @theta,@setheta ), "\n";
 378 }
 379
 380
 381 #####################################################################################################################################
 382 ###                                                         Subroutines                                                           ###
 383 #####################################################################################################################################
 384
 385 ### Retrieves the ofv from a NONMEM-output file using function "ofv" in PsN ###
 386 sub get_ofv {
 387   my ($n) = @_;
 388   my $lst = new output ('filename' => "se_model_$n.lst");
 389   my $obf = $lst -> ofv;
 390   return $obf -> [0][0];
 391 }
 392
 393 ### Retreives the estimated thetas from a NONMEM-output file using function "theta" in PsN ###
 394 sub get_theta {
 395   my ($n) = @_;
 396   my $lst = new output ('filename' => "se_model_$n.lst");
 397   my $theta = $lst -> thetas;
 398   if(defined $theta and defined $theta->[0] and defined $theta -> [0][0]){ return @{$theta -> [0][0]}; }
 399   else { return ('.','.','.'); }
 400 }
 401
 402 ### Retreives the SEs estimated thetas from a NONMEM-output file using function "sethetas" in PsN ###
 403 sub get_setheta {
 404   my ($n) = @_;
 405   my $lst = new output ('filename' => "se_model_$n.lst");
 406   my $setheta = $lst -> sethetas;
 407   if(defined $setheta and defined $setheta->[0] and defined $setheta -> [0][0]){ return @{$setheta -> [0][0]}; }
 408   else { return ('.','.','.'); }
 409 }
 410
 411
 412