/*the first i1 missing observations in case group are imputed as im11, the first i2 missing obs in control group are imputed 
as im01. Calutated the estimate for beta_x and delta_x and their standard error. Results are saved in dataset _est_*/
 %macro impute(im11=,/*the value impute for the ith value when it is missing and in case group*/ 
                           im12=,/*value for other observation in case group*/
                           im01=,/*the value impute for the first ith value when x is missing and in control group*/
                           im02=/*value for other observation in control group*/);
data _miss1new_;
set _miss1_;
if _N_<=&i1  then x=&im11;
else x=&im12;
run;
data _miss0new_;
set _miss0_;
if _N_<=&i2  then x=&im01;
else x=&im02;
run;
data _new_;
set _miss0new_ _miss1new_ _ob_;
_time_=2-y;
run;
proc phreg data=_new_ outest=_phout_ covout noprint;
 model _time_*y(0)=w x;
 strata mid;run;
data _phout0_; set _phout_; keep x _NAME_;run;
data _betax_ _sx_;
set _phout0_;
select (_NAME_);
when('x') output _sx_;
when('_time_') output _betax_;
otherwise;
end;
drop _NAME_;
run;
data _betax_; set _betax_; betax=x;drop x;run;
data _sx_; set _sx_; sx=sqrt(x); drop x;run;
proc logistic data=_new_ descending outest=_logitout_ covout noprint;
model r= x y;
run;
data _logit0_; set _logitout_; keep x _NAME_;run;
data _deltax_ _dsx_;
set _logit0_;
select (_NAME_);
when('x') output _dsx_;
when('r') output _deltax_;
otherwise;
end;
drop _NAME_;
run;
data _deltax_; set _deltax_; deltax=x;drop x;run;
data _dsx_; set _dsx_; dsx=sqrt(x); drop x;run;
data _est_; merge _betax_ _sx_ _deltax_ _dsx_; run;
proc datasets nolist; delete _sx_ _betax_ _deltax_ _dsx_ _phout0_ _phout_ _logit0_ _logitout_ _miss1new_ _miss0new_ _new_;quit;  
%mend impute;


/*repeat the macro impute throughout the extreme cases*/
%macro iter(m1=,/*number of missing in case group*/
                     m0=/*number of missing in control group*/);
data _miss1_ _miss0_ _ob_;
   set test0;
	if x=. and y=1 then output _miss1_;
	if x=. and y=0 then output _miss0_;
	else if x^=. then output _ob_;
	run;
%let i1=0;
%let i2=0;
%impute(im11=0,im12=1,im01=1,im02=0);
          data est; set _est_;run; 
 %do %while (&i2<=&m0);
    %do %while(&i1<=&m1);
     %impute(im11=0,im12=1,im01=1,im02=0);
          data est; set est _est_;run; 
    %let i1=%eval(&i1+1);
    %end;
     %let i1=0;
	 %let i2=%eval(&i2+1);
%end;
proc datasets nolist; delete _miss1_ _miss0_ _ob_;quit;
%mend iter;
%iter(m1=47,m0=18);

/*the plot of betax v.s. deltax*/
proc gplot data=est;
symbol i=none v=star;
plot betax*deltax;
run;

/*delete the extreme values and then plot*/
data est0; set est; if deltax>=-10 and deltax<=10; run;
proc gplot data=est0;
symbol i=none v=star;
plot betax*deltax;
run;