* Program survwtmacs.v3.sas ;
* Created by Nancy Cook (BWH, HMS, HSPH) ;
* Reclassification macros for WEIGHTED SURVIVAL data;
* Macro to compute net reclssification improvement (NRI) - SURV_NRIWT;
* Uses non-null variance;
* Macro for Infinite category version of NRI - SURV_CONTNRIWT;
* Macro to compute integrated discrimination improvement (IDI) - SURV_IDIWT;
* Note: Should bootstrap for better estimate of SE;
* Macro to compute H-L Calibration event for Reclassification - SURV_RECLASSWT;
* Note: the last is best for low censoring situations by time T;
* DETAIL = 2 for detailed printout, 1 for limited, 0 for none;

* Usage examples:
* %surv_nriwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy,4,0.05,0.10,0.20);
* %surv_contnriwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy);
* %surv_idiwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy);
* %surv_reclasswt(probs,1,wt,outxy,pyrs,10,pdx,pdxy, 4, 0.05, 0.10, 0.20);


****************************************************************;

%macro SURV_NRIWT(DSNAME,DETAIL,WT,EVENT,PYRS,T,PROB1,PROB2,NCAT,C1,C2,C3,C4,C5,C6,C7,C8,C9);
* Macro to compute Net Reclassification Index (NRI) of Pencina, Stat Med 2007; 
* Uses up to 10 categories with cutpoints c1-c9;
*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 2 for detailed printout, 1 for limited, 0 for none;
*  WT = sample weight for case-cohort analysis (indiv, not risk set);
*  EVENT = outcome variable (coded 0,1);
*  PYS = follow-up time;
*  T = time of predicted risk (eg., 10 years) = scalar;
*  PROB1 = probability for model 1;
*  PROB2 = probability for model 2;
*  NCAT = number of categories in classification;
*  C1-C9 = category cutpoints (should have ncat-1 cutpoints);
*** NOTE: data need to be bootstrapped to get correct SE;

title2 Weighted Survival NRI for &NCAT Categories for Event &EVENT at Time &T : &PROB2 vs &PROB1 ;
title3 "Need to Bootstrap for SE";

data nri1; set &dsname;
_id_=_N_;
* compute diffs in probs;
event=&EVENT;
prob1=&PROB1;
prob2=&PROB2;
diffp=prob2-prob1;
pyrs=&pyrs;
wt=&wt;
if prob1>. and prob2>.;
keep _id_ wt event pyrs prob1 prob2 ;
run;

data nri1; set nri1;
* compute categories - set to above cutpoints;
%do i=1 %to &ncat-1;
  cut&i=&&c&i;
%end;
cut&ncat=1.01;
    
array cut {10} cut1-cut10;
array pcat {10} pcat1-pcat10;
if .<prob1<cut1 then pcat1=1;
if .<prob2<cut1 then pcat2=1;
do i=2 to &ncat;
   if cut(i-1)<=prob1<cut(i) then pcat1=i;
   if cut(i-1)<=prob2<cut(i) then pcat2=i;
end;

if pcat2 ne . and pcat1 ne . then do;
  if pcat2>pcat1 then disc=1;
  else if pcat2=pcat1 then disc=0;
  else if pcat2<pcat1 then disc=-1;
end;  
run;

%if &detail=2 %then %do;
proc freq data=nri1;
tables event*disc ;
weight wt;
title4 "Crude Reclassification Ignoring Time";
run;
title4; run;
%end;

ods listing exclude all;
ods html exclude all;
ods output OneWayFreqs=freqs;
run;
proc freq data=nri1;
tables disc ;
weight wt;
run;
ods listing exclude none;
ods html exclude none;
run;

* Compute KM estimator for survival for those moving up or down;
data two;
_id_=0;
wt=1;
event=.; pyrs=&t; 
* Output for up and down;
disc=1; output;
disc=0; output;
disc=-1; output;
run;

data _allph_; set nri1 two;
run;

proc phreg data=_allph_(where=(disc ne .)) noprint;
   model pyrs*event(0) =  ;
   id _id_;
   strata disc;
   weight wt;
   output out=survdat survival=_survest_ ;
run;

data base; set survdat;
if _id_=0;
* survival estimate for referent person at time &t;
_basesurv_=_survest_;
inc=1-_basesurv_;
keep disc _basesurv_ inc;
run;

proc sort data=base; by disc;
proc sort data=freqs; by disc;

data calc ; merge freqs base; by disc;
keep disc _basesurv_ inc Frequency Percent;
run;

data nri; set calc end=eof;
retain n_up n_same n_down inc_up inc_same inc_down ecase_up ecase_same ecase_down
       econt_up econt_same econt_down 0;
if disc=1 then do; 
   n_up=Frequency; inc_up=inc; ecase_up=n_up*inc_up; econt_up=n_up-ecase_up; 
   end;
else if disc=0 then do;
   n_same=Frequency; inc_same=inc; ecase_same=n_same*inc_same; econt_same=n_same-ecase_same; 
   end;
else if disc=-1 then do; 
   n_down=Frequency; inc_down=inc; ecase_down=n_down*inc_down; econt_down=n_down-ecase_down;  
   end;

if eof then do;
  totn=sum(n_up,n_same,n_down);
  ecase=sum(ecase_up,ecase_same,ecase_down);
  econt=sum(econt_up,econt_same,econt_down);
  %if ecase=0 or econt=0 %then %do;
     nri=.; znri=.; p2nri=.;
     %goto exit;
  %end;
  pup_case=ecase_up/ecase;
  pdown_case=ecase_down/ecase;
  pup_cont=econt_up/econt;
  pdown_cont=econt_down/econt;
  ri_case=pup_case-pdown_case;
  ri_contr=pdown_cont-pup_cont;
  nri=ri_case+ri_contr;

  *** Null variance (use non-null);
  evri_case_null=(pup_case+pdown_case)/ecase;
  evri_contr_null=(pup_cont+pdown_cont)/econt;
  evnri_null=evri_case_null+evri_contr_null;
  *** Non-null variance - use this;
  evri_case=(pup_case+pdown_case-(pup_case-pdown_case)**2)/ecase;
  evri_contr=(pup_cont+pdown_cont-(pup_cont-pdown_cont)**2)/econt;
  evnri=evri_case+evri_contr;

  esenri=sqrt(evnri);
  enricil=nri-1.96*esenri;
  enriciu=nri+1.96*esenri;
  ezri_case=ri_case/sqrt(evri_case);
  ezri_contr=ri_contr/sqrt(evri_contr);
  eznri=nri/sqrt(evnri);
  ep2ri_case=2*(1-probnorm(abs(ezri_case)));
  ep2ri_contr=2*(1-probnorm(abs(ezri_contr)));
  ep2nri=2*(1-probnorm(abs(eznri)));
  output;
  keep totn ecase econt n_up n_same n_down inc_up inc_same inc_down 
       pup_case pdown_case pup_cont pdown_cont
       ri_case ri_contr nri evri_case evri_contr evnri esenri enricil enriciu
       ezri_case ezri_contr eznri ep2ri_case ep2ri_contr ep2nri
       ;
end;

data _all_; set nri ; 
* compute categories again - set to above cutpoints (to print);
%do i=1 %to &ncat-1;
  cut&i=&&c&i;
%end;
if &ncat<10 then do;
  %do i=&ncat %to 9;
    cut&i=.;
  %end;
end;
ncat=&ncat;
run;

%if &detail>0 %then %do;
proc print data=_all_;
var    ncat cut1-cut9 totn ecase econt n_up n_same n_down inc_up inc_same inc_down 
       pup_case pdown_case pup_cont pdown_cont
       ri_case ri_contr nri evri_case evri_contr evnri esenri enricil enriciu
       ezri_case ezri_contr eznri ep2ri_case ep2ri_contr ep2nri
       ;
run;
%end;
title2; run;

* Usage:
* %surv_nriwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy,4,0.05,0.10,0.20);
%exit: %mend SURV_NRIWT;

****************************************************************;

%macro SURV_CONTNRIWT(DSNAME,DETAIL,WT,EVENT,PYRS,T,PROB1,PROB2);
* Macro to compute NRI with any up or down (no categories); 
*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 2 for detailed printout, 1 for limited, 0 for none;
*  WT = sample weight for case-cohort analysis (indiv, not risk set);
*  EVENT = outcome variable (coded 0,1);
*  PYS = follow-up time;
*  T = time of predicted risk (eg., 10 years) = scalar;
*  PROB1 = probability for model 1;
*  PROB2 = probability for model 2;
*** NOTE: data need to be bootstrapped to get correct SE;

title2 Weighted Continuous NRI for Event &EVENT at Time &T : &PROB2 vs &PROB1 ;
title3 "Need to Bootstrap for SE";

data nri1; set &dsname;
_id_=_N_;
* compute diffs in probs;
event=&EVENT;
pyrs=&pyrs;
prob1=&PROB1;
prob2=&PROB2;
diffp=prob2-prob1;
wt=&wt;
if prob1>. and prob2>.;
keep _id_ wt event pyrs prob1 prob2 diffp;
run;

data nri1; set nri1;
if prob2>prob1 then disc=1;
  else if prob2=prob1 then disc=0;
  else if prob2<prob1 then disc=-1;
run;

*** Control output as needed ***;
%if &detail=2 %then %do;
proc means data=nri1 n sumwgt sum mean median stddev min max;
var prob1 prob2 diffp disc event;
weight wt;
run;
%end;

%if &detail=2 %then %do;
proc freq data=nri1;
tables event*disc ;
weight wt;
title4 "Crude Reclassification Ignoring Time";
run;
title4; run;
%end;

ods listing exclude all;
ods html exclude all;
ods output OneWayFreqs=freqs;
run;
proc freq data=nri1;
tables disc ;
weight wt;
run;
ods listing exclude none;
ods html exclude none;
run;

* Compute KM estimator for survival for those moving up or down;
data two;
_id_=0;
wt=1;
event=.; pyrs=&t; 
* Output for up and down;
disc=1; bd=1; output;
disc=0; bd=1; output;
disc=-1; bd=1; output;
* Output for overall;
disc=.; bd=.; output;
run;

data _all_; set nri1 two;
run;

proc phreg data=_all_(where=(disc=1 or disc=-1)) noprint;
   model pyrs*event(0) =  ;
   id _id_;
   strata disc;
   weight wt;
   output out=survdat survival=_survest_ ;
run;

proc phreg data=_all_(where=(bd ne 1)) noprint;
   model pyrs*event(0) =  ;
   id _id_;
   weight wt;
   output out=survdat0 survival=survall ;
run;

data base; set survdat;
if _id_=0;
* survival estimate for referent person at time &t;
_basesurv_=_survest_;
inc=1-_basesurv_;
keep disc _basesurv_ inc;
run;

data base0; set survdat0;
* overall survival estimate;
if _id_=0;
* survival estimate for referent person at time &t;
inc0=1-survall;
disc=.;
keep disc survall inc0;
run;

proc sort data=base; by disc;
proc sort data=freqs; by disc;

data calc ; merge freqs base base0; by disc;
keep disc _basesurv_ inc survall inc0 Frequency Percent;
run;

data _all_; set calc end=eof;
retain n_up n_down totn inc_up inc_down inc_all ecase_up ecase_down
       econt_up econt_down 0;
if disc=1 then do; 
   n_up=Frequency; inc_up=inc; ecase_up=n_up*inc_up; econt_up=n_up-ecase_up; 
   totn=totn+Frequency;
   end;
else if disc=0 then do;
   n_same=Frequency; 
   totn=totn+Frequency;
   end;
else if disc=-1 then do; 
   n_down=Frequency; inc_down=inc; ecase_down=n_down*inc_down; econt_down=n_down-ecase_down;  
   totn=totn+Frequency;
   end;
else if disc=. then inc_all=inc0;

if eof then do;
  ecase=totn*inc_all;
  econt=totn-ecase;
  %if ecase=0 or econt=0 %then %do;
     nri=.; znri=.; p2nri=.;
     %goto exit;
  %end;
  pup_case=ecase_up/ecase;
  pdown_case=ecase_down/ecase;
  pup_cont=econt_up/econt;
  pdown_cont=econt_down/econt;
  ri_case=pup_case-pdown_case;
  ri_contr=pdown_cont-pup_cont;
  nri=ri_case+ri_contr;
  * corrected variance to reflect binomial;
  evri_case=pup_case*(1-pup_case)/ecase;
  evri_contr=pup_cont*(1-pup_cont)/econt;
  evnri=4*(evri_case+evri_contr);
  esenri=sqrt(evnri);
  enricil=nri-1.96*esenri;
  enriciu=nri+1.96*esenri;
  ezri_case=ri_case/sqrt(evri_case);
  ezri_contr=ri_contr/sqrt(evri_contr);
  eznri=nri/sqrt(evnri);
  ep2ri_case=2*(1-probnorm(abs(ezri_case)));
  ep2ri_contr=2*(1-probnorm(abs(ezri_contr)));
  ep2nri=2*(1-probnorm(abs(eznri)));
  output;
  keep totn ecase econt inc_all n_up n_same n_down inc_up inc_down
       pup_case pdown_case pup_cont pdown_cont
       ri_case ri_contr nri evri_case evri_contr evnri esenri enricil enriciu
       ezri_case ezri_contr eznri ep2ri_case ep2ri_contr ep2nri
       ;
end;

%if &detail>0 %then %do;
proc print data=_all_;
var    totn ecase econt inc_all n_up n_same n_down inc_up inc_down
       pup_case pdown_case pup_cont pdown_cont
       ri_case ri_contr nri evri_case evri_contr evnri esenri enricil enriciu
       ezri_case ezri_contr eznri ep2ri_case ep2ri_contr ep2nri
       ;
run;
%end;
title2; run;
* Usage:
* %surv_contnriwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy);
%exit: %mend SURV_CONTNRIWT;

****************************************************************;

%macro SURV_IDIWT(DSNAME,DETAIL,WT,EVENT,PYRS,T,PROB1,PROB2);
*  Macro to compute difference in Yates slopes or
     integrated discrimination improvement (IDI) from Pencina, 2007;
*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 1 or 2 for limited printout, 0 for none;
*  WT = sample weight for case-cohort analysis (indiv, not risk set);
*  EVENT = outcome variable (coded 0,1) (if 1,2 alter signs);
*  PYS = follow-up time;
*  T = time of predicted risk (eg., 10 years) = scalar;
*  PROB1 = probability for model 1;
*  PROB2 = probability for model 2;
*** NOTE: data need to be bootstrapped to get correct SE;

title2 Weighted Difference in Yates Slope (IDI) for Event &EVENT at Time &T : &PROB2 vs &PROB1 ;
title3 "Need to Bootstrap for SE";

data _ididat_; set &dsname;
_id_=_N_;
diffprob=&prob2-&prob1;
prob1=&prob1;
prob2=&prob2;
prob1sq=prob1*prob1;
prob2sq=prob2*prob2;
event=&EVENT;
pyrs=&PYRS;
wt=&wt;
if prob1>. and prob2>.;
run;

%if &detail<2 %then %do;
ods listing exclude all;
ods html exclude all;
%end;
proc ttest data=_ididat_;
class event;
var prob1 prob2 diffprob;
weight wt;
title4 "Crude Difference in Probs by Any Event - Ignoring Time";
run;
title4; run;

proc means data=_ididat_;
var prob1 prob2 prob1sq prob2sq;
weight wt;
output out=idistat mean=meanp1 meanp2 meansqp1 meansqp2 ;
title4 "Mean Predicted Probs";
run;
title4; run;

* Compute KM estimator for survival for those moving up or down;
data two;
_id_=0;
wt=1;
event=.; pyrs=&t; 
run;

data _all_; set _ididat_ two;
run;

proc phreg data=_all_ noprint;
   model pyrs*event(0) =  ;
   id _id_;
   weight wt;
   output out=survdat survival=_survest_ ;
run;

data base; set survdat;
if _id_=0;
* survival estimate for referent person at time &t;
_basesurv_=_survest_;
inc=1-_basesurv_;
keep _basesurv_ inc;
run;

data _idi_; merge idistat base;
* Note: compute population var from mean and squared probs;
varp1=meansqp1-(meanp1)**2;
varp2=meansqp2-(meanp2)**2;
yates1=varp1/(inc*(1-inc));
yates2=varp2/(inc*(1-inc));
idi=yates2-yates1;
rel_idi=yates2/yates1 - 1;
drop _TYPE_ _FREQ_;
run;

ods listing exclude none;
ods html exclude none;
%if &detail>0 %then %do;
proc print data=_idi_;
run;
%end;
title2; run;

* Usage:
* %survidiwt(probs,1,wt,outxy,pyrs,10,pdx,pdxy);
%mend SURV_IDIWT;

****************************************************************;

%macro SURV_RECLASSWT(DSNAME,DETAIL,WT,EVENT,PYRS,T,PROB1,PROB2,NCAT,C1,C2,C3,C4,C5,C6,C7,C8,C9);
*  Macro to compute WEIGHTED calibration eventistics for SURVIVAL DATA;
*  KxK (K=NCAT) table formed using categs of predicted probs (eg 0, 5, 10, 20%);
*  Counts categs for DF and computes for n>=20;
*  Allows up to 10 categories (usually 3 or 4);

*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 2 for detailed printout, 1 for limited, 0 for none;
*  WT = sampling weights;
*  EVENT = outcome variable (coded 0,1);
*  PYRS = person-time for Cox model;
*  T = desired follow-up time for probability estimate;
*  PROB1 = probability for model 1 (at time t);
*  PROB2 = probability for model 2 (at time t);
*  NCAT = number of categories in classification;
*  C1-C9 = category cutpoints (should have ncat-1 cutpoints);
*  Reclassification calibration test for cells with n>=20;
*  Also computes the percent reclassified correctly;


title2 Weighted Survival Reclassification for &NCAT Categories with CutPoints = %do j = 1 %to &ncat; &&c&j %end;;
title3 Event &EVENT at Time &T : &PROB2 vs &PROB1 ;

   data calc; set &DSNAME;
   _id_=_N_;
   event=&EVENT;
   prob1=&PROB1;
   prob2=&PROB2;
   wt=&wt;
   if prob1>. and prob2>.;

   * compute product terms for J eventistic;
   var1=prob1*(1-prob1);
   var2=prob2*(1-prob2);
   k=&ncat;

   * compute categories - set to above cutpoints;
   %do i=1 %to &ncat-1;
     cut&i=&&c&i;
   %end;
   cut&ncat=1.01;
      
   array cut {10} cut1-cut10;
   if .<prob1<cut1 then pcat1=1;
   if .<prob2<cut1 then pcat2=1;
   do i=2 to &ncat;
     if cut(i-1)<=prob1<cut(i) then pcat1=i;
     if cut(i-1)<=prob2<cut(i) then pcat2=i;
   end;

   * look at number of categories moved;
   diffcat=pcat2-pcat1;
   * truncate diffcat;
   if diffcat>2 then diffcat2=2;
      else if .<diffcat<-2 then diffcat2=-2;
      else diffcat2=diffcat;
   run;


   *** check output needed;
   %if &detail=2 %then %do;
   proc freq data=calc;
   tables pcat1*pcat2;
   weight wt;
   title3 "Reclassification table for &PROB1 and &PROB2";
   run;
   %end;


* Compute KM estimator for survival within strata;
data two;
_id_=0;
event=.; &pyrs=&t; 
* Set strata variables in categories;
pcat1=0;
pcat2=0;
do i=1 to &ncat;
  pcat1=pcat1+1;
  do j=1 to &ncat;
     pcat2=pcat2+1;
     output;
  end;
  pcat2=0;
end;
run;

data _all_; set calc two;
proc phreg data=_all_ noprint;
   model &pyrs*event(0) =  ;
   id _id_;
   strata pcat1 pcat2;
   weight wt;
   output out=survdat survival=_survest_ ;
run;

data base; set survdat;
if _id_=0;
* survival estimate for referent person;
_basesurv_=_survest_;
keep pcat1 pcat2 _basesurv_ ;
run;

proc sort data=base; by pcat1 pcat2;
proc sort data=calc; by pcat1 pcat2;

data calc ; merge calc base; by pcat1 pcat2;
if _id_ ne 0;
surv=_basesurv_;
obsp=1-surv;
run;
 
   ods listing exclude all;
   ods html exclude all;
   ods output CrossTabFreqs=classfreqs;
   proc freq data=calc;
   tables pcat1*diffcat2;
   weight wt;
   title3 "Reclassification table for &PROB1 and &PROB2";
   run;

   proc sort data=classfreqs; by pcat1;
   data reclass; set classfreqs; by pcat1;
   retain pctup1 pctup2 pctd1 pctd2 pctsame;
   if first.pcat1 then do;
      pctup1=0; pctup2=0; pctd1=0; pctd2=0; pctsame=0;
   end;
   if diffcat2=0 then pctsame=RowPercent;
     else if diffcat2=1 then pctup1=RowPercent;
     else if diffcat2=2 then pctup2=RowPercent;
     else if diffcat2=-1 then pctd1=RowPercent;
     else if diffcat2=-2 then pctd2=RowPercent;
   if pcat1=. then do;
     if diffcat2=0 then pctsame=Percent;
     else if diffcat2=1 then pctup1=Percent;
     else if diffcat2=2 then pctup2=Percent;
     else if diffcat2=-1 then pctd1=Percent;
     else if diffcat2=-2 then pctd2=Percent;
   end;
   if last.pcat1 then do;
     pctup=pctup1+pctup2;
     pctdown=pctd1+pctd2;
     pctreclass=100-pctsame;
     * compute categories again - set to above cutpoints;
     ncat=&ncat;
     %do i=1 %to &ncat-1;
         cut&i=&&c&i;
     %end;
     if &ncat<10 then do;
       %do i=&ncat %to 9;
         cut&i=.;
       %end;
     end;
     keep ncat cut1-cut9 pcat1 pctup1 pctup2 pctd1 pctd2 pctup pctdown pctsame pctreclass;
     output;
   end;  
   run;
   ods listing exclude none;
   ods html exclude none;
   run;
   %if &detail>0 %then %do;
   proc print data=reclass;
   run;
   %end;
      
   %if &detail=2 %then %do;
    proc freq data=calc;
    tables pcat1*pcat2*obsp/ list;
    weight wt;
    title3 "Calibration table for &PROB1 and &PROB2";
    run;
   %end;
   
   ods listing exclude all;
   ods html exclude all;
   ods output Summary=calibs;
   proc means data=calc n sumwgt sum mean median min max;
     * Eliminate missing from all variables;
     where obsp ne . and prob1 ne . and prob2 ne .;
     class pcat1 pcat2;
     types pcat1*pcat2;
     var prob1 prob2 var1 var2 obsp k; 
     weight wt;
   run;
   
   ods listing exclude none;
   ods html exclude none;
   run;
   %if &detail>0 %then %do;
   proc print data=calibs;
   var pcat1 pcat2 obsp_N obsp_SumWgt obsp_Sum prob1_Mean prob2_Mean obsp_Mean
       var1_Sum var2_Sum;
   title3 "Estimated Risks";
   run;
   %end;
   

   data caltermsb calstatsb(keep= ncat maxcat dfmax nccat df chisq1 pcalib1
            chisq_adj1 pcal_adj1 chisq2 pcalib2 chisq_adj2 pcal_adj2 
            dfj jsq1 pj1 jsq_adj1 pj_adj1 jsq2 pj2 jsq_adj2 pj_adj2 
            nreclass ncorr pctcorr) ; 
     set calibs end=eof;
     retain chisq1 chisq_adj1 chisq2 chisq_adj2
            jsq1 jsq_adj1 jsq2 jsq_adj2 nccat nreclass ncorr 0;
     * only save cells with count at least 20;
    if obsp_N >=20 then do;
     num1=(obsp_Sum - prob1_Sum)**2;
     denom1=prob1_Sum*(1-prob1_Mean);
     denom12=(prob1_Sum+1)*(1-prob1_Mean+1/prob1_SumWgt);
     chisq1=chisq1 + num1/denom1;
     chisq_adj1=chisq_adj1 + num1/denom12;
     phi1=var1_Sum/denom1;
     jsq1=jsq1 + num1/(phi1*denom1);
     * Note: evaluate adjusted Jsq, but do not print;
     jsq_adj1=jsq_adj1 + num1/(phi1*denom12);
     num2=(obsp_Sum - prob2_Sum)**2;
     denom2=prob2_Sum*(1-prob2_Mean);
     denom22=(prob2_Sum+1)*(1-prob2_Mean+1/prob2_SumWgt);
     chisq2=chisq2 + num2/denom2;
     chisq_adj2=chisq_adj2 + num2/denom22;
     phi2=var2_Sum/denom2;
     jsq2=jsq2 + num2/(phi2*denom2);
     jsq_adj2=jsq_adj2 + num2/(phi2*denom22);
     nccat=nccat+1;

     OEratio1=obsp_Sum/prob1_Sum;
     OEratio2=obsp_Sum/prob2_Sum;
     if obsp_Sum>0 then do;
       sdo=sqrt(1/obsp_Sum);
       OEcil1=OEratio1*exp(-1.96*sdo);
       OEciu1=OEratio1*exp(1.96*sdo);
       OEcil2=OEratio2*exp(-1.96*sdo);
       OEciu2=OEratio2*exp(1.96*sdo);
     end;

     * compute those reclass correctly;
     if pcat2 ne pcat1 then nreclass=nreclass+obsp_SumWgt;
     %do i=1 %to &ncat-1;
        cut&i=&&c&i;
        if pcat2>pcat1 and pcat2=&i+1 then do;
            * put pcat2= pcat1= cut&i= obsp_Sum= obsp_Mean= obsp_N= obsp_SumWgt= ;
            if obsp_Mean>=cut&i then ncorr=ncorr+obsp_SumWgt;
        end;
        else if pcat2<pcat1 and pcat2=&i then do;
            * put pcat2= pcat1= cut&i= obsp_Sum= obsp_Mean= obsp_N= obsp_SumWgt= ;
            if obsp_Mean<cut&i then ncorr=ncorr+obsp_SumWgt;
        end;
     %end;
     
     output caltermsb;
    end;
     if eof then do;
       dfmax=k_Mean*k_Mean-2;
       df=nccat-2;
       pcalib1=1-probchi(chisq1,df);
       pcal_adj1=1-probchi(chisq_adj1,df);
       pcalib2=1-probchi(chisq2,df);
       pcal_adj2=1-probchi(chisq_adj2,df);
       dfj=nccat-1;
       pj1=1-probchi(jsq1,dfj);
       pj_adj1=1-probchi(jsq_adj1,dfj);
       pj2=1-probchi(jsq2,dfj);
       pj_adj2=1-probchi(jsq_adj2,dfj);
       ncat=&ncat;
       maxcat=ncat**2;
       pctcorr=ncorr/nreclass;
       output calstatsb;
     end;

   %if &detail=2 %then %do;
   proc print data=caltermsb;
   var pcat1 pcat2 num1 denom1 denom12 phi1 num2 denom2 denom22 phi2
       nreclass ncorr OEratio1 OEcil1 OEciu1 OEratio2 OEcil2 OEciu2 ;
   title3 "Weighted Calibration Terms for &PROB1 and &PROB2 in Predicting &EVENT";
   title4 "Using KM Incidence Estimates and Cross-Classified Cells with N >= 20";
   run;
   %end;
   %if &detail>0 %then %do;
   proc print data=calstatsb;
   var ncat maxcat dfmax nccat df chisq1 pcalib1 chisq_adj1 pcal_adj1
       chisq2 pcalib2 chisq_adj2 pcal_adj2 dfj jsq1 pj1 jsq2 pj2
       nreclass ncorr pctcorr;
   title3 "Weighted Calibration Statistics for &PROB1 and &PROB2 in Predicting &EVENT";
   title4 "Using KM Incidence Estimates and Cross-Classified Cells with N >= 20";
   run;
   %end;
   title2; 
   run;

* Usage:
* %surv_reclasswt(probs,1,wt,outxy,pyrs,10,pdx,pdxy, 4, 0.05, 0.10, 0.20);
%mend SURV_RECLASSWT;

****************************************************************;