subroutine erckcmp(iopt,iropt,zval,nref1,nref2) implicit double precision(a-h, o-z) parameter (kmax=70, nbmax=15, nqmax=5, nmax=75000, nevmax=2000) parameter (mxitr=100, tol=1.0d-5,didl=1.0d-5) parameter (hllow=1.0d-8) external dlinrg, dlinds, dmxtyf, dmrrrr, dmurrv, dmxytf, dsrch dimension x(kmax,nbmax), a(kmax,kmax), f(nqmax,kmax,kmax) dimension biggam(nqmax,nqmax), kv(nmax), del(nmax) dimension t0(nmax), t(nmax), mm(nmax) dimension b(kmax,kmax), pih(kmax), cnt(nmax), iwt(nmax) dimension xtil(kmax,nbmax), xtilex(kmax,nbmax) dimension ilow(nmax), iup(nmax) dimension sst(kmax), dcl(nmax,kmax) dimension rsk(kmax), ev(kmax), grn(kmax) dimension sstar(nmax,kmax), green(nmax,kmax), gstar(nmax,kmax) dimension s(nmax,kmax), svar(nmax,kmax), hl(kmax) dimension q(nevmax,kmax,kmax) dimension omega(kmax,kmax), phi(kmax,nqmax) dimension ominv(kmax,kmax), xttomi(nbmax,kmax) dimension desmat(nbmax,nbmax), v(nbmax,nbmax), rhs(nbmax) dimension beta(nbmax), diff(kmax), bdiff(nbmax), evec(kmax) dimension dlhl(kmax), omnew(kmax,kmax), xtinv(nbmax,nbmax) dimension xtinvom(nbmax,nbmax), omni(kmax,kmax) dimension xttomni(nbmax,kmax) dimension evnt(nmax), ikix(kmax) dimension xtfin(kmax,nbmax), omfin(kmax,kmax) dimension tevnt(nevmax), mmev(nevmax) dimension maux(nevmax,nmax), numaux(nevmax) common /data0/ n, k, nb, nq common /data1/ kv, del, t0, t common /data2/ cnt, ilow, iup common /data3/ x, a, f, biggam common /data4/ beta, v, ifail dn = n nb1 = nb + 1 ipth = 1 ifail=0 nref1 = 0 nref2 = 0 C COMPUTE B (INVERSE OF A) call dlinrg(k,a,kmax,b,kmax) C SET UP EXTENDED DESIGN MATRIX do l = 1, k xtil(l,1) = 1.0d0 do ir = 1, nb xtil(l,ir+1) = x(l,ir) enddo enddo C TALLY COUNTS IN EACH CONFIG do l = 1, k pih(l) = 0.0d0 enddo do i = 1, n pih(kv(i)) = pih(kv(i))+1.0d0 enddo C SORT THE DATA AND IDENTIFY UNIQUE F-U TIMES call sort(n,ndist) C INITIALIZATIONS do l = 1, k sst(l) = 1.0d0 grn(l) = 0.0d0 hl(l) = 0.0d0 enddo C LOOP OVER F-U TIMES do m = 1, ndist evnt(m) = 0.0d0 tcurr = t(ilow(m)) C COMPUTATION OF KAPLAN-MEIER AND CUML HAZ ESTIMATES do l = 1, k rsk(l) = 0.0d0 ev(l) = 0.0d0 enddo do i = ilow(m), n if (t0(i) .lt. tcurr) rsk(kv(i))=rsk(kv(i))+1.0d0 enddo do i = ilow(m), iup(m) ev(kv(i)) = ev(kv(i)) + del(i) evnt(m) = evnt(m) + del(i) enddo do l = 1, k p = 0.0d0 if (rsk(l) .gt. 0.0d0) then p = ev(l)/rsk(l) sst(l) = sst(l)*(1.0d0-p) grn(l) = grn(l) + p/(rsk(l)*(1.0d0-p)) endif sstar(m,l) = sst(l) green(m,l) = grn(l) dcl(m,l) = p gstar(m,l) = rsk(l)/pih(l) enddo C END OF COMPUTATION OF K-M AND CUML HAZ ESTIMATES C COMPUTE S = B*SSTAR AND ITS VARIANCE do l1 = 1, k s(m,l1) = 0.0d0 svar(m,l1) = 0.0d0 do l2 = 1, k addtrm = b(l1,l2)*sstar(m,l2) s(m,l1) = s(m,l1) + addtrm svar(m,l1) = svar(m,l1) + (addtrm**2)*green(m,l2) enddo enddo enddo C END OF LOOP OVER F-U TIMES C IDENTIFY FIRST AND LAST TIMES AT WHICH BOUNDS ON S ARE MET istflg=0 do m = 1, ndist iwflg = 1 do l = 1, k chwid = zval*dsqrt(svar(m,l)) shi = s(m,l) + chwid slo = s(m,l) - chwid if ((iropt .eq. 1) .and. (shi .gt. 1.0d0)) then iwflg=0 nref1 = nref1 + 1 endif if (slo .le. 0.0d0) then iwflg=0 nref2 = nref2 + 1 endif enddo iwt(m) = iwflg if ((istflg .eq. 0) .and. (iwflg .eq. 1)) then mlow=m istflg=1 endif if (iwflg .eq. 1) mhi=m enddo if (istflg .eq. 0) then write(10,*) 'Computation Failed - S Bounds Violated' ifail=1 return endif write(6,*) 'FOR EACH LEFT' C IDENTIFY EVENT TIMES nevtim = 0 do m = 1, ndist if (evnt(m) .gt. 0.0d0) then nevtim = nevtim + 1 tevnt(nevtim) = t(ilow(m)) mmev(nevtim) = m endif enddo C FOR EACH LEFT TRUNCATION TIME IDENTIFY THE EVENT TIME (FROM BELOW) do i = 1, n t0i = t0(i) incx = 1 call dsrch(nevtim,t0i,tevnt,incx,index) if (index .ge. 1) mm(i)=mmev(index) if (index .eq. -1) mm(i)=0 if (index .le. -2) mm(i)=mmev(-index-1) enddo write(3,*) 'COMPUTE L' C COMPUTE L'S do m = mlow, mhi if (iwt(m) .eq. 1) then do l = 1, k hlkm = -dlog(s(m,l)) hl(l) = hl(l) + cnt(m)*hlkm enddo endif enddo do i = 1, n mmi = mm(i) if (mmi .ge. 1) then if (iwt(mmi) .eq. 1) then do l = 1, k hlkm = -dlog(s(mmi,l)) hl(l) = hl(l) - hlkm enddo endif endif enddo do l = 1, k hl(l) = hl(l)/dn enddo write(3,*) 'COMPUTE Q' C SET UP AUXILIARY ARRAY FOR Q COMPUTATION iev = 1 m1 = mmev(1) do while ((iev .le. nevtim) .and. (m1 .le. mhi)) mstart = max(m1,mlow) ixm = 0 do i = 1, n mmi = mm(i) if (mmi .ge. mstart) then ixm = ixm + 1 maux(iev,ixm) = mmi endif enddo numaux(iev) = ixm iev = iev + 1 m1 = mmev(min(iev,nevtim)) enddo C COMPUTE Q do ir = 1, k do is = 1, k iev = 1 m1 = mmev(1) do while ((iev .le. nevtim) .and. (m1 .le. mhi)) q(iev,ir,is) = 0.0d0 mstart = max(m1,mlow) do m2 = mstart, mhi if (iwt(m2) .eq. 1) then addtrm = cnt(m2)*sstar(m2,is)/s(m2,ir) q(iev,ir,is) = q(iev,ir,is) + addtrm endif enddo do indaux = 1, numaux(iev) mmi = maux(iev,indaux) if (iwt(mmi) .eq. 1) then addtrm = sstar(mmi,is)/s(mmi,ir) q(iev,ir,is) = q(iev,ir,is) - addtrm endif enddo q(iev,ir,is) = q(iev,ir,is)/dn iev = iev+1 m1 = mmev(min(iev,nevtim)) enddo enddo enddo C COMPUTE OMEGA ASSUMING KNOWN CLASSIFICATION RATES do ir = 1, k do is = 1, k ixm = 0 omega(ir,is) = 0.0d0 do m = 1, mhi if (evnt(m) .gt. 0.0d0) then ixm = ixm + 1 do l = 1, k if (dcl(m,l) .gt. 0.0d0) then addtrm = q(ixm,ir,l)*q(ixm,is,l)/gstar(m,l) addtrm = addtrm*b(ir,l)*b(is,l)*dcl(m,l)/(pih(l)/dn) omega(ir,is) = omega(ir,is) + addtrm endif enddo endif enddo enddo enddo C ADJUST OMEGA FOR ESTIMATION OF CLASSIFICATION RATES if (nq .gt. 0) then do l = 1, k do iu = 1, nq phi(l,iu) = 0.0d0 enddo enddo do m = mlow, mhi if (iwt(m) .eq. 1) then do l = 1, k do iu = 1, nq do ir = 1, k do is = 1, k addtrm = (1.0d0/s(m,l))*b(l,ir) $ *f(iu,ir,is)*s(m,is)*cnt(m) phi(l,iu) = phi(l,iu) + addtrm/dn enddo enddo enddo enddo endif enddo do 90 i = 1, n mmi = mm(i) if (mmi .lt. mlow) goto 90 if (mmi .gt. mhi) goto 90 if (iwt(mmi) .ne. 1) goto 90 do l = 1, k do iu = 1, nq do ir = 1, k do is = 1, k addtrm = (1.0d0/s(mmi,l))*b(l,ir) $ *f(iu,ir,is)*s(mmi,is) phi(l,iu) = phi(l,iu) - addtrm/dn enddo enddo enddo enddo 90 continue do ir = 1, k do is = 1, k do it = 1, nq do iu = 1, nq addtrm = phi(ir,it)*phi(is,iu)*biggam(it,iu) omega(ir,is) = omega(ir,is) + dn*addtrm enddo enddo enddo enddo endif write(3,*) 'COMPUTE ESTIMATES' C COMPUTE WEIGHTED LEAST SQUARES ESTIMATES C SETUP UP LOG-TRANSFORMED RESPONSE VECTOR C AND CORRESPONDING COVARIANCE MATRIX kgud = 0 do l = 1, k if (hl(l) .gt. hllow) then kgud = kgud + 1 ikix(kgud) = l dlhl(kgud) = dlog(hl(l)) endif enddo if (kgud .lt. nb1) then ifail=1 write(3,*) 'Computation Failed - ', $ 'Insufficient Nonzero L Values' return endif do l = 1, kgud do j = 1, nb1 xtfin(l,j) = xtil(ikix(l),j) enddo enddo do ir = 1, kgud ikr = ikix(ir) omfin(ir,ir) = omega(ikr,ikr) omnew(ir,ir) = omega(ikr,ikr)/(hl(ikr)**2) do is = ir+1, kgud iks = ikix(is) omfin(ir,is) = omega(ikr,iks) omfin(is,ir) = omega(ikr,iks) omnew(ir,is) = omega(ikr,iks)/(hl(ikr)*hl(iks)) omnew(is,ir) = omnew(ikr,iks) enddo enddo if (kgud .ne. k) then write(3,*) 'Some strata eliminated, number of strata ', $ 'remaining = ', kgud write(3,*) 'Labels of retained strata are as follows:' write(3,*) (ikix(ir), ir = 1, kgud) endif C SQUARE CASE if (nb1 .eq. kgud) then call dlinrg(nb1,xtfin,kmax,xtinv,nbmax) call dmurrv(nb1,nb1,xtinv,nbmax,nb1,dlhl,ipth,kgud,beta) call dmrrrr(nb1,nb1,xtinv,nbmax,nb1,nb1,omnew,kmax,nb1,nb1, $ xtinvom,nbmax) call dmxytf(nb1,nb1,xtinvom,nbmax,nb1,nb1,xtinv,nbmax, $ nb1,nb1,v,nbmax) goto 100 endif C NON-ITERATIVE ESTIMATE call dlinds(kgud,omnew,kmax,omni,kmax) call dmxtyf(kgud,nb1,xtfin,kmax,kgud,kgud,omni,kmax,nb1,kgud, $ xttomni,nbmax) call dmrrrr(nb1,kgud,xttomni,nbmax,kgud,nb1,xtfin,kmax, $ nb1,nb1,desmat,nbmax) call dlinrg(nb1,desmat,nbmax,v,nbmax) call dmurrv(nb1,kgud,xttomni,nbmax,kgud,dlhl,ipth,nb1,rhs) call dmurrv(nb1,nb1,v,nbmax,nb1,rhs,ipth,nb1,beta) C IF REQUESTED, PROCEED TO ITERATIVE ESTIMATE. OTHERWISE, FINISH. if (iopt .ne. 1) goto 100 call dlinds(kgud,omfin,kmax,ominv,kmax) C WE USE HERE A NAIVE IMPLEMENTATION OF GAUSS-NEWTON FITTING C IF THERE ARE PROBLEMS WITH CONVERGENCE, A MORE SOPHISTICATED C ALGORITHM MAY BE NEEDED itr = 0 icflg = 0 do while ((itr .lt. mxitr) .and. (icflg .eq. 0)) itr = itr + 1 do l = 1, kgud evec(l) = 0.0d0 do ib = 1, nb1 evec(l) = evec(l) + xtfin(l,ib)*beta(ib) enddo evec(l) = dexp(evec(l)) enddo do l = 1, kgud diff(l) = hl(ikix(l)) - evec(l) do ib = 1, nb1 xtilex(l,ib) = evec(l)*xtfin(l,ib) enddo enddo call dmxtyf(kgud,nb1,xtilex,kmax,kgud,kgud,ominv,kmax,nb1, $ kgud, xttomi,nbmax) call dmrrrr(nb1,kgud,xttomi,nbmax,kgud,nb1,xtilex,kmax, $ nb1,nb1,desmat,nbmax) do ib = 1, nb1 desmat(ib,ib) = desmat(ib,ib)+didl enddo call dlinds(nb1,desmat,nbmax,v,nbmax) call dmurrv(nb1,kgud,xttomi,nbmax,kgud,diff,ipth,nb1,rhs) call dmurrv(nb1,nb1,v,nbmax,nb1,rhs,ipth,nb1,bdiff) icflg = 1 do ib = 1, nb1 beta(ib) = beta(ib) + bdiff(ib) rdiff = dabs(bdiff(ib))/dmax1(1.0d0,dabs(beta(ib))) if (rdiff .gt. tol) icflg=0 enddo enddo ifail = 1-icflg 100 continue if (ifail .ne. 0) $ then write(10,*) 'Computation Failed - G-N Did Not Converge' else do ib1 = 1, nb1 v(ib1,ib1) = v(ib1,ib1)/dn do ib2 = ib1+1, nb1 v(ib1,ib2) = v(ib1,ib2)/dn v(ib2,ib1) = v(ib1,ib2) enddo enddo endif return end C********************************************************************** subroutine sort(nobs,ndis) implicit double precision (a-h,o-z) parameter (nmax=75000) external dsvrgp dimension del(nmax), kv(nmax), indx(nmax), wksp(nmax) dimension ilow(nmax), iup(nmax), cnt(nmax) dimension t0(nmax), t(nmax) common /data1/ kv, del, t0, t common /data2/ cnt, ilow, iup C SORT DATA do i = 1, nobs indx(i) = i enddo call dsvrgp(nobs,t,t,indx) do i = 1, nobs wksp(i) = t0(i) enddo do i = 1, nobs t0(i) = wksp(indx(i)) enddo do i = 1, nobs wksp(i) = del(i) enddo do i = 1, nobs del(i) = wksp(indx(i)) enddo do i = 1, nobs wksp(i) = kv(i) enddo do i = 1, nobs kv(i) = wksp(indx(i)) enddo C IDENTIFY UNIQUE F-U TIMES ix = 1 told = t(1) ilow(1) = 1 do i = 2, nobs if (t(i) .ne. told) then iup(ix) = i-1 ix = ix+1 ilow(ix) = i told = t(i) endif if (i .eq. nobs) iup(ix) = nobs enddo ndis = ix do m = 1, ndis ncnt = iup(m) - ilow(m) + 1 cnt(m) = ncnt enddo return end