Survival Curve for Interval Censored Data
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IMNLPEX6 */
/* TITLE: Survival Curve for Interval Censored Data */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: IML */
/* DATA: */
/* */
/* SUPPORT: saswmh/sasghg UPDATE: frwick 5/2015 */
/* REF: Technical Report Examples */
/* MISC: */
/* */
/****************************************************************/
data carcin;
input id l r @@;
datalines;
1 20 22 11 30 32 21 22 24 31 34 36
2 22 24 12 32 34 22 22 24 32 34 36
3 26 28 13 32 34 23 28 30 33 36 38
4 26 28 14 32 34 24 28 30 34 38 40
5 26 28 15 34 36 25 32 34 35 38 40
6 26 28 16 36 38 26 32 34 36 42 44
7 28 30 17 42 44 27 32 34 37 42 44
8 28 30 18 30 50 28 32 34 38 46 48
9 30 32 19 34 50 29 32 34 39 28 50
10 30 32 20 20 22 30 32 34 40 48 50
;
proc iml;
use carcin;
read all var{l r};
close carcin;
nobs= nrow(l);
/*********************************************************
construct the nonoverlapping intervals (Q,P) and
determine the number of pseudo-parameters (NPARM)
*********************************************************/
pp= unique(r); npp= ncol(pp);
qq= unique(l); nqq= ncol(qq);
q= j(1,npp, .);
do;
do i= 1 to npp;
do j= 1 to nqq;
if ( qq[j] < pp[i] ) then q[i]= qq[j];
end;
if q[i] = qq[nqq] then goto lab1;
end;
lab1:
end;
if i > npp then nq= npp;
else nq= i;
q= unique(q[1:nq]);
nparm= ncol(q);
p= j(1,nparm, .);
do i= 1 to nparm;
do j= npp to 1 by -1;
if ( pp[j] > q[i] ) then p[i]= pp[j];
end;
end;
/********************************************************
generate the X-matrix for the likelihood
********************************************************/
_x= j(nobs, nparm, 0);
do j= 1 to nparm;
_x[,j]= choose(l <= q[j] & p[j] <= r, 1, 0);
end;
/********************************************************
log-likelihood function (LL)
********************************************************/
start LL(theta) global(_x,nparm);
xlt= log(_x * theta`);
f= xlt[+];
return(f);
finish LL;
/********************************************************
gradient vector (GRAD)
*******************************************************/
start GRAD(theta) global(_x,nparm);
g= j(1,nparm,0);
tmp= _x # (1 / (_x * theta`) );
g= tmp[+,];
return(g);
finish GRAD;
/*************************************************************
estimate the pseudo-parameters using quasi-newton technique
*************************************************************/
/* options */
optn= {1 2};
/* constraints */
con= j(3, nparm + 2, .);
con[1, 1:nparm]= 1.e-6;
con[2:3, 1:nparm]= 1;
con[3,nparm + 1]=0;
con[3,nparm + 2]=1;
/* initial estimates */
x0= j(1, nparm, 1/nparm);
/* call the optimization routine */
call nlpqn(rc,rx,"LL",x0,optn,con,,,,"GRAD");
/*************************************************************
survival function estimate (SDF)
************************************************************/
tmp1= cusum(rx[nparm:1]);
sdf= tmp1[nparm-1:1];
/*************************************************************
covariance matrix of the first nparm-1 pseudo-parameters (SIGMA2)
*************************************************************/
mm= nparm - 1;
_x= _x - _x[,nparm] * (j(1, mm, 1) || {0});
h= j(mm, mm, 0);
ixtheta= 1 / (_x * ((rx[,1:mm]) || {1})`);
if _zfreq then
do i= 1 to nobs;
rowtmp= ixtheta[i] # _x[i,1:mm];
h= h + (_freq[i] # (rowtmp` * rowtmp));
end;
else do i= 1 to nobs;
rowtmp= ixtheta[i] # _x[i,1:mm];
h= h + (rowtmp` * rowtmp);
end;
sigma2= inv(h);
/*************************************************************
standard errors of the estimated survival curve (SIGMA3)
*************************************************************/
sigma3= j(mm, 1, 0);
tmp1= sigma3;
do i= 1 to mm;
tmp1[i]= 1;
sigma3[i]= sqrt(tmp1` * sigma2 * tmp1);
end;
/*************************************************************
95% confidence limits for the survival curve (LCL,UCL)
*************************************************************/
/* confidence limits */
tmp1= probit(.975);
*print tmp1;
tmp1= tmp1 * sigma3;
lcl= choose(sdf > tmp1, sdf - tmp1, 0);
ucl= sdf + tmp1;
ucl= choose( ucl > 1., 1., ucl);
/*************************************************************
print estimates of pseudo-parameters
*************************************************************/
q= q`;
p= p`;
theta= rx`;
ParamEst = q || p || theta;
print ParamEst[L="Parameter Estimates"
colname={Q P THETA} format=BEST9.];
/*************************************************************
print survival curve estimates and confidence limits
*************************************************************/
left= {0} // p;
right= q // p[nparm];
sdf= {1} // sdf // {0};
lcl= {.} // lcl //{.};
ucl= {.} // ucl //{.};
SurvEst = left || right || sdf || lcl || ucl;
print SurvEst[L="Survival Curve Estimates and 95% Confidence Intervals"
colname={left right estimate lower upper}];
quit;