options  nocenter  ps = 65 ls = 75;

goptions device=png gsfname=graphout;

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

log2bili=log2(bilirubin);
/* nedenfor laves variable til lineære splines, baseret på Q1, median og Q3 */
bili25=(bilirubin-15)*(bilirubin>15);
bili50=(bilirubin-31.5)*(bilirubin>31.5);
bili75=(bilirubin-74)*(bilirubin>74);

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;

/* SPM 1 */
title 'SPM1: log-rank test';

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

title 'SPM1: Kaplan-Meier kurver';

ods graphics on / imagename="kaplan_meier";
proc phreg plots(overlay=row cl)=survival data=pbc;
  model months*dead(0)= ;
  strata azathiop;
run;
ods graphics off;

/* 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;
run;

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

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

proc sort data=tegn;
  by bilirubin;
run;

ods graphics / imagename="spline";
proc sgplot data=tegn;
series Y=spline X=bilirubin;
run;
ods graphics off;

/* 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;
  model months*dead(0)=azathiop log2bili / rl;
run;

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

ods graphics on / imagename="modelkontrol";
proc phreg plots(overlay=row cl)=cumhaz data=pbc;
  model months*dead(0)=log2bili;
  strata azathiop;
  baseline out=ud1 loglogs=logCumrate;
run;
ods graphics off;

ods graphics on / imagename="logCumrate";
proc sgplot data=ud1;
series Y=logCumrate X=months / group=azathiop;
run;
ods graphics off;

ods graphics on / imagename="modelkontrol_bili";
proc phreg plots(overlay=row cl)=cumhaz data=pbc;
  model months*dead(0)=log2bili azathiop;
  strata bilirubin(37,92);
  baseline out=ud2 loglogs=logCumrate;
run;
ods graphics off;

ods graphics on / imagename="logCumrate_bili";
proc sgplot data=ud2;
series Y=logCumrate X=months / group=bilirubin;
run;
ods graphics off;

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;
ods graphics off;

title 'SPM6: numerisk check, bilirubin';

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;
ods graphics off;