/*************************************************************************/ /* */ /* TITLE: PS_FRAIL.MACRO */ /* AUTHOR: Youyi SHU */ /* ( yshu@post.its.mcw.edu ) */ /* VERSION: 1.0 -- April 20, 1997 */ /* 1.01 -- February 14, 2001 */ /* ( minor changes to adapt to SAS V8 ) */ /* PURPOSE: Estimation of the dependence parameter and Regression */ /* Coefficients in the Positive Stable Frailty Model */ /* */ /*************************************************************************/ /* */ /* Reference: Shu, Y. and Klein, J.P. (1999). A SAS Macro for the */ /* Positive Stable Frailty Model. American Statistical */ /* Association Proceedings of the Statistical Computing */ /* Section, p.47-52. */ /* */ /*************************************************************************/ %macro ps_frail( dsn, /* data set name */ greport=, /* grouping information option */ itprint=, /* iteration history option */ itsumm=, /* summary of the iteration history option */ alpha=0.05, /* Confidence level option */ rl=, /* risk limits option */ cov=, /* estimated variance-covariance matrix option */ outest=, /* output data set */ baseline=, /* output data set */ output=, /* output data set */ convbeta=0.001, /* convergence criterion based on beta */ lambda=0.25, /* blending constant in Marquardt's method */ gridsize=0.05, /* step size in the grid search method */ goldleng=0.001 /* tolerance length for golden section search */ ); options ls=80 ps=59 center; /*************************************************************************/ /***************************** PART (A) *******************************/ /* Data Cleaning and Sorting */ /*************************************************************************/ /* Check if data set name has null value */ %if %length(&dsn)=0 %then %do; %put ==============================================; %put !!! THE DATA SET HAS NOT BEEN SPECIFIED !!!; %put ==============================================; %goto over; %end; /* Macro to determine whether a SAS data set exists */ %macro exist(dsn); %global exist; %if &dsn ne %then %str( data _null_; if 0 then set &dsn; stop; run; ); %if &syserr=0 %then %let exist=yes; %else %let exist=no; %mend exist; %exist(&dsn) %if &exist=no %then %do; %put ==================================================; %put !!! THE DATA SET YOU SPECIFIED DOES NOT EXIST !!!; %put ==================================================; %goto over; %end; /* Macro to determine the number of observations in a SAS data set */ %macro numobs(dsn); %global num; data _null_; if 0 then set &dsn nobs=count; call symput('num',left(put(count,12.))); stop; run; %mend numobs; %numobs(&dsn) %if &num=0 %then %do; %put ======================================================; %put !!! THE DATA SET YOU SPECIFIED HAS 0 OBSERVATIONS !!!; %put ======================================================; %goto over; %end; /* Obtain useful information about the data set you specified */ proc contents data=&dsn out=contout noprint; run; data contout; set contout (keep=name type varnum); proc sort; by varnum; run; /* Check if there are at least 4 variables in the data set you specified */ %numobs(contout) %if &num<4 %then %do; %put ==================================================================; %put %str( )!!! THE DATA SET SHOULD HAVE AT LEAST 4 VARIABLES !!!; %put %str( )(i.e. group, time, censoring indicator and at least 1 covariate); %put ==================================================================; %goto over; %end; /* Check if all the variables are of numeric type */ %let findchar=; data _null_; set contout; if type=2 then do; call symput('findchar','yes'); stop; end; run; %if &findchar ne %then %do; %put ============================================================; %put !!! All THE VARIABLES IN THE DATA SET SHOULD BE NUMERIC !!!; %put %str( )(including the variable for group effect ); %put ============================================================; %goto over; %end; /* Put original variables into macro variables */ data _null_; set contout; if _n_=1 then call symput('group',trim(name)); else if _n_=2 then call symput('t',trim(name)); else if _n_=3 then call symput('cens',trim(name)); else call symput('covl'||left(put(_n_-3,4.)),trim(name)); run; /* Macro to Generate the List of Covariate Names */ %macro namelist(number,prefix=covl); %local i; %do i=1 %to &number; "&&&prefix&i" %end; %mend namelist; %let covlist=%namelist(%eval(&num-3)); /*****************************************************************/ /* IML program to discard the observations containing either */ /* missing values or non-positive on study times or both */ /*****************************************************************/ proc iml; reset noname; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Data Cleaning'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; use &dsn; read all into dsnm; n=nrow(dsnm); *n: total number of observations in the raw data; call symput('nraw',trim(left(char(n,12,0)))); nmissrow=loc((dsnm=.)[,+]=0); missrow=loc((dsnm=.)[,+]^=0); tposrow=loc(dsnm[,2]>0); tnposrow=loc(dsnm[,2]^=. & dsnm[,2]<=0); start modifdsn; goodrow=xsect(nmissrow,tposrow); dsnm=dsnm[goodrow,]; close &dsn; create &dsn from dsnm [colname={"&group" "&t" "&cens" &covlist}]; append from dsnm; nrgood=ncol(goodrow); print ,; print "==>" nrgood "out of &nraw observations used."; finish modifdsn; start progends; print 'All observations will be excluded from the analysis due to having either'; print 'missing values or non-positive on study times or both. '; print ,; print '!!! THE DATA SET HAS NO USABLE OBSERVATIONS !!!'; print ' *** Program Ends *** '; call symput('deleall','yes'); finish progends; %let deleall=; mattrib badrow rowname='Case#'; if ncol(missrow)>0 & ncol(tnposrow)=0 then do; badrow=missrow; nrbad=ncol(badrow); *nrbad: total number of bad observations; if nrbad0 then do; badrow=tnposrow; nrbad=ncol(badrow); if nrbad0 & ncol(tnposrow)>0 then do; badrow=unique(missrow||tnposrow); nrbad=ncol(badrow); if nrbad All &nraw observations will be used in the analysis."; end; quit; %if &deleall ne %then %do; %put ================================================; %put !!! THE DATA SET HAS NO USABLE OBSERVATIONS !!!; %put ================================================; %goto over; %end; /*********************************************************************/ /* IML program to check if all the data are censored, or if */ /* there are any illegal codings of the censoring status variable */ /*********************************************************************/ proc iml; use &dsn; read all into dsnm; cens=dsnm[,3]; n=nrow(cens); *n: total number of observations (after data cleaning); nc=ncol(loc(cens=0)); *nc: total number of censored observations (after data cleaning); %let allcens=; if n=nc then do; call symput('allcens','yes'); print '!!! All OBSERVATIONS (AFTER DATA CLEANING) ARE CENSORED !!!'; print ' NO analysis can be performed. '; print ' *** Program Ends *** '; end; else do; /* If censoring status is neither 1 nor 0, then replace it with 1 */ check_c=loc(cens^=1 & cens^=0); if ncol(check_c)>0 then do; cens=choose((cens^=1 & cens^=0),1,cens); ncheck_c=ncol(check_c); print '************************************************************'; print '!!! WARNING MESSAGE !!!'; print '************************************************************'; print "For the data set after cleaning, if an observation's censoring status is"; print 'neither 1 nor 0, then it will be replaced with 1 automatically. '; print ' ( Note: 1-event 0-censored ) '; mattrib ncheck_c colname='' label=''; print '==> There are' ncheck_c 'such cases happened.'; dsnm[,3]=cens; close &dsn; create &dsn from dsnm [colname={"&group" "&t" "&cens" &covlist}]; append from dsnm; end; end; quit; %if &allcens ne %then %do; %put ============================================================; %put !!! All OBSERVATIONS (AFTER DATA CLEANING) ARE CENSORED !!!; %put %str( ) NO analysis can be performed.; %put ============================================================; %goto over; %end; /*****************************************/ /* IMPORTANT: Sorting the Data */ /*****************************************/ proc sort data=&dsn out=_cox_; by descending &t &cens; run; /*************************************************************************/ /***************************** PART (B) *******************************/ /* Statistical Analysis for the Positive Stable Frailty Model */ /*************************************************************************/ /* SAS/IML Program */ /* Description of parameters: */ /* n total number of observations */ /* gb number of groups */ /* gm vector of the size of each group */ /* gd vector of the number of events within each group */ /* p dimension of covariates */ /* gindex vector used to keep track of everybody's group number */ /* t vector of survival times of length n */ /* cens indicator vector for event or censored */ /* ( 1-event 0-censored ) */ /* z A p by n covariate matrix (each column corresponds */ /* to the covariate vector of an individual on study) */ /* nrd number of distinct event times (ignoring group) */ /* ip index vector for the distinct event times after */ /* sorting the whole data by survival time */ /* d vector of the number of events at each distinct */ /* event time */ /* s A p by nrd matrix (each column corresponds to the */ /* sum of the covariate vectors over the individuals */ /* who fail at each distinct event time) */ /*************************************************************************/ proc iml symsize=150; /* Read in matrix of covariate names */ use contout; read all var{name} into covlist; covlist=covlist[4:nrow(covlist)]; /* Read data into IML matrices */ use _cox_; read all into coxm; n=nrow(coxm); col=ncol(coxm); p=col-3; group=coxm[,1]; t=coxm[,2]; cens=coxm[,3]; z=t(coxm[,4:col]); free coxm; /* Pick up the unique grouping value and obtain # of groups */ gvalue=unique(group)`; gb=nrow(gvalue); /* Compute size and # of events per group and assign group index */ gm=j(gb,1,0); gd=j(gb,1,0); gindex=j(n,1,0); do i=1 to gb; do j=1 to n; if group[j]=gvalue[i] then do; gindex[j]=i; gm[i]=gm[i]+1; if cens[j]=1 then gd[i]=gd[i]+1; end; end; end; /* Data profile printing */ print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Data Profile'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; reset nocenter; %if %index(&dsn,.)=0 %then %let dsname=work.&dsn; %else %let dsname=&dsn; print " Data set: %upcase(&dsname)"; print " Grouping Variable: &group"; print " Dependent Variable: &t"; print " Censoring Variable: &cens"; print " Censoring Value: 0"; print " Ties Handling: BRESLOW"; reset center; label=''; mattrib n colname='Total' label=label nevent colname='Event' label=label ncens colname='Censored' label=label percentc colname='%Censored' label=label format=6.2; nevent=sum(cens); ncens=n-nevent; percentc=100*ncens/n; print ,; print ' Summary of the Number of '; print ' Event and Censored Values'; print n nevent ncens percentc; print ,; call symput('mac_gb',trim(left(char(gb,12,0)))); reset nocenter; print " Total number of groups: &mac_gb"; reset center; %if %upcase(&greport)=Y %then %str( print /; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Summary of Group Size and the Number of Events Per Group'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; gseq=(1:gb)`; mattrib gseq colname='Sequence#' label=label gvalue colname='Group' label=label gm colname='Size' label=label gd colname='Events' label=label; Print gseq gvalue gm gd; ); /* Pick up the Distinct Event Time Indices */ ip0=j(n,1,0); k=0; do j=1 to n-1; if t[j]>t[j+1] & cens[j]=1 then do; k=k+1; ip0[k]=j; end; end; if cens[n]=1 then do; k=k+1; ip0[k]=n; end; nrd=k; ip=ip0[1:nrd]; free ip0; ipip={0}//ip; /* Count the # of events at each distinct event time */ d=j(nrd,1,0); s=j(p,nrd,0); do k=1 to nrd; do j=1 to n; if t[j]=t[ip[k]] & cens[j]=1 then do; d[k]=d[k]+1; s[,k]=s[,k]+z[,j]; end; end; end; /******************************************************************/ /* Compute S0 and LL for the Modified Cox Model */ /******************************************************************/ /* LL log partial likelihood */ /******************************************************************/ start s0LL; s0=j(nrd,1,0); do k=1 to nrd; do j=ipip[k]+1 to ipip[k+1]; s0[k]=s0[k]+w[gindex[j]]*exp(t(beta)*z[,j]); end; end; s0=cusum(s0); LL=beta`*s[,+]-sum(d#log(s0)); finish s0LL; /******************************************************************/ /* Compute S1, S2, U, and H for the Modified Cox Model */ /******************************************************************/ /* U score vector of length p */ /* H Hessian matrix (p by p) */ /******************************************************************/ start s12uh; u=j(p,1,0); h=j(p,p,0); s1=j(p,nrd,0); s2=j(nrd*p,p,0); do k=1 to nrd; do j=ipip[k]+1 to ipip[k+1]; y=w[gindex[j]]*exp(t(beta)*z[,j]); s1[,k]=s1[,k]+z[,j]*y; s2[((k-1)*p+1):(k*p),]=s2[((k-1)*p+1):(k*p),]+z[,j]*t(z[,j])*y; end; if k>=2 then do; s1[,k]=s1[,k]+s1[,k-1]; s2[((k-1)*p+1):(k*p),]=s2[((k-1)*p+1):(k*p),] +s2[((k-2)*p+1):((k-1)*p),]; end; u=u+s[,k]-d[k]*s1[,k]/s0[k]; h=h+d[k]*(s1[,k]*t(s1[,k])-s2[((k-1)*p+1):(k*p),]*s0[k])/s0[k]**2; end; g=j(p,p,0); do i=1 to p; g[i,i]=1/sqrt(abs(h[i,i])); end; finish s12uh; /************************************************************/ /* Marquardt's method to estimate regression coefficients */ /* in the modified cox model for a given frailty (W) */ /************************************************************/ start coxmar; maxiter=25; *beta=j(p,1,0); beta=currbeta; run s0LL; run s12uh; oldbeta=beta; oldLL=LL; %if %upcase(&itprint)=Y %then %str( print "initial value: " u LL beta; ); lambda=λ signal=0; do iter=2 to maxiter until(signal=1); /* Beginning of the Beta Loop */ newbeta=oldbeta-g*inv(g*h*g-lambda*i(p))*g*u; beta=newbeta; run s0LL; newLL=LL; if newLL>=oldLL then do; run s12uh; %if %upcase(&itprint)=Y %then %str( print iter lambda u LL beta; ); stoprule=abs(newbeta-oldbeta); if stoprule<=&convbeta then signal=1; else do; oldbeta=newbeta; oldLL=newLL; lambda=lambda/4; end; end; else do; %if %upcase(&itprint)=Y %then %str( print iter lambda u LL ' *** LL decrease ***'; ); lambda=10*lambda; end; end; /* End of the Beta Loop */ finish coxmar; /***************************************************************/ /* compute the cumulative baseline hazard rate and WH matrix */ /***************************************************************/ start cbhr_wh; alpha=d/s0; cbhr=j(n,1,0); do j=1 to n; do k=1 to nrd; if t[ip[k]]<=t[j] then cbhr[j]=cbhr[j]+alpha[k]; end; end; *print cbhr; wh=j(gb,1,0); do i=1 to gb; do j=1 to n; if gindex[j]=i then wh[i]=wh[i]+cbhr[j]*exp(beta`*z[,j]); end; end; *print wh; finish cbhr_wh; /*************************************************************/ /* Compute the expectation of W (frailty) given data */ /*************************************************************/ start w_hat; jnum=j(gb,1,0); jden=j(gb,1,0); do i=1 to gb; if gd[i]=0 then do; jnum[i]=1; jden[i]=1; end; else do; q=gd[i]+1; do m=1 to q; jnum[i]=jnum[i]+omega[q,m]*wh[i]**(-(m-1)*theta); end; q=gd[i]; do m=1 to q; jden[i]=jden[i]+omega[q,m]*wh[i]**(-(m-1)*theta); end; end; end; w=theta*(wh##(theta-1))#(jnum/jden); *print jnum jden w; finish w_hat; /**************************************************************/ /* Full Likelihood Evaluation based on the EM Algorithm */ /**************************************************************/ start likehood; %if %upcase(&itprint)=Y %then %str( print ,; print '######################' theta '######################'; ); /* Compute the Polynomial Omega */ q=max(gd)+1; omega=j(q,q,0); omega[,1]=1; do i=2 to q; omega[i,i]=(theta**(1-i))*gamma(i-theta)/gamma(1-theta); end; do i=3 to q; do j=2 to i-1; omega[i,j]=omega[i-1,j]+omega[i-1,j-1]*((i-1)/theta-(i-(j-1))); end; end; *print omega; maxwit=25; wstop=100; wbeta=j(p,maxwit,0); wbeta[,1]=beta; %if %upcase(&itprint)=Y %then %str( print 'Initial Value (for W_iter):' beta; ); %if %upcase(&itsumm)=Y %then %str( it_hist={1 .}; ); mattrib wit colname='W_iter' label=''; do wit=2 to maxwit while(wstop>&convbeta); /* Beginning of the W Loop */ %if %upcase(&itprint)=Y %then %str( print wit; ); run w_hat; currbeta=beta; run coxmar; run cbhr_wh; wbeta[,wit]=beta; wstop=max(abs(wbeta[,wit]-wbeta[,wit-1])); %if %upcase(&itsumm)=Y %then %str( it_hist=it_hist//(wit||iter); ); end; /* End of the W Loop */ run w_hat; jdh=jden; Lfull=log(theta)*sum(gd)+(theta-1)*sum(gd#log(wh))-sum(wh##theta) +sum(log(jdh))+sum(cens#(z`*beta))+sum(d#log(alpha)); %if %upcase(&itprint)=Y %then %str( print '--->' theta beta Lfull; ); %if %upcase(&itsumm)=Y %then %str( histelem=char(((theta||Lfull)//j(nrow(it_hist)-1,2,.))||it_hist)//dashline; if rungold=1 & starline=0 then do; histelem=j(1,4,' * ')//dashline//histelem; starline=1; end; history=history//histelem; ); finish likehood; /**********************************************/ /* Golden section method of searching MLE */ /**********************************************/ start goldsect; %if %upcase(&itprint)=Y %then %str( PRINT ,; PRINT ,; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; PRINT '+ Start the Golden Section Search ! +'; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; print 'Bracketing triplet: ' a1 a2 a3; ); ratio=0.61803; cratio=1-ratio; x0=a1; x3=a3; if ( abs(a3-a2)>abs(a2-a1) ) then do; x1=a2; x2=a2+cratio*(a3-a2); fu1=a2_like; theta=x2; run likehood; fu2=Lfull; end; else do; x2=a2; x1=a2-cratio*(a2-a1); fu2=a2_like; theta=x1; run likehood; fu1=Lfull; end; do while (abs(x3-x0)>&goldleng); if fu1fu2 then thetahat=x1; else thetahat=x2; theta=thetahat; run likehood; betahat=beta; Lfullhat=Lfull; %if %upcase(&itprint)=Y %then %str( print ,; print '/** final result **/' thetahat betahat Lfullhat; ); finish goldsect; /****************************************/ /* Bi-Section Method */ /****************************************/ start bisect; %if %upcase(&itprint)=Y %then %str( PRINT ,; PRINT ,; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; PRINT '+ Start the Bi-Sectional Search ! +'; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; ); x=1-stepsize; bimax=1-(1E-7)/2; xstop=0; do while (xLfull0 then do; if x>=1-&goldleng/2 then do; thetahat=x; betahat=beta; Lfullhat=xlike; %if %upcase(&itprint)=Y %then %str( print ,; print '/** final result **/' thetahat betahat Lfullhat; ); end; else do; a1=2*x-1; a2=x; a3=1; a2_like=xlike; %if %upcase(&itsumm)=Y %then %str(rungold=1;); run goldsect; end; xstop=1; end; end; if xstop=0 then do; thetahat=x; betahat=beta; Lfullhat=xlike; %if %upcase(&itprint)=Y %then %str( print ,; print '/** final result **/' thetahat betahat Lfullhat; print '*******************************************'; print '==> This is the independence model !'; print '*******************************************'; ); indmodel=1; *indicator of independence model; end; finish bisect; /*****************************************/ /* GRID search method */ /*****************************************/ start gridmod; %if %upcase(&itprint)=Y %then %str( PRINT ,; PRINT ,; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; PRINT '+ Start the Grid Search Method ! +'; PRINT '+++++++++++++++++++++++++++++++++++++++++++++'; ); stepsize=&gridsize; ngrid=ceil(1/stepsize); grid=j(ngrid,1,0); grid[ngrid]=Lfull0; xgrid=1-stepsize; maxig=ngrid-1; gridstop=0; %if %upcase(&itsumm)=Y %then %str( rungold=0; starline=0; ); do ig=maxig to 1 by -1 until(gridstop=1); theta=xgrid; run likehood; grid[ig]=Lfull; if grid[ig]>grid[ig+1] then xgrid=xgrid-stepsize; else do; if ig=maxig then run bisect; else do; a1=xgrid; a2=xgrid+stepsize; a3=xgrid+2*stepsize; a2_like=grid[ig+1]; %if %upcase(&itsumm)=Y %then %str(rungold=1;); run goldsect; end; gridstop=1; end; end; if gridstop=0 then do; a1=0; a2=1-maxig*stepsize; a3=a2+stepsize; a2_like=grid[1]; %if %upcase(&itsumm)=Y %then %str(rungold=1;); run goldsect; end; finish gridmod; /*************************************************************************/ /* For the Frailty Model: Calculation of the Observed Information Matrix */ /*************************************************************************/ START INFORM; %if %upcase(&itprint)=Y %then %str( print ,; print ,; print 'Calculation of the Observed Information Matrix'; print '=========================================================='; ); /* Define the 1st derivative of Gamma function */ start dgamma(x); y=digamma(x)*gamma(x); return (y); finish dgamma; /* Define the 2nd derivative of Gamma function */ start ddgamma(x); y=(trigamma(x)*(gamma(x)**2)+dgamma(x)**2)/gamma(x); return (y); finish ddgamma; /* Define the 1st derivative of Gamma(i-theta)/Gamma(1-theta) */ start deriv(i,theta); y=(gamma(i-theta)*dgamma(1-theta) -dgamma(i-theta)*gamma(1-theta))/(gamma(1-theta)**2); return (y); finish deriv; /* Define the 2nd derivative of Gamma(i-theta)/Gamma(1-theta) */ start dderiv(i,theta); y=( ( ddgamma(i-theta)*gamma(1-theta)-gamma(i-theta) *ddgamma(1-theta) )*(gamma(1-theta)**2) +( gamma(i-theta)*dgamma(1-theta)-dgamma(i-theta)*gamma(1-theta) ) *2*gamma(1-theta)*dgamma(1-theta) )/(gamma(1-theta)**4); return (y); finish dderiv; /* Compute 1st and 2nd derivatives of the polynomial omega */ /* Domega - 1st derivative DDomega - 2nd derivative */ q=max(gd); domega=j(q,q,0); *domega[,1]=0; do i=2 to q; domega[i,i]=(1-i)*(theta**(-i))*gamma(i-theta)/gamma(1-theta) +(theta**(1-i))*deriv(i,theta); end; do i=3 to q; do j=2 to i-1; domega[i,j]=domega[i-1,j]+domega[i-1,j-1]*((i-1)/theta-(i-(j-1))) +omega[i-1,j-1]*(-(i-1))/(theta**2); end; end; *print domega; ddomega=j(q,q,0); *ddomega[,1]=0; do i=2 to q; ddomega[i,i]=(1-i)*(-i)*(theta**(-i-1))*gamma(i-theta)/gamma(1-theta) +2*(1-i)*(theta**(-i))*deriv(i,theta) +(theta**(1-i))*dderiv(i,theta); end; do i=3 to q; do j=2 to i-1; ddomega[i,j]=ddomega[i-1,j]+ddomega[i-1,j-1]*((i-1)/theta-(i-(j-1))) +2*domega[i-1,j-1]*(-(i-1))/(theta**2) +omega[i-1,j-1]*2*(i-1)/(theta**3); end; end; *print ddomega; /* Compute function J[q,s] 's 1st and 2nd derivatives */ jdh_ds=j(gb,1,0); jdh_dss=j(gb,1,0); jdh_dt=j(gb,1,0); jdh_dtt=j(gb,1,0); jdh_dst=j(gb,1,0); do i=1 to gb; if gd[i]>=2 then do; q=gd[i]; swh=wh[i]; do m=1 to q; jdh_ds[i]=jdh_ds[i]+omega[q,m]*(-(m-1)*theta) *swh**(-(m-1)*theta-1); jdh_dss[i]=jdh_dss[i]+omega[q,m]*(-(m-1)*theta) *(-(m-1)*theta-1)*swh**(-(m-1)*theta-2); jdh_dt[i]=jdh_dt[i]+domega[q,m]*swh**(-(m-1)*theta) +omega[q,m]*(-(m-1)*log(swh))*swh**(-(m-1)*theta); jdh_dtt[i]=jdh_dtt[i]+ddomega[q,m]*swh**(-(m-1)*theta) +2*domega[q,m]*(-(m-1)*log(swh))*swh**(-(m-1)*theta) +omega[q,m]*(((m-1)*log(swh))**2)*swh**(-(m-1)*theta); jdh_dst[i]=jdh_dst[i]+domega[q,m]*(-(m-1)*theta) *swh**(-(m-1)*theta-1) +omega[q,m]*(-(m-1))*swh**(-(m-1)*theta-1) *(1+(-(m-1)*theta)*log(swh)); end; end; end; *print jdh jdh_ds jdh_dss jdh_dt jdh_dtt jdh_dst; free omega domega ddomega; /* Calculate the matrices of Ha, Hb, Hba, Hbb */ Ha=j(gb,nrd,0); do i=1 to gb; do k=1 to nrd; do j=1 to n; if gindex[j]=i & t[ip[k]]<=t[j] then Ha[i,k]=Ha[i,k]+exp(beta`*z[,j]); end; end; end; *print Ha; Hb=j(gb,p,0); do i=1 to gb; do m=1 to p; do j=1 to n; if gindex[j]=i then Hb[i,m]=Hb[i,m]+cbhr[j]*z[m,j]*exp(beta`*z[,j]); end; end; end; *print Hb; Hba=j(gb*p,nrd,0); do i=1 to gb; do m=1 to p; do k=1 to nrd; do j=1 to n; if gindex[j]=i & t[ip[k]]<=t[j] then Hba[(i-1)*p+m,k]=Hba[(i-1)*p+m,k]+z[m,j]*exp(beta`*z[,j]); end; end; end; end; *print Hba; Hbb=j(gb*p,p,0); do i=1 to gb; do m=1 to p; do r=1 to p; do j=1 to n; if gindex[j]=i then Hbb[(i-1)*p+m,r]=Hbb[(i-1)*p+m,r]+cbhr[j]*z[m,j]*z[r,j] *exp(beta`*z[,j]); end; end; end; end; *print Hbb; /* Fisher information matrix */ Ftt=0; do i=1 to gb; Ftt=Ftt+gd[i]/theta**2+(wh[i]**theta)*(log(wh[i])**2) +jdh_dt[i]**2/jdh[i]**2-jdh_dtt[i]/jdh[i]; end; Ftb=j(1,p,0); do m=1 to p; do i=1 to gb; Ftb[1,m]=Ftb[1,m]+Hb[i,m]*( wh[i]**(theta-1) +theta*(wh[i]**(theta-1))*log(wh[i])-gd[i]/wh[i] +jdh_dt[i]*jdh_ds[i]/jdh[i]**2-jdh_dst[i]/jdh[i] ); end; end; Fta=j(1,nrd,0); do k=1 to nrd; do i=1 to gb; Fta[1,k]=Fta[1,k]+Ha[i,k]*( wh[i]**(theta-1) +theta*(wh[i]**(theta-1))*log(wh[i])-gd[i]/wh[i] +jdh_dt[i]*jdh_ds[i]/jdh[i]**2-jdh_dst[i]/jdh[i] ); end; end; Fbb=j(p,p,0); do m=1 to p; do r=m to p; do i=1 to gb; Fbb[m,r]=Fbb[m,r]+Hb[i,m]*Hb[i,r]*( gd[i]*(theta-1)/wh[i]**2 +theta*(theta-1)*(wh[i]**(theta-2)) +jdh_ds[i]**2/jdh[i]**2-jdh_dss[i]/jdh[i] ) -Hbb[(i-1)*p+m,r]*( gd[i]*(theta-1)/wh[i] -theta*(wh[i]**(theta-1))+jdh_ds[i]/jdh[i] ); end; end; end; Fba=j(p,nrd,0); do m=1 to p; do k=1 to nrd; do i=1 to gb; Fba[m,k]=Fba[m,k]+Hb[i,m]*Ha[i,k]*( gd[i]*(theta-1)/wh[i]**2 +theta*(theta-1)*(wh[i]**(theta-2)) +jdh_ds[i]**2/jdh[i]**2-jdh_dss[i]/jdh[i] ) -Hba[(i-1)*p+m,k]*( gd[i]*(theta-1)/wh[i] -theta*(wh[i]**(theta-1))+jdh_ds[i]/jdh[i] ); end; end; end; Faa=j(nrd,nrd,0); do u=1 to nrd; do v=u to nrd; do i=1 to gb; Faa[u,v]=Faa[u,v]+Ha[i,u]*Ha[i,v]*( gd[i]*(theta-1)/wh[i]**2 +theta*(theta-1)*(wh[i]**(theta-2)) +jdh_ds[i]**2/jdh[i]**2-jdh_dss[i]/jdh[i] ); end; Faa[u,v]=Faa[u,v]+(u=v)*d[u]/alpha[u]**2; end; end; free jdh jdh_ds jdh_dss jdh_dt jdh_dtt jdh_dst Ha Hb Hba Hbb; fisher=(Ftt||Ftb||Fta)//(j(p,1,0)||Fbb||Fba)//(j(nrd,1+p,0)||Faa); free Ftt Ftb Fta Fbb Fba Faa; f=1+p+nrd; do i=2 to f; do j=1 to i-1; fisher[i,j]=fisher[j,i]; end; end; *print fisher; store _all_; free / fisher; invf=inv(fisher); load _all_; cov=invf[1:(1+p),1:(1+p)]; stderr=sqrt(vecdiag(cov)); p_Wald=j(1+p,1,0); p_Wald[1]=1-probchi((thetahat-1)**2/cov[1,1],1); do i=1 to p; p_Wald[i+1]=1-probchi(betahat[i]**2/cov[i+1,i+1],1); end; %if %upcase(&itprint)=Y %then %str( print cov; print stderr p_Wald; ); FINISH INFORM; /**************************************************************/ /**************** Main Control Program ****************/ /**************************************************************/ /* Compute the Full Likelihood for the Independence Model (theta=1) */ %if %upcase(&itprint)=Y %then %str( print /; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Iteration History'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; iternote={"* THETA: dependence parameter ", "* BETA: regression coefficients vector ", "* U: score vector ", "* LL: log partial likelihood ", "* LFULL0: log full likelihood for Cox model ", "* LFULL: log full likelihood for frailty model ", "* LAMBDA: blending constant in Marquardt's method"}; reset noname; print iternote; reset name; print ,; print 'Regression Coefficients for the Independence Model'; print '===================== ( THETA=1 ) ====================='; ); theta=1; w=j(gb,1,1); currbeta=j(p,1,0); run coxmar; run cbhr_wh; indbeta=beta; indcov=-inv(h); indse=sqrt(vecdiag(indcov)); indp_W=1-probchi(indbeta##2/indse##2,1); indrr=exp(indbeta); indalpha=alpha; indcbhr=cbhr; indxbeta=z`*indbeta; Lfull0=-sum(wh)+sum(cens#indxbeta)+sum(d#log(alpha)); indLpart=Lfull0-sum(d#(log(d)-1)); %if %upcase(&itprint)=Y %then %str( print '--->' theta beta Lfull0; ); %if %upcase(&itsumm)=Y %then %str( dashline=j(1,4,'-------------'); history=dashline //{" THETA " " LFULL " " W_iter " " Marq_ITER "} //dashline//char(theta||Lfull0||{.}||iter)//dashline; ); /* Estimation of the dependence parameter and regression coefficients */ indmodel=0; *indicator of independence model; run gridmod; /* Likelihood Ratio Test of the Dependence Parameter */ mattrib df_LR colname='DF' label=label format=1.0 Chisq_LR colname='Chi-square' label=label format=9.5 p_LR colname='p-value' label=label format=9.5; df_LR=1; if indmodel=0 then Chisq_LR=(-2)*(Lfull0-Lfullhat); else Chisq_LR=0; p_LR=1-probchi(Chisq_LR,df_LR); %if %upcase(&itprint)=Y %then %str( print 'Likelihood Ratio Test (H0: THETA=1) ==> ' df_LR Chisq_LR p_LR; ); /* Standard errors of the estimates */ if indmodel=0 then RUN INFORM; ***show space; /********************************************************/ /* Summary of interation history */ /********************************************************/ %if %upcase(&itsumm)=Y %then %str( print /; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Summary of Iteration History'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; histnote={"* THETA: dependence parameter ", "* LFULL: log full likelihood ", "* W_iter: iteration sequence number for the frailty (W) ", "* Marq_ITER: total number of iterations using Marquardt's", " method for a given frailty (W) "}; reset noname; print histnote; print history; reset name; if rungold=1 then do; reset nocenter; print " Note: The starting triplet for the golden section search is"; reset center noname; print '(' a1 ',' a2 ',' a3 ')'; reset name; end; ); /*******************************************************/ /* Summary of final results */ /*******************************************************/ print /; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print 'Summary of Final Results'; print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'; print '/*******************************************/'; print '/* Usual Cox Independence Model */'; print '/*******************************************/'; reset noname nocenter; print ' Log partial likelihood =' indLpart; print ' Log full likelihood =' Lfull0; reset name center; print ,; print 'Analysis of Maximum Likelihood Estimates'; mattrib variable rowname='' colname='Variable' label=label format=$12. inddf colname='DF' label=label format=1.0 indbeta colname='Estimate' label=label format=9.5 indse colname='Stderr' label=label format=9.5 indp_W colname='p_Wald' label=label format=9.5 indrr colname='RR' label=label format=9.5; variable=covlist; inddf=j(p,1,1); %if %upcase(&rl)^=Y %then %do; print variable inddf indbeta indse indp_W indrr; %end; %else %do; print variable inddf indbeta indse indp_W; zalpha=probit(1-&alpha/2); UCLindrr=indrr#exp(zalpha*indse); LCLindrr=indrr#exp(-zalpha*indse); cc=100*(1-&alpha); mattrib cc format=2.0 LCLindrr colname='Lower' label=label format=9.5 UCLindrr colname='Upper' label=label format=9.5; print ,; print 'Conditional Risk Ratio and'; reset noname; print cc '% Confidence Limits'; reset name; print variable indrr LCLindrr UCLindrr; %end; %if %upcase(&cov)=Y %then %str( print ,; print 'Estimated Covariance Matrix'; mattrib indcov rowname=covlist colname=covlist label=label /*format=9.5*/; reset noname; print indcov; reset name; ); *print ,; *print ,; print /; print '/*******************************************/'; print '/* Positive Stable Frailty Model */'; print '/*******************************************/'; if indmodel=0 then do; tau=1-thetahat; effect={'Depend Parm'}//covlist; df=j(p+1,1,1); estimate=thetahat//betahat; rr_w=exp(betahat); rr_b=rr_w##thetahat; augrr_w={.}//rr_w; augrr_b={.}//rr_b; reset noname nocenter; mattrib thetahat format=7.5 tau format=7.5; print ' Dependence Parameter: THETA =' thetahat; print " Kendall's TAU =" tau; print ' Log full likelihood =' Lfullhat; print ,; reset center spaces=4; print 'Likelihood Ratio Test of Independence Model (H0: THETA=1)'; print df_LR Chisq_LR p_LR; reset name spaces=1; print ,; print 'Analysis of Maximum Likelihood Estimates'; mattrib effect rowname='' colname='Effect' label=label format=$12. df colname='DF' label=label format=1.0 estimate colname='Estimate' label=label format=9.5 stderr colname='Stderr' label=label format=9.5 p_Wald colname='p_Wald' label=label format=9.5 augrr_w colname='RR_within' label=label format=9.5 augrr_b colname='RR_between' label=label format=9.5; %if %upcase(&rl)^=Y %then %do; print effect df estimate stderr p_Wald augrr_w augrr_b; %end; %else %do; print effect df estimate stderr p_Wald; se_beta=stderr[2:(1+p)]; UCLrr_w=rr_w#exp(zalpha*se_beta); LCLrr_w=rr_w#exp(-zalpha*se_beta); se_thb=j(p,1,0); do i=1 to p; se_thb[i]=sqrt( (betahat[i]**2)*cov[1,1] +(thetahat**2)*cov[i+1,i+1] +2*thetahat*betahat[i]*cov[1,i+1] ); end; UCLrr_b=rr_b#exp(zalpha*se_thb); LCLrr_b=rr_b#exp(-zalpha*se_thb); mattrib rr_w colname='RR_within' label=label format=9.5 LCLrr_w colname='Lower' label=label format=9.5 UCLrr_w colname='Upper' label=label format=9.5 rr_b colname='RR_between' label=label format=9.5 LCLrr_b colname='Lower' label=label format=9.5 UCLrr_b colname='Upper' label=label format=9.5; print ,; print 'Conditional Risk Ratio and'; reset noname; print cc '% Confidence Limits'; reset name; print variable rr_w LCLrr_w UCLrr_w rr_b LCLrr_b UCLrr_b; %end; %if %upcase(&cov)=Y %then %str( print ,; print 'Estimated Covariance Matrix'; mattrib cov rowname=effect colname=effect label=label /*format=9.5*/; reset noname; print cov; reset name; ); end; else print '<==> Cox Independence Model.'; /*****************************************/ /* Producing the OUTEST Data Set */ /*****************************************/ %if &outest ne %then %do; if indmodel=1 then do; _est_m=(t(indbeta)//indcov)||j(1+p,1,indLpart)||j(1+p,1,Lfull0); create _est_ from _est_m [colname={&covlist "_LNLIKE_" "_LLFULL_"}]; append from _est_m; _estch_m=j(1+p,1,"BRESLOW")||({"PARMS"}//j(p,1,"COV")) ||({"ESTIMATE"}//covlist); create _estch_ from _estch_m [colname={"TIES" "_TYPE_" "_NAME_"}]; append from _estch_m; end; else do; _est_m=(t(estimate)//cov)||j(2+p,1,Lfullhat); create _est_ from _est_m [colname={"DEP_PARM" &covlist "_LLFULL_"}]; append from _est_m; _estch_m=j(2+p,1,"BRESLOW")||({"PARMS"}//j(1+p,1,"COV")) ||({"ESTIMATE","DEP_PARM"}//covlist); create _estch_ from _estch_m [colname={"TIES" "_TYPE_" "_NAME_"}]; append from _estch_m; end; %end; /*******************************************/ /* Producing the BASELINE Data Set */ /*******************************************/ %if &baseline ne %then %do; if indmodel=1 then /* RECOMMEND: SAS/Proc PHREG */ do; se_alpha=j(nrd,1,.); /* No standard error of ALPHA is available in this case */ _base_m=t[ip]||indalpha||se_alpha||j(nrd,1,1); end; else do; covalpha=invf[(2+p):(1+p+nrd),(2+p):(1+p+nrd)]; se_alpha=sqrt(vecdiag(covalpha)); _base_m=t[ip]||alpha||se_alpha||j(nrd,1,thetahat); end; create &baseline from _base_m [colname={"_time_" "_bhr_" "_sebhr_" "_theta_"}]; append from _base_m; %end; /*******************************************/ /* Producing the OUTPUT Data Set */ /*******************************************/ %if &output ne %then %do; if indmodel=1 then _outp_m=group||t||cens||z`||indxbeta||indcbhr||j(n,1,1); else do; xbeta=z`*betahat; _outp_m=group||t||cens||z`||xbeta||cbhr||j(n,1,thetahat); end; create &output from _outp_m [colname={"&group" "&t" "&cens" &covlist "_xbeta_" "_cbhr_" "_theta_"}]; append from _outp_m; %end; QUIT; %if &outest ne %then %do; data &outest; merge _estch_ _est_; run; %end; %if &baseline ne %then %do; proc sort data=&baseline; by _time_; run; %end; %over: %mend ps_frail;