%include "MLGA.sas" ; /*Replication instructions----------------------------------------------------------------------------*/ /*Execute this program in SAS. Please ensure MLGA.sas code is placed on the Code folder---------------*/ /*Reported numbers in paper were obtained using SAS 9.2 64bits on a SunOS 5.10 Unix machine-----------*/ /*Exact numbers can differ in Windows or other platforms----------------------------------------------*/ /*----------------------------------------------------------------------------------------------------*/ /*-----------------------------------------------------------------------------------------------------*/ /*Code 1 Michaelis Menten Kinetics---------------------------------------------------------------------*/ /*Referenced in Section 4.1----------------------------------------------------------------------------*/ /*Used to produce paragraph that describes how PROC NLMIXED misses the restimation and also Table 1----*/ /*-----------------------------------------------------------------------------------------------------*/ data sIMU_DATA_MM ; do day = 1 to 30 ; Vmax = 158 ; Km = 0.07 ; S = 1 * day + ranuni(10) ; y = normal(13) + (Vmax*S)/(Km+S) ; output ; end ; run ; data coefi ; length coefficient $50. ; coefficient = "Vmax" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "Km" ; lbound = 0.1 ; ubound = 500 ; output ; coefficient = "sd2" ; lbound = 0.2 ; ubound = 10 ; output ; run ; data list_var ; length var $50. ; var = "S"; output ; run ; %MLGA( data = SIMU_DATA_MM , var_list = list_var , bounds = coefi , expr = "y = ( Vmax * S ) / ( Km + S )" , out_table = XDG , seed = 11 ); proc nlmixed data = SIMU_DATA_MM ; parms Vmax = 40 Km = 40 s2 = 1.5; pred = ( Vmax * S )/( Km + S ); model y~ normal(pred , s2) ; run ; /*-------------------------------------------------------------------------------------------------*/ /*Code 2 Emax model--------------------------------------------------------------------------------*/ /*Referenced in Section 4.2------------------------------------------------------------------------*/ /*Used to produce paragraph that describes how PROC NLMIXED misses the estimation and Table 2------*/ /*-------------------------------------------------------------------------------------------------*/ data SIMU_DATA_EM ; do D = 1 to 1000 ; Emax = 37 ; E0 = 183 ; beta = 6.22 ; ED50 = 5.06 ; dose = 80 * ranuni(13) ; predx = normal(13) + E0 + (( Emax * ( Dose**beta ))/(( ED50**beta ) + ( Dose**beta ))) ; output ; end ; run ; data coefi ; length coefficient $40.; coefficient = "beta" ; lbound = 0.1 ; ubound = 50 ; output; coefficient = "E0" ; lbound = 0.1 ; ubound = 1000 ; output; coefficient = "ED50" ; lbound = 0.1 ; ubound = 1000 ; output; coefficient = "Emax" ; lbound = 0.1 ; ubound = 1000 ; output; coefficient = "sd2" ; lbound = 0.2 ; ubound = 10 ; output; run; data list_var ; length var $20. ; var = "Dose" ; output ; run ; %MLGA( data = SIMU_DATA_EM , var_list = list_var , bounds = coefi , expr = "predx = E0 + (( Emax * ( Dose**beta ))/(( ED50**beta ) + ( Dose**beta )))" , out_table = XDG , seed = 11 ); proc nlmixed data = SIMU_DATA_EM ; parm e0 = 250 beta = 0.1 ed50 = 2.06 emax = 130 s2 = 1 ; pred = E0 + (( Emax * ( Dose**beta ))/(( ED50**beta ) + ( Dose**beta ))) ; model predx~ normal(pred , s2) ; run ; proc nlmixed data = SIMU_DATA_EM ; parm e0 = 183 beta = 6.22 ed50 = 5.06 emax = 130 s2 = 1 ; pred = E0 + (( Emax * ( Dose**beta ))/(( ED50**beta ) + ( Dose**beta ))) ; model predx~ normal(pred , s2) ; run ; /*-------------------------------------------------------------------------------------------------*/ /*Code3 Restrictions in Biexponential regression---------------------------------------------------*/ /*Referenced in Section 5--------------------------------------------------------------------------*/ /*Used to produce Table 3--------------------------------------------------------------------------*/ /*-------------------------------------------------------------------------------------------------*/ data data_biexp(keep = predx dose) ; do D = 1 to 1000 ; dose = 0.1*D + 2*ranuni(13) ; A = 70 ; B = 50 ; alpha1 = 0.5 ; alpha2 = 0.05 ; predx = normal(13) + A*exp(-alpha1*Dose) + B*exp(-alpha2*Dose) ; output ; end ; run ; data coefi ; length coefficient $40. ; coefficient = "alpha1" ; lbound = 0.0 ; ubound = 100 ; output ; coefficient = "alpha2" ; lbound = 0.0 ; ubound = 100 ; output ; coefficient = "A" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "B" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "sd2" ; lbound = 0.2 ; ubound = 10 ; output ; run ; data list_var ; length var $20. ; var = "Dose" ; output ; run ; data restrictions ; restriction = "alpha1>alpha2" ; output ; run ; %MLGA( data = data_biexp, var_list = list_var , bounds = coefi , expr = "predx = A * exp(-alpha1 * Dose) + B * exp(-alpha2 * Dose)" , out_table = XDG , restrict = restrictions , gen_size = 800 , seed = 11 ); proc nlmixed data = data_biexp; parms A = 3 B = 3 alpha1 = 0.01 alpha2 = 0.01 s2 = 1 ; pred = A * exp(-alpha1 * Dose) + B * exp(-alpha2 * Dose) ; model predx~ normal(pred , s2) ; run ; /*-------------------------------------------------------------------------------------------------*/ /*Code4 Genetic Pressure---------------------------------------------------------------------------*/ /*Referenced in Section 6--------------------------------------------------------------------------*/ /*Used to produce Table 4--------------------------------------------------------------------------*/ /*-------------------------------------------------------------------------------------------------*/ data SIMU_DATA_EM ; do D = 1 to 1000 ; Emax = 37 ; E0 = 183 ; beta = 6.22 ; ED50 = 5.06 ; dose = 80*ranuni(13) ; predx = normal(13) + E0 + ((Emax*(Dose**beta))/((ED50**beta) + (Dose**beta))) ; output ; end ; run ; data coefi ; length coefficient $40. ; coefficient = "beta" ; lbound = 0.1 ; ubound = 50 ; output ; coefficient = "E0" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "ED50" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "Emax" ; lbound = 0.1 ; ubound = 1000 ; output ; coefficient = "sd2" ; lbound = 0.2 ; ubound = 10 ; output ; run; data list_var ; length var $20. ; var = "Dose" ; output ; run ; %MLGA( data = SIMU_DATA_EM , var_list = list_var , bounds = coefi , expr = "predx = E0 + (( Emax * ( Dose ** beta ))/(( ED50**beta ) + ( Dose**beta )))" , out_table = XDG , seed = 11 ); %MLGA( data = SIMU_DATA_EM , var_list = list_var , bounds = coefi , expr = "predx = E0 + (( Emax * ( Dose**beta ))/(( ED50**beta ) + ( Dose**beta )))" , out_table = XDG , seed = 11 , gapressure = 2 );