## 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;

/********************************************************
*******************************************************/
g= j(1,nparm,0);
tmp= _x # (1 /  (_x * theta`) );
g= tmp[+,];
return(g);

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

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

```