options ls=80 ps = 5000;
libname lib 'c:\Documents and Settings\Owner\My Documents\Mouse Leukemia Data';

data lib.tmp;
input DoB mmddyy. @8 DoD mmddyy. @15 J 1. @17 status 1. @19 z2 1. @21 z6  5.;
cards;
121175 122377 2 2 2 10000
121175 122377 3 2 . 02400
121175 060677 3 1 . 00000
121175 121476 2 1 . 10000
121175 122377 3 2 2 00000
121175 010577 2 1 . 08800
121175 111677 5 1 2 08000
121175 101477 1 1 . 10000
121175 122377 3 2 2 00000
121175 122377 3 2 . 00080
121175 122377 1 2 2 00000
121175 040777 1 1 1 10000
121175 081076 1 1 . .    
121175 050577 1 1 . 10000
121175 081077 1 1 1 10000
121175 091277 2 1 . 00080
121175 021777 1 1 . 10000
121175 091476 1 1 . .    
121175 031877 1 1 2 10000
121275 100676 1 1 . 10000
121275 122377 2 2 2 03600
121275 122377 5 2 2 02000
121275 062877 2 1 . 00000
121275 122377 3 2 2 00000
121275 081776 1 1 . .    
121275 080676 3 1 . .    
121275 122377 3 2 2 10000
121275 041877 3 1 . .    
121275 122377 3 2 2 10000
121275 122377 3 2 1 10000
121275 081977 2 1 1 10000
121275 122377 3 2 . 01000
121275 122377 3 2 1 10000
121275 102076 1 1 . 10000
121275 122377 3 2 . 05200
121275 102077 3 1 2 03000
121275 041177 3 1 . 10000
121275 052477 5 1 1 10000
121275 122377 3 2 2 05600
121275 070677 1 1 . 10000
121575 083076 6 1 . .    
121575 083076 6 1 . .    
121575 083076 6 1 . .    
121575 083076 6 1 . .    
121575 100776 3 1 . .    
121575 082677 2 1 . 10000
121575 091476 1 1 . .    
121575 080476 1 1 . .    
121575 082677 1 1 2 00000
121575 031877 1 1 1 10000
121575 110477 4 1 . 00000
121575 050577 1 1 1 10000
121575 122377 3 2 . 00000
121575 122176 1 1 . 10000
121575 122377 3 2 2 00800
121575 122377 3 1 1 00360
121575 091277 3 1 . 10000
121575 072577 6 1 . 00160
121575 072577 6 1 . 00480
121575 063077 3 1 1 00160
121575 101276 1 1 . 10000
121575 100477 1 1 1 10000
121575 100476 1 1 . .    
121575 052077 1 1 1 10000
121575 100477 3 1 2 00000
121575 100477 3 1 2 02000
121575 072676 3 1 . 03100
121575 020977 1 1 . 10000
121675 072577 6 1 . 00280
121675 072577 6 1 . 10000
121675 072577 6 1 . 10000
121675 110876 3 1 . .    
121675 122377 3 2 . 02300
121675 122377 3 2 2 06400
121675 122377 3 2 2 10000
121775 031877 1 1 1 10000
121775 081877 1 1 . 00400
121775 091477 6 1 . 00240
121775 120276 3 1 . .    
121775 020878 2 2 2 00040
121775 122377 5 2 2 00080
121775 122377 2 2 2 01400
121775 031877 3 1 2 00000
121775 062077 1 1 . 10000
121775 060777 1 1 1 10000
121775 122377 3 2 2 00080
121775 102776 1 1 . 10000
122275 111577 3 2 2 10000
122275 111677 1 1 1 10000
122275 061477 1 1 . 10000
122275 111577 3 2 2 00960
122275 111577 3 2 2 00000
122275 092476 1 1 . .    
122275 012677 1 1 . 10000
122275 111577 3 1 2 10000
122275 122377 1 2 2 04000
122275 042677 2 1 . 00100
122275 033077 1 1 1 10000
122275 010577 1 1 . 10000
122275 122377 3 2 2 00720
122275 070677 1 1 . 10000
122275 122377 5 2 2 02400
122275 122377 3 2 2 02000
122275 122377 3 2 1 10000
122275 102776 1 1 . 10000
122275 101477 3 1 2 00880
122275 122377 3 2 2 10000
122275 122377 3 2 1 10000
122275 052777 1 1 . 10000
122275 091377 2 1 . 01500
122275 122377 3 2 . 00000
010576 012778 1 1 . 00000
010576 072677 1 1 . 10000
010576 083077 3 1 2 00560
010576 020178 3 2 . 10000
010576 060777 1 1 . 10000
011276 020178 3 2 2 00000
011276 091376 1 1 . .    
011276 020178 1 2 2 00360
011276 020178 2 2 2 03800
011276 100676 1 1 . 10000
011276 122977 2 1 2 02800
011276 020178 3 2 2 00000
011276 020178 3 2 2 01500
011276 020178 3 2 2 10000
011276 011277 3 1 . .    
011276 120676 3 1 . .    
011276 082377 3 1 . 10000
011576 020178 3 2 2 00040
012676 020178 3 2 1 00040
012676 010577 3 1 . .    
012676 030877 1 1 1 10000
012676 102076 1 1 . 10000
012676 122977 2 2 2 00280
012676 020178 3 2 2 00280
012676 100477 4 1 . 10000
012676 012677 6 1 . .    
012676 020178 1 2 . 00080
012676 120776 1 1 . 10000
012676 020178 3 2 2 00040
012676 020178 1 2 . 10000
012676 020178 3 2 2 00040
012676 021777 1 1 . 10000
012676 120877 1 2 2 09200
012676 020178 2 2 2 10000
012676 020178 3 2 2 00040
012676 031477 4 1 . .    
012676 082076 1 1 . .    
012676 122077 1 2 . 10000
012676 021777 1 1 . 00000
012676 022176 1 1 . 10000
012676 111176 1 1 . 10000
012676 020178 3 2 2 06000
012676 020178 1 2 . 00880
012676 020178 3 2 2 01800
012676 020178 3 2 . 00200
012676 011777 2 1 . .    
012676 091477 3 1 . 00400
012676 020178 3 2 1 10000
012676 111677 2 1 . 10000
012676 010177 4 1 . .    
012676 081177 3 1 1 10000
012676 020977 1 1 . 10000
012676 012778 5 2 2 01200
012676 112277 1 1 . 10000
012776 021777 1 1 . 10000
012776 111776 1 1 . 10000
012776 111176 1 1 . 10000
012776 092077 4 1 . 05200
012776 102676 1 1 . .    
012776 020178 3 2 2 04000
012776 020178 2 2 2 10000
012776 062077 1 1 . 10000
012776 020878 3 2 2 10000
012776 033077 1 1 2 10000
020276 020878 3 2 2 00040
020276 100477 3 1 . 10000
020276 030877 1 1 1 10000
020276 010178 4 1 2 00240
020276 090676 1 1 . .    
020276 020878 3 2 1 10000
020276 040777 1 1 2 10000
020276 050677 1 1 1 10000
020276 101876 2 1 . 10000
020276 101876 3 1 . .    
020276 020878 3 2 2 10000
020276 050977 3 1 . .    
020276 020878 3 2 1 10000
020276 020878 3 2 2 00040
020276 020878 3 2 2 03900
020276 050577 1 1 . 10000
020276 012177 1 1 . 10000
020276 090777 1 1 1 04200
020276 061677 1 1 . 10000
020276 072977 3 1 1 02900
020976 011078 1 1 2 09200
020976 020878 3 2 2 00600
020976 020878 3 2 2 00800
020976 020878 3 2 1 10000
020976 020878 3 2 2 10000
020976 020878 3 2 2 01000
020976 020878 3 2 1 10000
020976 020878 3 2 2 10000
020976 020878 3 2 2 02200
;
run;

data lib.a;
set lib.tmp;
time = DoD - DoB;
if z6 < 10000 then vl = 1;
if z6 >= 10000 then vl = 0;
if z6 = . then vl = .;
if z2 = 1 then gpd1 = 1;
if z2 = 2 then gpd1 = 0;
if z2 = . then gpd1 = .;
m_gpd1 = 0; 
if gpd1 = . then m_gpd1 = 1;
m_vl = 0; 
if vl = . then m_vl = 1;
tleuD = 0; 
if J = 1 and status = 1 then tleuD = 1;
leuD = 0; 
if J in (1,2) and status = 1 then leuD = 1;
if vl =  . and gpd1 =  . then m = 3;
if vl ne . and gpd1 =  . then m = 2;
if vl =  . and gpd1 ne . then m = 1;
if vl ne . and gpd1 ne . then m = 0;
output;
run;

proc tphreg data = lib.a outest=outest1 covout;
class gpd1 vl;
model time*leuD(0) = gpd1 vl;
run;

proc print data=outest1;
run;

proc tphreg data = lib.a outest=outest1 covout;
class vl;
model time*leuD(0) = vl /rl;
run;
proc print data=outest1;
run;
***************************** method ;
***************************** method ;
***************************** method ;
***************************** method ;
***************************** method ;

%macro bootstr(bootnumber = 10, rn = 123456);
 data boot;
   set lib.a;
   where gpd1 ne . and vl ne .; 
 run;
* data out;
*  input gpd1 vl gpd10 vl0;
*  cards;
* ;
* run;
%DO i = 1 %to &bootnumber;
DATA an;
   choice = INT(RANUNI(0) * n) + 1;
   SET boot POINT = choice NOBS = n;
   j1 + 1;
   IF j1 > n THEN STOP;
RUN;

proc tphreg data = an outest=outest1 noprint;
class gpd1 vl; where gpd1 ne . and vl ne .; model time*leuD(0) = gpd1 vl;
run;

data outest1; set outest1; 
gpd1 = gpd10; vl = vl0;output;
keep gpd1 vl; run;

proc tphreg data = an outest=outest2 noprint;
class gpd1 vl; where gpd1 ne . and vl ne .; model time*leuD(0) = gpd1;
run;

data outest2; set outest2; keep gpd10;
run;

proc tphreg data = an outest=outest3 noprint;
class gpd1 vl; where gpd1 ne . and vl ne .; model time*leuD(0) = vl;
run;

data outest3;  set outest3;  keep vl0;
run;

data outest;
merge outest1 outest2 outest3;
run;

print outest;

PROC APPEND BASE = out DATA = outest;quit;
%END; 
%mend;

%bootstr(bootnumber=500,rn=12);

***************************** method ;
***************************** method ;
***************************** method ;
***************************** method ;
***************************** method ;

%macro bootstr2(bootnumber = 10, rn = 123456);
 data boot;
   set lib.a;
   where vl ne .; 
 run;
%DO i = 1 %to &bootnumber;
DATA an;
   choice = INT(RANUNI(0) * n) + 1;
   SET boot POINT = choice NOBS = n;
   j1 + 1;
   IF j1 > n THEN STOP;
RUN;

proc tphreg data = an outest=outest noprint;
model time*leuD(0) = vl;
run;

PROC APPEND BASE = out DATA = outest;quit;
%END; 
%mend;

%bootstr2(bootnumber=1000,rn=12);

proc gplot data = lib.out;
 plot vl*vl0;
 plot gpd1*vl0;
* plot vl*gpd1;
run;

data out;
set out;
keep vl _LNLIKE_;
run;

proc iml; /* bootstrap results */
  use out;
    read all into boot;
  close out;
  n = nrow(boot);
  mean = mean(boot[:,]);
  var = (boot - j(n,1,1) * mean)` * (boot - j(n,1,1) * mean) / (n - 1);
  print mean;
  print var;
quit;

proc gplot data = out;
 plot vl*_LNLIKE_;
run;

data lib.out;
set out; keep gpd1 vl vl0;
if -6 <= gpd1 <= 2 and -2 <= vl <= 5 then output;
run;

proc iml; /* bootstrap results */
  use lib.out;
    read all into boot;
  close lib.out;
  n = nrow(boot);
  mean = mean(boot[:,]);
  var = (boot - j(n,1,1) * mean)` * (boot - j(n,1,1) * mean) / (n - 1);
  print mean var;
quit;

proc gplot data = lib.out;
 plot vl*gpd1;
run;

proc corr data = lib.out;
run;

proc iml;
  vl_tilde = 1.90263;
  var_vl_tilde = 0.10266;
  vl_tilde_boot = 1.980864;
  var_vl_tilde_boot =  0.1213734;
  vl_hat = 2.01730;
  var_vl_hat = 0.30978;
  theta_hat = {-1.45927 1.21888};
  var_theta_hat = {0.32112 0.18797,  
                   0.18797 0.42019};
  *print theta_hat var_theta_hat;
  theta_gpd1_vl = {-1.5578 1.5032531 -2.200564 2.3247258};/* bootstrap */
  var_theta_gpd1_vl = { 1.3022881   0.1159537   1.0587751 -0.243926,
                        0.1159537   4.4406863  -0.3370130  4.157202,
                        1.0587751  -0.3370130   1.0179679 -0.541224,
                       -0.2439260   4.1572020  -0.5412240  4.112789}; /* bootstrap var */
  theta_vl = {-1.508325 1.2639831 2.0791849}; /* selected bootstrap */
  var_theta_vl = {0.4175926 0.2590969 0.0400165,
                  0.2590969 0.4865435 0.3358946,
                  0.0400165 0.3358946 0.3253472}; /* selected bootstrap var */
  cov = {0.0400165 0.3358946};
  lambda = cov` / 0.4865435;
  print 'LAMBDA' lambda;
  ind = {1 2};
  print 'old no boot' theta_hat var_theta_hat;
  theta = theta_hat` - lambda * (vl_hat - vl_tilde);
  var_theta = var_theta_hat - lambda * var_vl_hat * lambda`
              + lambda * var_vl_tilde *  lambda`;
  theta = theta`;
  print 'improved' theta var_theta;
  theta_boot = theta_vl[ind]`;
  var_theta_boot = var_theta_vl[ind,ind];
  print 'old boot' theta_boot var_theta_boot;
  theta_boot = theta_vl[ind]` - lambda` * (theta_vl[3] - vl_tilde_boot);
  var_theta_boot = var_theta_vl[ind,ind] - lambda * var_theta_vl[2,2] * lambda`
                 + lambda * var_vl_tilde_boot *  lambda`;
  theta_boot = theta_boot`;
  print 'boot improved' theta_boot var_theta_boot;
quit;

proc iml;
  use lib.boot;
    read all into boot;
  close lib.boot;
  n = nrow(boot);
  print n;
  mean = mean(boot[:,]);
  var = (boot - j(n,1,1) * mean)` * (boot - j(n,1,1) * mean) / (n - 1);
  print mean;
  print var;
quit;