/****************************************************************/
/* SAS SAMPLE LIBRARY */
/* */
/* NAME: PHR610EX */
/* TITLE: Applications of the Counting Process Formulation */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: survival analysis, */
/* PROCS: PHREG IML */
/* DATA: */
/* */
/* REF: SAS/STAT Software: Changes and Enhancements, */
/* Release 6.10 (PROC PHREG chapter) */
/* MISC: SAS/GRAPH part of the code is commented out. */
/* Remove comments and add device= option to use */
/* the PROC GPLOT examples */
/* */
/****************************************************************/
/*--------------------------------------------------------------
Example 1: Multiple Failure Outcomes
-------------------------------------------------------------*/
/*
The cancer bladder data listed in Wei, Lin, and Weissfeld (1989,
JASA 84, 1065-71) are used in fitting the following models:
1) Anderson-Gill (AG) multiplicative hazards model
2) Wei, Lin, and Weissfeld (WLW) marginal model
3) Prentice, Williams, and Peterson (PWP) conditional models
The data consist of 86 patients with superficial bladder tumors,
which were removed when they entered the study. Forty eight
of these patients were randomized into the placebo group,
and 38 into the thiotepa group. Many patients have multiple
recurrences of tumors in the study, and new tumors were removed
at each visit. The data set contains the first four recurrences
of the tumor for each patient, and each recurrence time was
measured from the patient's entry time into the study.
The input data consist of eight variables:
- TRT (treatment group, where 1=placebo and 2=thiotepa)
- TIME (followup time)
- NUMBER (number of initial tumors)
- SIZE (initial tumor size)
- T1, T2, T3, and T4 (times of the four possible recurrences of
the bladder tumor. A patient with only two recurrences has a
blank value in T3 and T4.
In the data set BLADDER, four observations are created for each
patient, each corresponding to one of the four tumor recurrences.
In addition to values of TRT, NUMBER, and SIZE for the patient,
each observation contains the following variables:
- ID (patient's identification, which is the sequence number of
the input data)
- VISIT (event number, where 1=first recurrence, 2= second
recurrence, and so on)
- TSTART (time of the $(k-1)$st recurrence if VISIT=$k$, or the
the entry time if VISIT=1)
- TSTOP (time of the $k$th recurrence if VISIT=$k$)
- STATUS (event status, where 1=recurrence, 0=censored).
For instance, a patient with only one recurrence time at month
6 who was followed until month 10 will have values for VISIT,
TSTART, TSTOP, and STATUS of (1,0,6,1), (2,6,10,0), (3,10,10,0),
(4,10,10,0). If the followup time of a patient is beyond the time
of the fourth tumor recurrence, an extra observation is created
to represent this last followup period. For instance, a patient
with recurrences at months 2, 15, 34, and 50 who was followed
until month 61 will have five observations whose values for
VISIT, TSTART, TSTOP, and STATUS are (1,0,2,1), (2,2,15,1),
(3,15,34,1), (4,34,50,1), and (5,50,61,0). In the former
situation, the last two observations are redundant for the AG
model but they are important for the WLW analysis. In the
latter situation, the fifth observation is not needed for the
WLW model, but it is indispensible to the AG analysis.
*/
title 'Example 1: Multiple Failure Outcomes';
data bladder(keep=id tstart tstop status trt number size visit);
retain id tstart 0;
array tt t1-t4;
infile cards missover;
input trt time number size @27 t1 @31 t2 @35 t3 @39 t4;
id + 1;
tstart=0;
do over tt;
visit=_i_;
if tt = . then do;
tstop=time;
status=0;
end;
else do;
tstop=tt;
status=1;
end;
output;
tstart=tstop;
end;
if (tstart < time) then do;
tstop= time;
status=0;
visit=5;
output;
end;
cards;
1 0 1 1
1 1 1 3
1 4 2 1
1 7 1 1
1 10 5 1
1 10 4 1 6
1 14 1 1
1 18 1 1
1 18 1 3 5
1 18 1 1 12 16
1 23 3 3
1 23 1 3 10 15
1 23 1 1 3 16 23
1 23 3 1 3 9 21
1 24 2 3 7 10 16 24
1 25 1 1 3 15 25
1 26 1 2
1 26 8 1 1
1 26 1 4 2 26
1 28 1 2 25
1 29 1 4
1 29 1 2
1 29 4 1
1 30 1 6 28 30
1 30 1 5 2 17 22
1 30 2 1 3 6 8 12
1 31 1 3 12 15 24
1 32 1 2
1 34 2 1
1 36 2 1
1 36 3 1 29
1 37 1 2
1 40 4 1 9 17 22 24
1 40 5 1 16 19 23 29
1 41 1 2
1 43 1 1 3
1 43 2 6 6
1 44 2 1 3 6 9
1 45 1 1 9 11 20 26
1 48 1 1 18
1 49 1 3
1 51 3 1 35
1 53 1 7 17
1 53 3 1 3 15 46 51
1 59 1 1
1 61 3 2 2 15 24 30
1 64 1 3 5 14 19 27
1 64 2 3 2 8 12 13
2 1 1 3
2 1 1 1
2 5 8 1 5
2 9 1 2
2 10 1 1
2 13 1 1
2 14 2 6 3
2 17 5 3 1 3 5 7
2 18 5 1
2 18 1 3 17
2 19 5 1 2
2 21 1 1 17 19
2 22 1 1
2 25 1 3
2 25 1 5
2 25 1 1
2 26 1 1 6 12 13
2 27 1 1 6
2 29 2 1 2
2 36 8 3 26 35
2 38 1 1
2 39 1 1 22 23 27 32
2 39 6 1 4 16 23 27
2 40 3 1 24 26 29 40
2 41 3 2
2 41 1 1
2 43 1 1 1 27
2 44 1 1
2 44 6 1 2 20 23 27
2 45 1 2
2 46 1 4 2
2 46 1 4
2 49 3 3
2 50 1 1
2 50 4 1 4 24 47
2 54 3 4
2 54 2 1 38
2 59 1 3
;
title2 'Andersen-Gill Multiplicative Hazards Model';
proc phreg data=bladder;
model (tstart, tstop) * status(0) = trt number size;
where tstart < tstop;
run;
* Fitting the marginal models of Wei, Lin and Weissfeld;
data bladder2;
set bladder;
if visit < 5;
trt1= trt * (visit=1);
trt2= trt * (visit=2);
trt3= trt * (visit=3);
trt4= trt * (visit=4);
number1= number * (visit=1);
number2= number * (visit=2);
number3= number * (visit=3);
number4= number * (visit=4);
size1= size * (visit=1);
size2= size * (visit=2);
size3= size * (visit=3);
size4= size * (visit=4);
run;
title2 'Marginal Proportional Hazards Models';
proc phreg data=bladder2 out=est1;
model tstop*status(0)=trt1-trt4 number1-number4 size1-size4;
output out=out1 dfbeta=dt1-dt4 / order=data;
strata visit;
id id;
run;
title2;
proc means data=out1 noprint;
by id;
var dt1-dt4;
output out=out2 sum=dt1-dt4;
proc iml;
use out2;
read all var{dt1 dt2 dt3 dt4} into x;
v=x` * x;
reset noname;
vname={"TRT1","TRT2","TRT3","TRT4"};
print,"Estimated Covariance Matrix",, v[colname=vname
rowname=vname format=10.5];
proc iml;
use est1;
read all var{trt1 trt2 trt3 trt4} into trt;
b= trt`;
use out2;
read all var{dt1 dt2 dt3 dt4} into x;
v=x` * x;
nparm= nrow(b);
se=sqrt(vecdiag(v));
reset noname;
stitle={"Estimate", " Std Error"};
vname={"TRT1","TRT2","TRT3","TRT4"};
tmpprt= b || se;
print,tmpprt[colname=stitle rowname=vname format=10.5];
print,"Estimated Covariance Matrix",,
v[colname=vname rowname=vname format=10.5];
/* H0: beta11=beta12=beta13=beta14=0 */
chisq= b` * inv(v) * b;
df= nrow(b);
p= 1-probchi(chisq,df);
print ,,"Testing H0: no treatment effects", ,
"Wald Chi-Square = " chisq, "DF = " df,
"p-value = "p[format=5.4],;
/* Assume beta11=beta12=beta13=beta14 and
estimate the common value */
c= {1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1};
cb= c * b;
si= c * v * t(c);
e= j(4,1,1);
isi=inv(si);
h= inv(e` * isi * e) * isi * e;
b1= t(h) * cb;
se= sqrt(t(h) * si * h);
zscore= b1 / se;
p= 1- probchi( zscore * zscore, 1);
print ,"Estimation of the Common Parameter for Treatment",,
"optimal weights = "h,
"Estimate = " b1,
"Standard Error = " se,
"z-score =" zscore,
"2-sided p-value = " p[format=5.4];
quit;
* Prentice-Williams-Peterson Conditional Models;
data bladder3(drop=lstatus);
retain lstatus;
set bladder2;
by id;
if first.id then lstatus=1;
if (status=0 and lstatus=0) then delete;
lstatus=status;
gaptime=tstop-tstart;
run;
title2 'PWP Total Time Model with Noncommon Effects';
proc phreg data=bladder3;
model tstop * status(0) = trt1-trt4 number1-number4 size1-size4;
strata visit;
run;
title2 'PWP Gap Time Model with Noncommon Effects';
proc phreg data=bladder3;
model gaptime * status(0) = trt1-trt4 number1-number4 size1-size4;
strata visit;
run;
title2 'PWP Total Time Model with Common Effects';
proc phreg data=bladder3;
model tstop * status(0) = trt number size;
strata visit;
run;
title2 'PWP Gap Time Model with Common Effects';
proc phreg data=bladder3;
model gaptime * status(0) = trt number size;;
strata visit;
run;
/*---------------------------------------------------------------
Example 2. Time-Dependent Repeated Measurements
---------------------------------------------------------------*/
/*
Consider an experiment to study the dosing effect of a tumor-
promoting agent. Forty-five rodents initially exposed to a
carcinogen were randomly assigned to three dose groups. After
the first death of an animal, the rodents were examined every
week for the number of papillomas. Investigators were interested
in determining the effects of dose on the carcinoma incidence
after adjusting for the number of papillomas.
The input data set TUMOR consists of 20 variables:
- ID (subject ID)
- TIME (survival time of the subject)
- DEAD (censoring status where 1=dead and 0=censored)
- DOSE (dose of the tumor-promoting agent)
- P1 - P16 (number of papillomas at the 16 times that animals
died. These 16 times are weeks 27, 33, 34, 37, 41, 43, 45, 46,
47, 49, 50, 51, 53, 65, 67, and 71. For instance, subject 1
died at week 47; it had no papilloma at week 27, one papilloma
at week 33, 5 papillomas at week 34, 6 at weeks 37, 8 at week
41, and 10 at weeks 43, 45, 46, 47. For an animal that died
before week 71, the number of papillomas is missing for those
times beyond its death.)
*/
title 'Example 2. Time-Dependent Repeated Measurements';
data tumor;
retain id 0;
infile cards missover;
input time dead dose p1-p16;
id + 1;
label id='Subject Number';
cards;
47 1 1.0 0 1 5 6 8 10 10 10 10
71 1 1.0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
81 0 1.0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
65 1 1.0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
69 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
67 1 1.0 0 0 0 1 1 2 2 2 2 3 3 3 3 3 3
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
37 1 1.0 9 9 9 9
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
77 0 1.0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
81 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
54 0 2.5 0 1 1 1 1 2 2 2 2 2 2 2 2
53 0 2.5 0 0 0 0 0 0 0 0 0 0 0 0 0
38 0 2.5 5 13 13 14
54 0 2.5 2 5 6 6 6 6 6 6 6 6 6 6 6
51 1 2.5 15 15 15 15 16 16 17 17 17 17 17 17
47 1 2.5 13 20 20 20 20 20 20 20 20
27 1 2.5 22
41 1 2.5 6 13 13 13 13
49 1 2.5 0 3 3 3 3 3 3 3 3 3
53 0 2.5 0 0 0 1 1 1 1 1 1 1 1 1 1
50 1 2.5 0 0 0 2 3 4 6 6 6 6 6
37 1 2.5 3 14 15 15
49 1 2.5 2 3 3 3 3 3 4 4 4 4
46 1 2.5 4 6 6 7 9 9 9 9
48 0 2.5 15 26 26 26 26 26 26 26 26
54 0 10.0 12 13 14 15 15 15 15 15 15 15 15 15 15
37 1 10.0 12 16 16 17
53 1 10.0 3 6 6 6 6 6 6 6 6 6 6 6 6
45 1 10.0 4 10 12 15 20 20 20
53 0 10.0 6 8 10 13 13 13 15 15 15 15 15 15 20
49 1 10.0 0 1 2 2 2 2 2 2 2 2
39 0 10.0 7 8 8 8
27 1 10.0 17
49 1 10.0 0 2 6 9 14 14 14 14 14 14
43 1 10.0 14 18 18 20 20 20
28 0 10.0 8
34 1 10.0 11 18 18
45 1 10.0 10 12 12 16 16 16 16
37 1 10.0 0 0 1 1
43 1 10.0 9 19 19 19 19 19
;
title2 'Using Programming Statements';
proc phreg data=tumor;
model time*dead(0) = dose npap /rl;
array pp p1-p15;
array ll l1-l15;
array uu l2-l16;
l1= 27;
l2= 33;
l3= 34;
l4= 37;
l5= 41;
l6= 43;
l7= 45;
l8= 46;
l9= 47;
l10= 49;
l11= 50;
l12= 51;
l13= 53;
l14= 65;
l15= 67;
l16= 71;
if time < l1 then npapillo= 0;
else if time >= l16 then npap= p16;
else
do over uu;
if ll <= time < uu then npap= pp;
end;
run;
data tumor1(keep=id time dead dose t1 t2 npap status);
retain l1 27 l2 33 l3 34 l4 37 l5 41 l6 43 l7 45 l8 46 l9 47
l10 49 l11 50 l12 51 l13 53 l14 65 l15 67 l16 71;
array pp p1-p15;
array qq p2-p16;
array ll l1-l16;
set tumor;
t1 = 0;
t2 = 0;
status = 0;
if ( time = l1 ) then do;
t2= l1;
npap= p1;
status= dead;
output;
end;
else do over pp;
if ( ll = time ) then do;
t2= time;
npap= pp;
status= dead;
output;
end;
else if (ll < time ) then do;
if (pp ^= qq) then do;
if qq = . then t2= time;
else t2= ll;
npap= pp;
status= 0;
output;
t1= t2;
end;
end;
end;
if ( time >= l16 ) then do;
t2= time;
npap= p16;
status= dead;
output;
end;
run;
title2 'Using Counting Process Style of Input';
proc phreg data=tumor1;
model (t1,t2) * status(0) = dose npap;
output out= out3 dfbeta= db_dose db_npap /order=data;
id id time;
run;
proc means data=out3 noprint;
by id time;
var db_dose db_npap;
output out=out4 sum=db_dose db_npap;
run;
/*
title2;
symbol1 v=dot;
goptions ftext=swiss ftitle=swiss htitle=3.5pct htext=3.0pct;
axis1 label = (angle=-90 rotate=90 'DFBETA for DOSE')
minor = none
order = (-.04 to .04 by .01);
axis2 label = (angle=-90 rotate=90 'DFBETA for NPAP')
minor = none
order = (-.030 to .020 by .005);
proc gplot data=out4;
plot db_dose * id /frame hminor=0 vaxis=axis1;
plot db_npap * id /frame hminor=0 vaxis=axis2;
run;
*/
|
