data pbc;
infile "http://staff.pubhealth.ku.dk/~lts/basal/data/pbc.txt"
       url firstobs=2;
input dead months azathiop bilirubin;

log2bili=log2(bilirubin);

bili25=(bilirubin-30)*(bilirubin>30);
bili50=(bilirubin-54)*(bilirubin>54);
bili75=(bilirubin-125)*(bilirubin>125);

LABEL
    Dead = "Død (0=censur, 1=død)"
    Months = "Måneder fra randomisering"
    Azathiop = "Azathioprine (0~Placebo, 1~Azathioprine)"
    Bilirubin = "Baseline serum bilirubin"
  ;     
RUN;

ods graphics on;

/* SPM 1 */
title 'SPM1: Kaplan-Meier kurver';

proc phreg plots(overlay=row cl)=survival data=pbc;
  model months*dead(0)= ;
  strata azathiop;
run;

title 'SPM1: log-rank test';

proc phreg data=pbc;
  class azathiop(ref="0") / param=glm;  
  model months*dead(0)=azathiop / ties=discrete rl;
estimate "placebo vs azathiop" azathiop -1 1 / exp cl;
estimate "azathiop vs placebo" azathiop 1 -1 / exp cl;
run;

/* SPM 2 */
title 'SPM2: Serum Bilirubin';

proc sgplot data=pbc;
histogram bilirubin;
run;

proc sgplot data=pbc;
vbox bilirubin / group=azathiop;
run;

proc sgplot data=pbc;
histogram log2bili;
run;

proc sgplot data=pbc;
vbox log2bili / group=azathiop;
run;

proc phreg data=pbc;
  model months*dead(0)=bilirubin log2bili;
  Begge: TEST bilirubin=log2bili=0;
run;

proc phreg data=pbc;
  class azathiop(ref="0") / param=glm;  
  model months*dead(0)=azathiop bilirubin log2bili;
  Begge: TEST bilirubin=log2bili=0;
run;

proc means N Q1 median Q3 data=pbc; 
var bilirubin;
run;

proc means N Q1 median Q3 data=pbc; where dead=1;
var bilirubin;
run;

proc phreg data=pbc;
  model months*dead(0)=bilirubin bili25 bili50 bili75;
  output out=tegn1 xbeta=spline;
  Testlinearitet: TEST bili25=bili50=bili75=0;
run;

proc sort data=tegn1;
  by bilirubin;
run;

proc sgplot data=tegn1;
series Y=spline X=bilirubin / markers;
run;

proc phreg data=pbc;
  class azathiop(ref="0") / param=glm;  
  model months*dead(0)=azathiop bilirubin bili25 bili50 bili75;
  output out=tegn xbeta=spline;
  Testlinearitet: TEST bili25=bili50=bili75=0;
run;

proc sort data=tegn;
  by azathiop bilirubin;
run;

proc sgplot data=tegn;
series Y=spline X=bilirubin / group=azathiop markers;
run;

/* SPM 3 */
title 'SPM3: Serum Bilirubin';

proc phreg data=pbc;
  model months*dead(0)=log2bili / rl;
run;

/* SPM 4 */
title 'SPM3:Azathiop, med Serum Bilirubin';

proc phreg data=pbc;
  class azathiop(ref="0") / param=glm;  
  model months*dead(0)=azathiop log2bili / rl;
assess var=(log2bili) ph / resample seed=106165;
run;

/* SPM 6 */
title 'SPM6: log kumulerede rater';

proc phreg plots(overlay=row cl)=cumhaz data=pbc;
  model months*dead(0)=log2bili;
  strata azathiop;
  baseline out=ud1 loglogs=logCumrate;
run;

proc sgplot data=ud1;
series Y=logCumrate X=months / group=azathiop;
run;

proc univariate data=pbc; where dead=1;
    var bilirubin;
    output out=regn pctlpre=P_ pctlpts=33.33, 66.67;
run;

proc print data=regn;
run;

proc phreg plots(overlay=row cl)=cumhaz data=pbc;
  model months*dead(0)=azathiop;
  strata bilirubin(37,92);
  baseline out=ud2 loglogs=logCumrate;
run;

proc sgplot data=ud2;
series Y=logCumrate X=months / group=bilirubin;
run;

title 'SPM6: numerisk check, azathiop';

proc phreg data=pbc;
  model months*dead(0)=log2bili az1 az2 az3 / rl;
  if    months<=12 then az1=azathiop; else az1=0;
  if 12<months<=56 then az2=azathiop; else az2=0;
  if 56<months     then az3=azathiop; else az3=0;
NumCheckAza: test az1=az2=az3;
run;

title 'SPM6: numerisk check, bilirubin';

proc univariate data=pbc; where dead=1;
    var months;
    output out=regn pctlpre=P_ pctlpts=33.33, 66.67;
run;

proc print data=regn;
run;

proc phreg data=pbc;
  model months*dead(0)=azathiop log2bili1 log2bili2 log2bili3 / rl;
  if    months<=19 then log2bili1=log2bili; else log2bili1=0;
  if 19<months<=50 then log2bili2=log2bili; else log2bili2=0;
  if 50<months     then log2bili3=log2bili; else log2bili3=0;
NumCheckBili: test log2bili1=log2bili2=log2bili3;
run;