/*                   Video 1: Slide 1-10                            */

* slide 4: Data;
data a1;
	infile "http://publicifsv.sund.ku.dk/~sr/BasicStatistics/datasets/veteran.txt" URL firstobs=2;
	input id celltype $ time status karno diagtime $ age prior treat;
     
    log2time=log2(time);
run;	

* slide 7: Boxplots mm;
proc sgplot data=a1;
	vbox time / category=treat;
run;
proc sgplot data=a1;
	vbox log2time / category=treat;
run;
proc ttest data=a1;
	var log2time;
	class treat;
run;
proc npar1way data=a1;
	var time;
	class treat;
run;	


/*                   Video 2: Slide 11-20                           */

* Slide 19-20: KM kurver;
proc phreg plots(overlay=row cl)=survival data=a1;
	model time*status(0)= ;
	strata treat;
run;



/*                   Video 3: Slide 21-30                           */

* Slide 21: Kumulativ incidens;
proc phreg data=a1;
	model time*status(0)= ;
	strata treat;
	baseline out=ud survival=_ALL_;
run;

* Slide 24-25: Log-rank test;
	* Via phreg - p-vaerdi aflaeses under Score test;
	proc phreg data=a1;
		class treat(ref="0") / param=glm;
		model time*status(0)=treat / ties=discrete rl;
	run;
* Via lifetest - p-vaerdi aflaeses under Log-Rank;
proc lifetest data=a1;
	time time*status(0);
	strata treat;
run;



/*                   Video 4: Slide 28-36                           */

* Slide 31: Kumulerede hazards;
proc phreg data=a1 plots(overlay=row cl)=cumhaz;
	model time*status(0)= ;
	strata treat;
run;

* Slide 33: Estimat af behandlingseffekt;
proc phreg data=a1;
	class treat (ref="0") / param=glm;
	model time*status(0)= treat / rl ties=discrete;
run;	


/*                   Video 5: Slide 37-42                           */

* Slide 37;
proc phreg plots(overlay=row cl)=(survival cumhaz) data=a1;
	model time*status(0)= ;
	strata treat;	
	baseline out=ud loglogs=logCumrate;
run;
data ud;
	set ud;
	log2time=log2(time);
run;	
proc sgplot data=ud;
	series Y=logCumrate X=log2time / group=treat;
run;

* Slide 40-41;
proc phreg data=a1;
	class treat(ref="0") / param=glm;
	model time*status(0)= treat / ties=discrete rl;
	assess ph / resample seed=106165;
run;

/*                   Video 6: Slide 42-48                           */

* Ingen kode *;

/*                   Video 7: Slide 50-57                           */

* Slide 50;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)=treat celltype karno / rl;
	estimate "effekt af 10 karno" karno 10 / exp cl;
	hazardratio karno / unit=10;
	hazardratio celltype / diff=all;
	assess var=(karno) ph / resample seed=106165;
	baseline out=ud loglogs=logCumrate survival=_ALL_;
run;


/*                   Video 8: Slide 58-67                           */


* Slide 59;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)= celltype karno / rl;
	strata treat;
	baseline out=ud loglogs=logCumrate;
run;
data ud;
	set ud;
	log2time=log2(time);
run;	
proc sgplot data=ud;
	series Y=logCumrate X=log2time / group=treat;
run;

* Slide 60;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)= treat karno / rl;
	strata celltype;
	baseline out=ud loglogs=logCumrate;
run;
data ud;
	set ud;
	log2time=log2(time);
run;	
proc sgplot data=ud;
	series Y=logCumrate X=log2time / group=celltype;
run;

* Slide 61: PH for alle variable;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)= treat celltype karno / rl;
	assess ph / resample seed=106165;
run;


/*                   Video 9: Slide 65-75                           */

* Slide 66: Linearitetscheck karno;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)= treat celltype karno / rl;
	assess var=(karno) / resample seed=106165;
run;

* Slide 67;
data a1;
	set a1;
	karno40=(karno>40)*(karno-40);
	karno80=(karno>80)*(karno-80);
run;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)=karno karno40 karno80 treat celltype / rl;
	output out=tegn xbeta=spline;
	Testlinearitet: TEST karno40=karno80=0;
run;
proc sort data=tegn; by karno; run;
proc sgplot data=tegn;
	series x=karno y=spline / group=celltype;
    where treat=0;
run;	
proc sgplot data=tegn;
	series x=karno y=spline / group=celltype;
    where treat=1;
run;	

* Slide 72: Inddeling af karno;
data a1;
	set a1;
	karnogrp="under 40";
	if karno>40 then karnogrp="40-60";
	if karno>60 then karnogrp="60-75";
	if karno>75 then karnogrp="over 75";
run;
proc phreg data=a1;
  class treat(ref="0") karnogrp(ref="over 75") celltype / param=glm;  
  model time*status(0)=treat celltype / rl;
  strata karnogrp;
  baseline out=ud loglogs=logCumrate;
run;
data ud;
	set ud;
	log2time=log2(time);
run;	
proc sgplot data=ud;
	series Y=logCumrate X=log2time / group=karnogrp;
run;


* Slide 73;
data a1;
 	set a1;
 	karno10=karno/10;
run; 	
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)=treat celltype karno1 karno2 karno3 karno4 / rl;
		if time<=23 then karno1=karno10; else karno1=0;
		if 23<time<=62 then karno2=karno10; else karno2=0;
		if 62<time<=147.5 then karno3=karno10; else karno3=0;
		if 147.5<time then karno4=karno10; else karno4=0;
	NumCheck: test karno1=karno2=karno3=karno4;
run;

/*                   Video 10: Slide 75-81                           */


* Slide 76: Interaktion;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)=treat celltype treat*celltype karno1 karno2 karno3 karno4 / rl;
		if time<=23 then karno1=karno10; else karno1=0;
		if 23<time<=62 then karno2=karno10; else karno2=0;
		if 62<time<=147.5 then karno3=karno10; else karno3=0;
		if 147.5<time then karno4=karno10; else karno4=0;
	NumCheck: test karno1=karno2=karno3=karno4;
run;
	
/*                   Video 10: Slide 76-slut                         */
	
* Slide 76: Interaktion;
proc phreg data=a1;
	class treat(ref="0") celltype / param=glm;
	model time*status(0)=celltype treat*celltype karno1 karno2 karno3 karno4 / rl;
		if time<=23 then karno1=karno10; else karno1=0;
		if 23<time<=62 then karno2=karno10; else karno2=0;
		if 62<time<=147.5 then karno3=karno10; else karno3=0;
		if 147.5<time then karno4=karno10; else karno4=0;
	NumCheck: test karno1=karno2=karno3=karno4;
	* EKSTRA: Hvis vi oensker at bestemme HR for hver behandlingsgruppe;
	hazardratio celltype / at (treat="0" "1");
run;