/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~* * MACRO 'EMICM' ------------------------------------------------------------------------------------------------------------------------------------------ The %EMICM macro computes the nonparametric maximum likelihood estimates of the survival curve (i.e., cumulative probabilites) for interval-censoring data. You can choose between the EM algorithm of Turnbull (1976), the ICM algorithm of Groeneboom & Wellner (1992) and the EM-ICM algorithm of Wellner & Zhan (1997). The estimated variance is computed based on the generalized Greenwood formula (Sun 2001). A plot of the survival curves can be displayed using ODS Graphics. The following arguments (separated by commas) can be specified in the %EMICM macro. All arguments are optional except LEFT= and RIGHT=. DATA= SAS data set to be analyzed. Default is DATA=_last_. LEFT= A numeric variable representing the left endpoint of the time interval. Left-censored observations have a missing value or a value of 0. RIGHT= A numeric variable representing the right endpoint of the time interval. Right-censored observations have a missing value. FREQ= A single numeric variable whose values represent the frequency of occurrence of the observations. GROUP= A variable identifying different treatment groups. A separate NPMLE is computed for each value of the GROUP= variable. METHOD= Optimization techique for computing the NPMLE: EM -- Turnbull's self-consistency algorithm ICM -- Iterative Convex Minorant algorithm EMICM -- EM-ICM algorithm Default is METHOD=EMICM. OUT= The output data set containing the NPMLE. OUTITER= The output data set containg the history of iterations. Variables ERROR1, ERROR2, ERROR3, and ERROR4 represent the convergence measures of ERRORTYPE=1, ERRORTYPE=2, ERRORTYPE=3, and ERRORTYPE=4, respectively. ERRORTYPE= Convergence criterion to be used: 1 -- The maximum of the closeness of consecutive estimates, 2 -- The closeness of the log likelihood function, 3 -- The gradient of the log likelihood function, 4 -- The maximum measures of ERRORTYPE=1, ERRORTYPE=2, and ERRORTYPE=3. Default is ERRORTYPE=1. RATECONV= Rate of convergence for the selected ERRORTYPE. Default is RATECONV=1e-7. mRS= Number of resampling for the generalized Greenwood formula. Default is mRS=50. SEED= Random seed for resampling for the generalized Greenwood formula. Default is SEED=8375. TITLE= Primary title for the NPMLE plot. TITLE2= Secondary title for the NPMLE plot. TIMELABEL= Label of time axis in the NPMLE plot. OPTIONS= List of display options (separated by blanks): NOTABLE -- Supressing the table of the NPMLE. PLOT -- Graphical display of the estimated survival curve. *~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/ %macro EMICM (version, data=_last_, left= , right= , group= , freq=, method=EMICM, out= , outiter= , errortype=1, rateconv=1e-7, mRS= 50, seed= 8375, title= , title2= , Timelabel= , options= ) ; %let _version=2.2; %if &version ne %then %put EMICM macro Version &_version; %let _opts = %sysfunc(getoption(notes)); options nonotes; /* Check for newer version */ %if %sysevalf(&sysver >= 8.2) %then %do; %let _notfound=0; filename _ver url 'http://ftp.sas.com/techsup/download/stat/versions.dat' termstr=crlf; data _null_; infile _ver end=_eof; input name:\$15. ver; if upcase(name)="&sysmacroname" then do; call symput("_newver",ver); stop; end; if _eof then call symput("_notfound",1); run; %if &syserr ne 0 or &_notfound=1 %then %put &sysmacroname: Unable to check for newer version; %else %if %sysevalf(&_newver > &_version) %then %do; %put &sysmacroname: A newer version of the &sysmacroname macro is available.; %put %str( ) You can get the newer version at this location:; %put %str( ) http://support.sas.com/ctx/samples/index.jsp; %end; %end; options &_opts; %if %INDEX(%str(&data), .) = 0 %then %let dataname= WORK.&data ; %else %let dataname = &data ; %if %INDEX(%str(&out), .) = 0 %then %let outdn = WORK.&out ; %else %let outdn = &out; %if %INDEX(%str(&outiter), .) = 0 %then %let outiterdn = WORK.&outiter ; %else %let outiterdn = &outiter; ************** EXISTENCE for the required variables in LEFT=, RIGHT= ***************; %if %sysfunc(EXIST(&data))=0 %then %do; %put ERROR: The data set %UPCASE(&dataname) does not exist. ; %goto EXIT; %end; %if (%bquote(&left)=) %then %do; %put ERROR: Expecting a name in 'LEFT= '.; %goto EXIT; %end; %if (%bquote(&right)=) %then %do; %put ERROR: Expecting a name in 'RIGHT= '. ; %goto EXIT; %end; %local savenote; %let savenote=%sysfunc(getoption(notes)); ************** Check if LEFT=, RIGHT=, FREQ=, and GROUP= are in &DATA ***************; Options nonotes; data _datacheck; set &data; run; Options &savenote; %let dsid = %sysfunc(open(%str(_datacheck))); %let xid = %sysfunc(varnum(&dsid, &left)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&left) does not exist.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&left) is non-numeric.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %let xid = %sysfunc(varnum(&dsid, &right)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&right) does not exist.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&right) is non-numeric.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %if (%bquote(&group)^=) %then %do; %let xid = %sysfunc(varnum(&dsid, &group)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&group) does not exist.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %end; %if (%bquote(&freq)^=) %then %do; %let xid = %sysfunc(varnum(&dsid, &freq)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&freq) does not exist.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&freq) is non-numeric.; %let rc1 = %sysfunc(close(&dsid)); %goto EXIT2; %end; %end; %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let _nvars = %sysfunc( attrn(&dsid, NVARS) ); %put NOTE: There were &_nobs observations read from the data set %UPCASE(&dataname).; %let rc1 =%sysfunc( close(&dsid) ); ************** Check if mRS > 1 ***************; %if &mRS<=1 %then %do; %put ERROR: The value for 'mRS=' should be greater than one. ; %goto EXIT; %end; Options nonotes; /* ************** process FREQ= variable ***************; %if (%bquote(&freq)=) %then %let checked=&data; %if (%bquote(&freq)^=) %then %do; data checked; set &data; do i=1 to &freq; output; end; run; %let checked= checked; %end; */ /* Delete obs with invalid time values (valid if L=0, R>=0) or invalid FREQ values */ data checked; drop _rm _zF _zT _bad; retain _rm _zF _zT 0; set &data END=eof; _bad=0; %if %bquote(&freq) ^= %then %do; if &freq < 1 then do; _bad= 1; _zF= 1; end; %end; if (missing(&left)) then &left=0; if (&left < 0) then do; _bad=1; _zT= 2; end; else if missing(&right)=0 and &right < &left then do; _bad=1; _zT= 2; end; _rm= _rm + _bad; if eof then do; call symput('_rmN', trim(left(put(_rm,12.)))) ; call symput('_missType', put(_zF+_zT,1.)); end; if _bad=1 then delete; %if %bquote(&freq) ^= %then %do; do i=1 to &freq; output; end; %end; %else %do; output; %end; run; %let checked=checked; %let dsid = %sysfunc(open(%str(checked))); %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let rc1 =%sysfunc( close(&dsid) ); %if &_nobs=0 %then %do; Options &savenote; %put ERROR: There were no valid observations in the data set %UPCASE(&data).; %goto EXIT; Options nonotes; %end; %else %if &_missType=3 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the LEFT=, RIGHT=, or FREQ= variables.; Options nonotes; %end; %else %if &_missType=2 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the LEFT= or Right= variables.; Options nonotes; %end; %else %if &_missType=1 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the FREQ= variable.; Options nonotes; %end; %if (%bquote(&group)=) %then %do; %let nGrp = 0; %let checked = &checked; %Param_Est (indata= &checked); %if (%bquote(&out)=) %then %let outfile= _data1_; %if (%bquote(&out)^=) %then %do; %let outfile= &out; data &outfile; set _data1_; run; %end; %if (%bquote(&outiter)=) %then %let outiterfile= _data2_; %if (%bquote(&outiter)^=) %then %do; %let outiterfile= &outiter; data &outiterfile; set _data2_; run; %end; %end; %else %do; proc sort data=&checked out=indata; by &group; run; data indata _GrpName_(keep= _Grp_ &group N) ; set indata END=eof; by &group; retain _Grp_ N 0; if first.&group=1 then do; N = 1; _Grp_+1; end; else do; N + 1; end; output indata; if last.&group=1 then output _GrpName_; if eof then call symput('nGrp', trim(left(put(_Grp_,12.)))) ; run; proc print data=_GrpName_ noobs ; var &group N ; title 'Number of Observations by Group'; run; goption reset=title; %do gg=1 %to &nGrp; data checked; set indata; by &group; if _Grp_ = ≫ run; %Param_Est (indata= checked); data _data1_; set _data1_; _Grp_ = ≫ run; data _data2_; set _data2_; _Grp_ = ≫ run; proc append base=_dataoutput_ data=_data1_ force; run; proc append base=_dataiter_ data=_data2_ force; run; %end; %if (%bquote(&out)=) %then %let outfile= _dataoutput_; %else %let outfile= &out; data &outfile; merge _GrpName_ _dataoutput_; by _Grp_; drop _Grp_ N; run; %if (%bquote(&outiter)=) %then %let outiterfile= _dataiter_; %else %let outiterfile= &outiter; data &outiterfile; merge _GrpName_ _dataiter_; by _Grp_; drop _Grp_ N; run; %end; ************** process options ************; %let print= 1; %let plot= 0; %let nn=1; %let token=%scan(&options,&nn,%str( )); %do %while(%bquote(&token)^= ); %let token=%upcase(&token); %if %qsubstr(&token,1,4)=PLOT %then %let plot=1; %else %if %qsubstr(&token,1,4)=NOTA %then %let print=0; %else %do; %let _xrc_= Unrecognized option &token; %put ERROR: &_xrc_..; %end; %let nn=%eval(&nn+1); %let token=%scan(&options,&nn,%str( )); %end; %if &print=1 %then %do; proc print data=&outfile; format Probability Cumulative_Probability Survival_Probability Standard_Error_Survival 8.4; run; %end; %if &plot=1 %then %do; %let L=Lower; %let R=Upper; %let jump=Probability; %if (%bquote(&TimeLabel)=) %then %let TimeLabel= Time; * The _Phantom data set is for plotting only; data _Phantom_(drop=lasts lastt); label t=%bquote(&TimeLabel) s="Survival Probability"; retain lastt lasts; set &outfile END=eof; %if (%bquote(&group)^=) %then %do; by &group; %let _start= first.&group; %let _end= last.&group; %end; %if (%bquote(&group)=) %then %do; %let _start=_n_; %let _end= eof; %end; if &_start=1 then do; /* at time 0 */ t=0; s=1; lb=1; ub=1; output; lastt= 0; lasts= 1; /* first L */ t=&L; s=lasts; lb=s - &Jump; ub=s; output; end; else do; if (&L=lastt) then do; t=&L; s= lasts; lb= s - &Jump; ub= s; output; end; else do; /* previous R */ t=lastt; s= lasts; lb= s; ub= s; output; /* current L */ t=&L; s=lasts; lb=s - &Jump; ub=s; output; end; end; if &_end=1 then do; if (&R <= .) then do; t=&L; s=lasts; lb=s; ub=s; output; end; else do; t=&R; s=lasts - &Jump; lb=s; ub=s; output; end; end; lastt=&R; lasts=lasts - &Jump; run; %if (%bquote(&title)^=) %then %do; %let _char1 = %bquote(%substr(&title,1,1)); %if (%bquote(&_char1) = %str(%') or %bquote(&_char1) = %str(%")) %then %let _titleprint= &title; %else %let _titleprint= %str(%"&title%"); %end; %else %let _titleprint=; %if (%bquote(&title2)^=) %then %do; %let _char2 = %bquote(%substr(&title2,1,1)); %if (%bquote(&_char2) = %str(%') or %bquote(&_char2) = %str(%")) %then %let _titleprint2= &title2; %else %let _titleprint2= %str(%"&title2%"); %end; %else %let _titleprint2=; proc template; define statgraph _ICEPlot; begingraph; %if (%bquote(&title)^=) %then %do; entrytitle %unquote(&_titleprint); %end; %if (%bquote(&title2)^=) %then %do; entrytitle %unquote(&_titleprint2)/ textattrs=GraphValueText; %end; layout overlay; %if (%bquote(&group)^=) %then %do; bandplot x=t limitupper=ub limitlower=lb / type=step group=&group datatransparency=0.5; seriesplot x=t y=s / group= &group name="Survival"; discretelegend "Survival" / location=outside title="&group"; %end; %if (%bquote(&group)=) %then %do; bandplot x=t limitupper=ub limitlower=lb / type=step; seriesplot x=t y=s; discretelegend "Survival" / location=outside ; %end; endlayout; endgraph; end; run; proc sgrender data=_Phantom_ template=_ICEPlot; run; %end; proc datasets NOLIST NOWARN NODETAILS; delete checked indata _data1_ _data2_ _phantom_ _dataoutput_ _dataiter_ _grpname_ ; quit; run; Options &savenote; %if (%bquote(&group)^= ) %then %put NOTE: There are &nGrp groups by the variable %UPCASE(&group).; %if (%bquote(&out)^= ) %then %do; %let dsid = %sysfunc( open(&out, i) ); %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let _nvars = %sysfunc( attrn(&dsid, NVARS) ); %put NOTE: The data set %UPCASE(&outdn) has &_nobs observations and &_nvars variables. ; %let rc2 =%sysfunc( close(&dsid) ); %end; %if (%bquote(&outiter)^= ) %then %do; %let dsid = %sysfunc( open(&outiter, i) ); %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let _nvars = %sysfunc( attrn(&dsid, NVARS) ); %put NOTE: The data set %UPCASE(&outiterdn) has &_nobs observations and &_nvars variables. ; %let rc3 =%sysfunc( close(&dsid) ); %end; %put NOTE: EMICM macro completed.; %EXIT2: Options nonotes; proc datasets NOLIST NOWARN NODETAILS; delete _datacheck; quit; run; Options &savenote; %EXIT: %mend EMICM; %macro Param_Est (indata= ); proc iml; use &indata; read all var{&Left &Right}; L= &Left; R=&Right; * handle exact event times ; ExactT= sum(L=R); if ExactT>0 then do; ExactTime = loc(L=R); L[ExactTime, ] = L[ExactTime, ] - 0.00001; end; * assign max(R)*2 for . right-censoring ; MissR = sum(R <= . ); if MissR> 0 then do; locR = loc(R <= .); maxR= max(R); maxL= max(L); if maxR > maxL then maxR2= maxR*2; else maxR2= maxL*2; R[locR, ] = maxR2 ; end; nobs=nrow(L); /********************************* construct the non-overlapping intervals (q, p] and determine the number of pseudo-parameters m **********************************/ 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]); m= ncol(q); p= J(1, m, .); do i=1 to m; do j= npp to 1 by -1; if (pp[j] > q[i]) then p[i] = pp[j]; end; end; sj = unique(p); m=ncol(sj); alpha= J(nobs, m, 0); do i=1 to nobs; alpha[i, ]= (sj > L[i,1])#(sj <= R[i,1]); end; /******************************************************/ /* Module to calculate isotonic estimator of x /******************************************************/ start isotonic(x,weight); n_X = nrow(x); phat=x; wt=weight; plevel= 1:n_X ; dummy=max(phat)*1000; is_vlt = ((phat//dummy) <(-dummy//phat)) # t(1:(n_X +1)); vlt_loc =min( (^is_vlt )*10*n_X + is_vlt); do while (vlt_loc < 10*n_X); idx= t(loc(phat[1:vlt_loc, ]=phat[vlt_loc-1]) || vlt_loc) ; phat[idx, ] = sum(wt[idx] #phat[idx]) / sum(wt[idx]) ; is_vlt = ((phat//dummy) < (-dummy//phat)) # t(1:(n_X +1)); vlt_loc =min( (^is_vlt )*10*n_X + is_vlt); end; return(phat); finish isotonic; if m>1 then do; error&errortype=1; epsilon= 0.2; * STEP 1: set up the initials for beta; iter=1; beta_old=J(m-1,1,1/m)#t(1:m-1); phat_old= (beta_old//1) - (0//beta_old); loglik= sum(log(alpha*phat_old)) ; alpha_dff = (alpha[ ,1:m-1]-alpha[ ,2:m]) ; %if (%bquote(%upcase(&method))^= EM) and (%bquote(%upcase(&method))^= ICM) and (%bquote(%upcase(&method))^= EMICM) %then %put WARNING: METHOD=&method is unrecognized option. METHOD=EMICM is used.; /******************************************************/ /* Module for EM-Iterative Convex Minorant Algorithm /******************************************************/ do while (error&errortype > &rateconv); beta_up = beta_old ; %if (%bquote(%upcase(&method))= EM) %then %goto SKIP_ICM; * STEP 2: Apply ICM ; ap= alpha*phat_old ; if sum(ap =0) > 0 then ap [loc(ap =0)] =0.00001; d1l_n= alpha_dff / repeat(ap, 1, m-1); d1l= t(d1l_n[+, ]); d2l_n= - (d1l_n##2); d2l=t(d2l_n[+, ]); W=diag(-d2l); x_new= beta_old + d1l/vecdiag(W); if sum(x_new<0) then x_new[loc(x_new<0)]=0.00001; beta_up = isotonic(x_new, vecdiag(W)) ; phat_up = (beta_up//1) - (0//beta_up); ap_up = alpha*phat_up; if sum(ap_up<=0) > 0 then ap_up[loc(ap_up<=0)] =0.00001; loglik_up = sum(log(ap_up)) ; /*** line search ***/ if loglik_up >= loglik then beta_new= beta_up; if loglik_up < loglik then do; lambda=1; gamma= 0.5; zz= beta_up; phat_old= (beta_old//1) - (0//beta_old) ; ap_old = alpha*phat_old ; if sum(ap_old=0) > 0 then ap_old[loc(ap_old=0)] =0.00001; d1l_n_old= alpha_dff / repeat(ap_old , 1, m-1); d1l_old= t(d1l_n_old[+, ]); loglik_old= sum(log(alpha*(phat_old))) ; loglik_zz= loglik_up ; lb = loglik_old + epsilon*t(d1l_old) *(zz-beta_old); ub= loglik_old + (1-epsilon)*t(d1l_old) *(zz-beta_old) ; do while ( (loglik_zz < lb) | (loglik_zz > ub) ); if loglik_zz < lb then lambda= lambda-gamma; if loglik_zz > ub then lambda= lambda+gamma; zz= beta_old + lambda*(beta_up - beta_old); gamma = gamma/2 ; ap_zz = alpha*( (zz//1) - (0//zz) ); if sum(ap_zz=0) > 0 then ap_zz[loc(ap_zz=0)] =0.00001; loglik_zz= sum( log(ap_zz) ) ; lbtemp = loglik_old + epsilon*t(d1l_old) *(zz-beta_old); ubtemp= loglik_old + (1-epsilon)*t(d1l_old) *(zz-beta_old) ; lb= min(lbtemp, ubtemp); ub=max(lbtemp, ubtemp); end; beta_new = zz; end; beta_up = beta_new; %SKIP_ICM: %if (%bquote(%upcase(&method))= ICM) %then %goto SKIP_EM; * STEP 3: Apply Self-consistency ; phat_old= (beta_up //1) - (0//beta_up ); ap = alpha*phat_old; if sum(ap =0) > 0 then ap[loc(ap =0)] =0.00001; p_temp= (alpha#t(phat_old))/repeat(ap , 1, m); p_new = t(p_temp[+, ]/nobs); beta_new=cusum(p_new)[1:m-1]; %SKIP_EM: * STEP 4: Convergence checking ; p_new = (beta_new //1) - (0//beta_new); ap_new = alpha*p_new ; if sum(ap_new =0) > 0 then ap_new[loc(ap_new =0)] =0.00001; loglik_new= sum( log(ap_new) ) ; d1l_n= alpha_dff / repeat(ap_new, 1, m-1); d1l= t(d1l_n[+, ]); error2= max(abs(loglik_new - loglik)); error1= max(abs(beta_new - beta_old)); error3 = sum(beta_new#d1l); error4 = max(error1, error2, abs(error3)); summary_itertp = iter || loglik_new || error1 || error2 || error3 ||error4 ; if iter=1 then summary_iter = summary_itertp ; else summary_iter = summary_iter // summary_itertp; if error&errortype > &rateconv then do; beta_old = beta_new; iter= iter +1; phat_old = p_new; loglik=loglik_new; end; end; beta_new = beta_new // 1; surv = 1-beta_new ; /******************************************************/ /* Generalized Greenwood formula based on resampling /******************************************************/ mRS= &mRS; * number of ReSampling; S_RS = J(mRS, 1, .); S_KM= J(mRS, m, .); do b=1 to mRS; Tib= J(nobs, 1, .); deltaib= J(nobs, 1, .); do ii=1 to nobs; if R[ii,1] > sj[ ,m] then do; deltaib[ii, ]= 0; Tib[ii, ]=L[ii , ]; end; else if R[ii,1] <= sj[ ,m] then do; deltaib[ii, ]= 1; Li= sum(sj <= L[ii ,1]) ; if Li=0 then SLi=1; else SLi= 1 - beta_new[Li, ]; Ri= sum(sj <= R[ii ,1]) ; SRi= 1 - beta_new[Ri, ]; fi = (sj>L[ii ,1])#(sj<= R[ii ,1])#t(p_new) / (SLi - SRi); positivef = loc(fi>1e-12); fi= fi[ ,positivef] / sum(fi[ ,positivef]); * to make sure sum(fi)=1::: required to run "RANDMULTINOMIAL" ; call randseed(&seed); Tib[ii, ] = sj[ ,sum(RANDMULTINOMIAL(1, 1, fi) # positivef)]; end; end; njb = repeat(Tib, 1, m) >= repeat(sj, nobs, 1); njb= njb[+, ]; djb= (repeat(deltaib#Tib, 1, m) = repeat(sj, nobs, 1)); djb = djb[+, ]; S_KMb = J(1, m, 0); nonzero = loc(njb > djb); if sum(njb > djb)>0 then S_KMb[1, nonzero] = (njb[1, nonzero]-djb[1, nonzero])/njb[1, nonzero]; S_KMb_temp = S_KMb[1,1]; do jj=2 to m; S_KMb_temp = S_KMb_temp*S_KMb[1,jj]; S_KMb[1,jj] = S_KMb_temp; end; S_KM[b, ] = S_KMb; end; *end of do b; S_KM_bar= S_KM[: ,]; phat= beta_new - (0//beta_new[1:m-1]); ap = alpha*phat; if sum(ap<=0) >0 then do; djp=J(nobs, m, 0); nonzeroi = loc(ap>0); djp[nonzeroi, ] = alpha[nonzeroi, ] #repeat(t(phat), nobs, 1)[nonzeroi, ] / repeat(ap[nonzeroi, ] , 1, m); end; else do; djp = alpha#repeat(t(phat), nobs, 1) / repeat(ap, 1, m); end; djp = djp[+, ]; njp = J(1, m, .); do jj=1 to m; njp[1, jj] = sum(djp[1, jj:m]); end; nonzero = loc(njp > djp); V_G0 = J(1, m, 0); if sum(njp > djp)>0 then V_G0[1, nonzero] = (djp[1, nonzero] / (njp[1, nonzero] #(njp[1, nonzero] - djp[1, nonzero] )) ); V_G = ((1-beta_new)##2) # t(cusum(V_G0)) ; *note: if nj=dj, then S(sj)=0 ; V_GG = V_G + t((S_KM - repeat(S_KM_bar, mRS, 1))[##, ])/(mRS-1); end; /******************************************************/ /* output /******************************************************/ tq = t(q); if ExactT > 0 then do; ExactUq = unique( L[ExactTime, ]); do u=1 to ncol(ExactUq); ET = ExactUq[1, u] ; if sum(tq = ExactUq[1, u])>0 then tq[loc(tq = ExactUq[1, u] )] = ExactUq[1, u] + 0.00001; end; end; tp= t(p) ; if (MissR > 0) then do; if sum(tp=maxR2) > 0 then tp[loc(tp=maxR2)] = . ; end; if m=1 then do; p_new = 1; beta_new=1; surv=0; V_GG=0; summary_iter = {0 . . . . .}; end; SE_GG = sqrt(V_GG); output = tq || tp|| p_new || beta_new || surv || SE_GG; coln= {'Lower' 'Upper' 'Probability' 'Cumulative Probability' 'Survival Probability' 'Standard Error Survival'}; colniter={'Iteration' 'LogLik' 'Error1' 'Error2' 'Error3' 'Error4'}; create _data1_ from output [colname=coln]; append from output; close _data1_; create _data2_ from summary_iter [colname=colniter]; append from summary_iter; close _data2_; quit; %mend Param_Est;