ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssc
c     This program is written by Xiaomei Liao.                         c
c     Last modified: Oct,16,2009                                       c
c                                                                      c
c     It applies on HPFS data the new method developed in the paper of c
c     "Survival analysis with error-prone time-varying covariates: a   c
c     risk set calibration approach" by Xiaomei Liao, David Zucker,    c
c     Yi Li, Donna Spiegelman.                                         c
c                                                                      c
c     The user can change the value of 'nex' and 'ncov' to define the  c
c     dimension of exposures and covariates. nex>0 and ncov>=0.        c 		
ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssc

      Module variable

      real*8 pi
      parameter (pi=3.14159265358979323846)

      integer nmain, nval, nage_val, nage_main, nageobs, ncov
      parameter (nmain=47759, nval=573, nageobs=36, nagedist=56)
      parameter (nage_val=10, nage_main=10, nex=2, ncov=0)
      integer j_val(nval), j_main(nmain)
      integer mainobs, valobs
      integer outfile1
      parameter (mainobs=390694,valobs=573)
c      parameter (mainobs=303692,valobs=573)	

      integer sub_id_val(nval)
      real*8 dtcalc_val(nval,nage_val,nex)
      real*8 ffqcalc_val(nval,nage_val,nex)
      real*8 age_val(nval,nage_val), bc_val(nval,nage_val)
      real*8 t_val(nval), del_val(nval)
      real*8 zt_depval(nval,nage_val,nex)
      real*8 zs_depval(nval,nage_val,nex+ncov)
cc    add in covariates
      real*8 covar_val(nval,nage_val,0:ncov)

      integer sub_id_main(nmain)
      real*8 age_main(nmain,nage_main), bc_main(nmain,nage_main)
      real*8 t_main(nmain), del_main(nmain)
      real*8 zs_depmain(nmain,nage_main,nex)
cc    add in covariates
      real*8 covar_main(nmain,nage_main,0:ncov)

      real*8 zs_common(nmain,nage_main,nex)
      real*8 covar_common(nmain,nage_main,0:ncov)
      
      integer ndist,nfa,nal,ilow(nmain), iup(nmain),fnum(nmain)
      integer fsub(nageobs,2), nfa2, ngrp
      real*8 evcnt(nmain), tf(nageobs), ageobs(nmain)

      integer ry_val(nageobs), ry_main(nageobs), case_cnt
      integer cntcom(nageobs), count1(nageobs), count2(nageobs)
      integer count3(nageobs), count4(60)

      real*8 zshat_depmain(nmain,nage_main,nex)
      real*8 covarhat_main(nmain,nage_main,ncov)
      real*8 xxcom(500,60,nex+ncov),yycom(500,60,nex)
      real*8 acom(60,1+nex+ncov,nex), r22(60)
      integer id_val(500,60)

      real*8 betaes_naiv(nex+ncov), serror_naiv(nex+ncov)
      real*8 betaes_rrc(nex+ncov), serror_rrc(nex+ncov)

      real*8 age_interval, age_low(60), age_up(60)

      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss


      program hpfs_data

      use variable
      implicit none

      integer i,j,k,k1,b, jtemp

      integer sub_id_valtmp(valobs)
      real*8 dtcalc_valtmp(valobs,nex), ffqcalc_valtmp(valobs,nex)
      real*8 age_valtmp(valobs), bc_valtmp(valobs)

      integer sub_id_maintmp(mainobs)
      real*8 ffqcalc_maintmp(mainobs,nex), ind_main(nmain,nage_main)

cc    add in covariates
      real*8 covar_maintmp(mainobs,0:ncov)
      real*8 covar_valtmp(valobs,0:ncov)

      real*8 age_maintmp(mainobs), bc_maintmp(mainobs)
      real*8 a, age_temp(43), age_start, age_end
      real*8 xtmp(nageobs), xhattmp(nageobs)
      real*8 xave(nageobs), bxave(nageobs)
      real*8 c_corr(nageobs), x_corr(nageobs)
      integer xcnt(nageobs), xhatcnt(nageobs), cnt2(nageobs)
      integer cnmiss(nageobs), xnmiss(nageobs), casecnt_age(nageobs)
      real*8 a0temp(nageobs), a1temp(nageobs)


c     outfile1 is the unit number for the output file
      outfile1=1011

      count2=-1
      count3=-1
      count4=-1


c     the new validation data set with both HPFS and EATS

      open(101, file='hpfsVal_cov_con_fin.dat', action='read')
      do i=1,valobs
         if (ncov.ge.1) then
            read(101,*) sub_id_valtmp(i),age_valtmp(i),
     *           (dtcalc_valtmp(i,k),ffqcalc_valtmp(i,k), k=1,nex),
     *           (covar_valtmp(i,j),j=1,ncov)
         else
            read(101,*) sub_id_valtmp(i),age_valtmp(i),
     *           (dtcalc_valtmp(i,k),ffqcalc_valtmp(i,k),k=1,nex)
         endif
      enddo

      close(101)


c     check the reading
c      do i=1,valobs
c         if (sub_id_valtmp(i).eq.165307) then
c            print *, sub_id_valtmp(i), age_valtmp(i), 
c     *           (dtcalc_valtmp(i,k),k=1,nex),
c     *           (ffqcalc_valtmp(i,k),k=1,nex)        
c         endif
c      enddo


c     way to check the correctness of reading
c      print *,'i=',i
c      if (ncov.ge.1) then
c         print *, sub_id_valtmp(1),age_valtmp(1),dtcalc_valtmp(1,1:nex),
c     *        ffqcalc_valtmp(1,1:nex),covar_valtmp(1,0:ncov)
c      endif


c     Main study data set for HPFS

      open(201,file='hpfsMain_cov_con_fin.dat', action='read')                                                              
      do i=1,mainobs
         if (ncov.ge.1) then
            read(201,*) sub_id_maintmp(i),age_maintmp(i),bc_maintmp(i),
     *           ffqcalc_maintmp(i,1:nex),covar_maintmp(i,1:ncov)  
         else
            read(201,*) sub_id_maintmp(i),age_maintmp(i),bc_maintmp(i),
     *           ffqcalc_maintmp(i,1:nex)
         endif  
      enddo      

      close(201)

c      way to check the correctness of reading
c      print *,'i=',i
c      print *,sub_id_maintmp(1),age_maintmp(1),bc_maintmp(1),
c     *     ffqcalc_maintmp(1,1:nex),covar_maintmp(1,0:ncov)



c     transformation of main study
      k=1
      j=1
      sub_id_main(k)=sub_id_maintmp(k)
      zs_depmain(k,j,1:nex)=ffqcalc_maintmp(k,1:nex)
      age_main(k,j)=age_maintmp(k)
      bc_main(k,j)=bc_maintmp(k)
      if (ncov.ge.1) then
         covar_main(k,j,1:ncov)=covar_maintmp(k,1:ncov)
      endif
      j_main(k)=j

      do i=2,mainobs
         if(sub_id_maintmp(i).eq.sub_id_maintmp(i-1)) then
            j=j+1
            sub_id_main(k)=sub_id_maintmp(i)
            zs_depmain(k,j,1:nex)=ffqcalc_maintmp(i,1:nex)
            age_main(k,j)=age_maintmp(i)
            bc_main(k,j)=bc_maintmp(i)
            if (ncov.ge.1) then
               covar_main(k,j,1:ncov)=covar_maintmp(i,1:ncov)
            endif
            j_main(k)=j
         else         
            k=k+1
            j=1
            j_main(k)=j
            sub_id_main(k)=sub_id_maintmp(i)
            zs_depmain(k,j,1:nex)=ffqcalc_maintmp(i,1:nex)
            age_main(k,j)=age_maintmp(i)
            bc_main(k,j)=bc_maintmp(i)
            if (ncov.ge.1) then
               covar_main(k,j,1:ncov)=covar_maintmp(i,1:ncov)
            endif
         endif
      enddo
      
c      print *, '# of subjects in the main study: nmain=', k
c      print *, 'The last subject in the main study:', sub_id_main(k)
c      print *, 'The last subject data:', zs_depmain(k,1:j,1:nex), 
c     *age_main(k,1:j)

c     To determine follow-up time and censoring indicator 
      do i=1,nmain
         t_main(i)=age_main(i,j_main(i))
         del_main(i)=bc_main(i,j_main(i))
      enddo


c     tesing for the main study to see there's still missing 
c     value in it
      do i=1,nmain
         do j=1,j_main(i)
            if (zs_depmain(i,j,1).eq.-1.0d0) then
               print *,'Main',sub_id_main(i), zs_depmain(i,j,1:nex),
     *              age_main(i,j),bc_main(i,j)
            endif
         enddo
      enddo


c     writing main study data out
      open(401, file='newMain1.dat', action='write')
      do k=1,nmain
         do j=1,j_main(k)
            if (ncov.ge.1) then
               write(401,15) sub_id_main(k),
     *              age_main(k,j),bc_main(k,j),
     *              (zs_depmain(k,j,i),i=1,nex),
     *              (covar_main(k,j,i),i=1,ncov)
            else
               write(401,16) sub_id_main(k),
     *              age_main(k,j),bc_main(k,j),
     *              (zs_depmain(k,j,i),i=1,nex)
            endif
         enddo
      enddo

      print *, '# of prostate cancer cases in the main study:', 
     *     sum(del_main)


 15   format(1X,I6,2X,F6.2,2X,F4.1,2X,1(F8.4,2X),1(F8.4,2X))
 16   format(1X,I6,2X,F6.2,2X,F4.1,2X,1(F8.4,2X))
       

      call sort_main(nmain,nage_main,sub_id_main,zs_depmain,
     *     del_main,t_main, age_main,j_main, bc_main, 
     *     covar_main)
      

      print *, 'The following results are for the naive cox method: '
      call naivecoxmethod(zs_depmain,covar_main)

c      stop

c     tranformation of validation study
      k=1
      j=1
      sub_id_val(k)=sub_id_valtmp(k)
      dtcalc_val(k,j,1:nex)=dtcalc_valtmp(k,1:nex)
      ffqcalc_val(k,j,1:nex)=ffqcalc_valtmp(k,1:nex)
      age_val(k,j)=age_valtmp(k)
      if (ncov.ge.1) then
         covar_val(k,j,1:ncov)=covar_valtmp(k,1:ncov)
      endif
      j_val(k)=j

      do i=2,valobs
         if (sub_id_valtmp(i).eq.sub_id_valtmp(i-1)) then
            j=j+1
            sub_id_val(k)=sub_id_valtmp(i)
            dtcalc_val(k,j,1:nex)=dtcalc_valtmp(i,1:nex)
            ffqcalc_val(k,j,1:nex)=ffqcalc_valtmp(i,1:nex)
            age_val(k,j)=age_valtmp(i)
            if (ncov.ge.1) then
               covar_val(k,j,1:ncov)=covar_valtmp(i,1:ncov)
            endif
            j_val(k)=j
         else
            k=k+1
            j=1
            j_val(k)=j
            sub_id_val(k)=sub_id_valtmp(i)
            dtcalc_val(k,j,1:nex)=dtcalc_valtmp(i,1:nex)
            ffqcalc_val(k,j,1:nex)=ffqcalc_valtmp(i,1:nex)
            age_val(k,j)=age_valtmp(i)
            if (ncov.ge.1) then
               covar_val(k,j,1:ncov)=covar_valtmp(i,1:ncov)
            endif
         endif
      enddo

            
c     To determine follow-up time and censoring indicator 
      do i=1,nval
         t_val(i)=age_val(i,j_val(i))
      enddo


c     writing the validation data out      
c     original data set
      open(501, file='newVal1.dat', action='write')
      do k=1,nval
         do j=1,j_val(k)
            if (ncov.ge.1) then
               write(501,17) sub_id_val(k),age_val(k,j),
     *              (dtcalc_val(k,j,i),i=1,nex), 
     *              (ffqcalc_val(k,j,i),i=1,nex),
     *              (covar_val(k,j,i),i=1,ncov)
            else
               write(501,18) sub_id_val(k),age_val(k,j),
     *              (dtcalc_val(k,j,i),i=1,nex), 
     *              (ffqcalc_val(k,j,i),i=1,nex)
            endif
         enddo
      enddo

      
 17   format(1X,I6,2X,F6.2,2X,1(F8.4,2X),1(F8.4,2X),1(F8.4,2X))
 18   format(1X,I6,2X,F6.2,2X,1(F8.4,2X),1(F8.4,2X))
 

      call sort_val(nval,nage_val,sub_id_val,dtcalc_val,
     *     ffqcalc_val,t_val,age_val,j_val,covar_val)


      zt_depval(:,:,1:nex)=dtcalc_val
      zs_depval(:,:,1:nex)=ffqcalc_val
      if (ncov.ge.1) then
         zs_depval(:,:,(1+nex):(nex+ncov))=covar_val(:,:,1:ncov)
      endif

      print *, 'The following results are for the RRC method: '
      call rrc
      call rrccoxmethod(zshat_depmain, covar_main)


c     summarizing the result

      call data_summary

c     corresponding to the program 
      end


ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss


      subroutine sort_main(nmax,ntm,subid,zs,del,t,t_age,j_ind,censor,
     *     cov_vec)

      use variable
      implicit none

c     on the entry variable
      integer nmax, ntm, j_ind(nmax), subid(nmax)
      real*8 zs(nmax,ntm,nex), del(nmax), censor(nmax,ntm)
      real*8 t(nmax), t_age(nmax,ntm)
      real*8 cov_vec(nmax,ntm,0:ncov)

c     dummy variable
      integer wksp(nmax)
      real*8 wksp1(nmax),wksp2(nmax,ntm)
      real*8 wksp3(nmax,ntm,0:ncov)
      real*8 wksp4(nmax,ntm,nex)
      real*8 told,rsknum,rm, rsk(nmax)
      integer indx(nmax), ix,i,j
 
      external dsvrgp
  
c     SORT DATA
      do i = 1, nmax
         indx(i) = i
      enddo
      call dsvrgp(nmax,t,t,indx)

      do i = 1, nmax
         wksp(i) = j_ind(i)
      enddo
      do i = 1, nmax
         j_ind(i) = wksp(indx(i))
      enddo
      do i = 1, nmax 
         wksp(i) = subid(i) 
      enddo 
      do i = 1, nmax 
         subid(i) = wksp(indx(i)) 
      enddo 

      do i = 1, nmax
         wksp1(i) = del(i)
      enddo
      do i = 1, nmax
         del(i) = wksp1(indx(i))
      enddo
      do i = 1, nmax
         wksp4(i,:,:) = zs(i,:,:)
      enddo
      do i = 1, nmax
         zs(i,:,:) = wksp4(indx(i),:,:)
      enddo
      do i = 1, nmax
         wksp2(i,:) = t_age(i,:)
      enddo
      do i = 1, nmax
         t_age(i,:) = wksp2(indx(i),:)
      enddo
      do i = 1, nmax
         wksp2(i,:) = censor(i,:)
      enddo
      do i = 1, nmax
         censor(i,:) = wksp2(indx(i),:)
      enddo
      if (ncov.ge.1) then
         do i = 1, nmax
            wksp3(i,:,1:ncov)= cov_vec(i,:,1:ncov)
         enddo
         do i = 1, nmax
            cov_vec(i,:,1:ncov) = wksp3(indx(i),:,1:ncov)
         enddo
      endif


c     IDENTIFY UNIQUE FOLLOW-UP TIME
      ix = 1
      told = t(1)
      ilow(1) = 1
      rsknum = nmax
      rsk(1) = rsknum
      rm = 1.0d0
      evcnt(1) = del(1)
      ageobs(1)=told
      do i  = 2, nmax
         if (t(i) .ne. told) then
            iup(ix) = i-1
            rsk(ix) = rsknum
            rsknum = rsknum - rm
            ix = ix+1
            ilow(ix) = i
            rm = 0.0d0
            evcnt(ix) = 0.0d0
            told = t(i)
            ageobs(ix)=told
         endif
         rm = rm + 1.0d0
         evcnt(ix) = evcnt(ix) + del(i)
         if (i .eq. nmax)  iup(ix) = nmax
      enddo
      ndist = ix
c      print *, 'ndist=', ndist

c     IDENTIFY UNIQUE FAILURE TIME
      nfa=0
      tf=0.0d0
      fnum=0
      do i=1,ndist
         if (evcnt(i).ne.0.0d0) then
            nfa=nfa+1
            fsub(nfa,1)=ilow(i)
            fsub(nfa,2)=iup(i)
            j=ilow(i)
            tf(nfa)=t(j)
c            fnum(nfa)=nint(t(j)/dt1)
         endif
      enddo
c      nal=sum(fnum)

c      print *, '# of unique failure time=', nfa

      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss

      subroutine sort_val(nmax,ntm,subid,zt,zs,t,t_age,j_ind,cov_vec)
c     this subroutine is used to sort the data in the validation study

      use variable

      implicit none
  
c     variables on the entry
      integer nmax, ntm, j_ind(nmax),subid(nmax)
      real*8 cov_vec(nmax,ntm,0:ncov)
      real*8 zt(nmax,ntm,nex),zs(nmax,ntm,nex),t(nmax)
      real*8 t_age(nmax,ntm)

c     dummy variables
      integer i, indx(nmax),wksp(nmax)
      real*8  wksp1(nmax),wksp2(nmax,ntm)
      real*8 wksp3(nmax,ntm,0:ncov), wksp4(nmax,ntm,nex) 
  
c     imsl subroutine and functions used here
c     dsvrgp
  

c     SORT DATA
      do i = 1,nmax
         indx(i) = i
      enddo
      call dsvrgp(nmax,t,t,indx)

      do i = 1, nmax
         wksp(i) = j_ind(i)
      enddo
      do i = 1, nmax
         j_ind(i) = wksp(indx(i))
      enddo
      do i = 1, nmax 
         wksp(i) = subid(i) 
      enddo 
      do i = 1, nmax 
         subid(i) = wksp(indx(i)) 
      enddo 

      do i = 1, nmax
         wksp4(i,:,:) = zt(i,:,:)
      enddo
      do i = 1, nmax
         zt(i,:,:) = wksp4(indx(i),:,:)
      enddo
      do i = 1, nmax
         wksp4(i,:,:) = zs(i,:,:)
      enddo
      do i = 1, nmax
         zs(i,:,:) = wksp4(indx(i),:,:)
      enddo
      
      do i = 1, nmax
         wksp2(i,:) = t_age(i,:)
      enddo
      do i = 1, nmax
         t_age(i,:) = wksp2(indx(i),:)
      enddo
    
      if (ncov.ge.1) then
         do i = 1, nmax 
            wksp3(i,:,1:ncov)= cov_vec(i,:,1:ncov) 
         enddo 
         do i = 1, nmax 
            cov_vec(i,:,1:ncov) = wksp3(indx(i),:,1:ncov) 
         enddo 
      endif
    

      return
      end

ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss


      subroutine naivecoxmethod(z,covar)

      use variable
      implicit none

      integer nd, ifail
      parameter(nd=nex+ncov)
      real*8 betaes(nd), st_error(nd), z(nmain,nage_main,nex)
      real*8 covar(nmain,nage_main,0:ncov)

      if (nd.ge.2) then
         call cox_bcoah(nd,z,covar,betaes,st_error,ifail)
      else if (nd.eq.1) then
         call cox_zbren(nd,z,betaes,st_error,ifail)
      else
         print *, 'The dimension of data is zero!'
      endif
         
      betaes_naiv=betaes
      serror_naiv=st_error

      if (ifail.eq.0) then
         print *,'naive','beta=', betaes(1), 'se=', st_error(1)
      else
         print *,'The method does not converge.'
      endif
      
      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 
      subroutine cox_bcoah(nd,z,covar,betaes,st_error,ifail)

      use variable
      implicit none

c     variables on the entry
      integer nd, ifail
      real*8 z(nmain,nage_main,nex),covar(nmain,nage_main,0:ncov)
      real*8 betaes(nd), st_error(nd)
      
c     dummy variables
      real*8 fvalue,fscale 
      real*8 xguess(nd),xlb(nd),xub(nd),xscale(nd),x(nd) 
      real*8 hes_fval(nd,nd),hesinv(nd,nd),eval(nd)
 
      integer iparam(7), rparam(7),ibtype
      integer isev, iwk(nd),i, ier 
      integer n1rty,iercd

      external rrccox_likelihood, rrcgradient_cox
      external rrchessian_cox, dbcoah


      xguess= 0.2d0 
c      print *, 'initial guess=', xguess
  
      ibtype=0
      xlb=-2.0d0 
      xub=2.0d0

      xscale=1.0d0
      fscale=1.0d0
      iparam(1)=0

      zs_common=z
      covar_common=covar
      call dbcoah(rrccox_likelihood,rrcgradient_cox, rrchessian_cox
     *     ,nd,xguess,ibtype,xlb,xub,xscale,fscale,iparam,rparam,x,
     *     fvalue)
 
      ier = iercd()
      isev = n1rty(1)
      ifail = 1
c      print *, 'isev=', isev

      if (isev .eq. 0) then
         
         betaes=x
c         print *,'bcoah estimates=', betaes

         call rrchessian_cox(nd,betaes,hes_fval,nd)
c         print *, 'hessian_ana=', hes_fval

         call devlsf(nd,hes_fval,nd,eval)
c         print *, 'eigenvalue=',eval
         if (minval(eval).gt.0.0d0) then
            ifail=0
            call dlinds(nd,hes_fval,nd,hesinv,nd)
            Do i=1,nd
               st_error(i)=dsqrt(hesinv(i,i))
            End Do    
         endif

c         print *, 's.e.=', st_error
         
      endif
 
      return
      end 

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss

      subroutine cox_zbren(nd,z,betaes,st_error,ifail)

      use variable
      implicit none

c     variables on the entry
      integer nd, ifail
      real*8 betaes(nd), st_error(nd)
      real*8 z(nmain,nage_main,nex)
      
c     dummy variables
      integer maxfn, isev, n1rty
      real*8 dx,abmx,za,zb,errabs,errrel
      real*8 hess(nd,nd), derv
      real*8 coxscore
      
      external dzbren, coxscore

      zs_common=z

      dx=0.10d0
      abmx=0.1d0
      za = abmx-dx
      zb = abmx+dx
      do while (coxscore(za)*coxscore(zb).gt.0.0)
         za = za-dx
         zb = zb+dx
      end do
      
      errabs = 0.0d0
      errrel = 0.00001d0
      maxfn = 300
      call dzbren(coxscore,errabs,errrel,za,zb,maxfn)

      isev = n1rty(1)
      ifail = 1
      st_error=0.0
      if (isev .eq. 0) then 
         betaes = zb
         call rrchessian_cox(nd,betaes,hess,nd)
         derv=hess(1,1)
         if (derv.gt.0.0) then
            ifail=0
            st_error=dsqrt(1.0d0/derv)
         endif
      endif


      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss

      double precision function coxscore(bval)

      implicit none

c     variables on the entry
      real*8 bval

c     dummy variables
      integer nd
      parameter (nd=1)
      real*8 bval(nd), g_fval(nd),coxscore

      call rrcgradient_cox(nd,bval,g_fval)

      coxscore=g_fval(1)

      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss


      subroutine rrc

      use variable
      implicit none

      integer pp,qq
      parameter(pp=nex,qq=nex+ncov)
      integer i,j,k,k1,i1,j1,j2
      integer ifail
      
      real*8 beta_rrc, st_rrc, st_rrcs, st_rrcs1
      real*8 gamma(qq+1,pp), temp1(pp,1), temp2(qq+1,1)


c     Define the risk process ry_main and ry_val for each observation time
      ry_main=0
      i=1
      do j=1,nfa
         do while (t_main(i).lt.tf(j)) 
            if (i.lt.nmain) then
               i=i+1
            else
               ry_main(j)=nmain+1
               go to 10
            endif              
         enddo
         ry_main(j)=i
 10   enddo
      
      ry_val=0
      i=1
      do j=1,nfa
         do while (t_val(i).lt.tf(j)) 
            if (i.lt.nval) then
               i=i+1
            else 
               ry_val(j)=nval+1
               go to 20
            endif
         enddo
         ry_val(j)=i
 20   enddo



c     Calculate the regression coefficients in the risk sets of 
c     the validation study.
c     Grouping idea is using here. Not exactly in this version.

      ngrp= 23
      age_low(1)=tf(1)
      age_up(1)=tf(1)
 
      do i=2,ngrp-1
         age_low(i)=tf(i)
         age_up(i)=age_low(i)
      enddo

      age_low(ngrp)=tf(ngrp)
      age_up(ngrp)=tf(nfa)

c     if the censoring age of a subject is not younger than age_up()
c     then this subject is in the risk set (age_low, age_up).

      id_val=0
c     special dealing with the first entry
      k=1
      i1=0
      do i=1,nval
         if (t_val(i).ge.age_up(k)) then
            do j1=1,j_val(i)
               if ((age_val(i,j1).ge.age_low(k))
     *              .and.(age_val(i,j1).le.age_up(k))
     *              .and.(zt_depval(i,j1,1).ne.-1.0d0)
     *              .and.(zs_depval(i,j1,1).ne.-1.0d0)
     *              ) then
                  i1=i1+1
                  xxcom(i1,k,1:(nex+ncov))=
     *                 zs_depval(i,j1,1:(nex+ncov))
                  yycom(i1,k,1:nex)=zt_depval(i,j1,1:nex)
                  id_val(i1,k)=sub_id_val(i)
               endif
            enddo
         endif
      enddo
      count4(k)=i1
      

      do k=2, ngrp-1
         i1=0
         do i=1,nval
            if (t_val(i).ge.age_up(k)) then
               do j1=1,j_val(i)
                  if ((age_val(i,j1).ge.age_low(k))
     *                 .and.(age_val(i,j1).le.age_up(k))
     *                 .and.(zt_depval(i,j1,1).ne.-1.0d0)
     *                 .and.(zs_depval(i,j1,1).ne.-1.0d0)
     *                 ) then
                     i1=i1+1
                     xxcom(i1,k,1:(nex+ncov))=
     *                    zs_depval(i,j1,1:(nex+ncov))
                     yycom(i1,k,1:nex)=zt_depval(i,j1,1:nex)
                     id_val(i1,k)=sub_id_val(i)
                  endif
               enddo
            endif
         enddo
         count4(k)=i1
      enddo


c     special dealing with the last entry
      k=ngrp
      i1=0
      do i=1,nval
         if (t_val(i).ge.age_low(k)) then
            do j1=1,j_val(i)
               if ((age_val(i,j1).ge.age_low(k))
     *              .and.(age_val(i,j1).le.age_up(k))
     *              .and.(zt_depval(i,j1,1).ne.-1.0d0)
     *              .and.(zs_depval(i,j1,1).ne.-1.0d0)
     *              ) then
                  i1=i1+1
                  xxcom(i1,k,1:(nex+ncov))=
     *                 zs_depval(i,j1,1:(nex+ncov))
                  yycom(i1,k,1:nex)=zt_depval(i,j1,1:nex)
                  id_val(i1,k)=sub_id_val(i)
               endif
            enddo
         endif
      enddo
      count4(k)=i1


      do j=1,ngrp
         if (count4(j).ge.2) then
            call mvlm2(pp,qq,count4(j),yycom(1:count4(j),j,:),
     *           xxcom(1:count4(j),j,:),gamma)

            acom(j,:,:)=gamma 
c            print *, 'j,acom(j,k1)=', j, acom(j,:,:)
         endif
      enddo
     

c     construct the x_hat value from the regression coefficients
c     acom(60,1+nex+ncov,nex)
      zshat_depmain=-100.0d0
      do j=1,nfa

         do k=1,ngrp
            if ((tf(j).ge.age_low(k)).and.(tf(j).le.age_up(k))) then
               do i=ry_main(j),nmain
                  do j1=1,j_main(i)
                     if (age_main(i,j1).eq.tf(j)) then
                        if (ncov.ge.1) then
                           temp2(:,1)=(/1.0d0,zs_depmain(i,j1,1:nex),
     *                          covar_main(i,j1,1:ncov)/)
                        else
                           temp2(1:(1+nex),1)=
     *                          (/1.0d0,zs_depmain(i,j1,1:nex)/)
                        endif


                        call dmxtyf(qq+1,pp,acom(k,:,:),qq+1,qq+1,1,
     *                       temp2,qq+1,pp,1,temp1,pp)
                        zshat_depmain(i,j1,:)=temp1(:,1)
                     endif
                  enddo
               enddo
            endif
         enddo
 
      enddo


      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
      
      subroutine rrccoxmethod(z, covar)
      
      use variable
      implicit none
         
      integer nd, ifail
      parameter(nd=nex+ncov)
      real*8 z(nmain,nage_main,nex), covar(nmain,nage_main,0:ncov)
      real*8 betaes(nd), st_error(nd)
      real*8 st_errors(nd), st_errors1(nd)

      if (nd.ge.2) then
         call cox_bcoah(nd,z,covar,betaes,st_error,ifail)
      else if (nd.eq.1) then
         call cox_zbren(nd,z,betaes,st_error,ifail)
      else
         print *, 'The dimension of data is zero!'
      endif
      betaes_rrc=betaes

      if (ifail.eq.0) then
         call rrc_sderror(betaes, st_error, st_errors, st_errors1)
         serror_rrc=st_errors
         print *, 'All variance=', serror_rrc
         print *, 'rrc','beta=', betaes_rrc(1), 'se=', serror_rrc(1)
      else     
         print *,'The method does not converge.'
      endif
      
      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss

      subroutine rrc_sderror(bval,st_error, st_errors, st_errors1)
c     this subroutine is used to calculate the variance estimator of
c     rrc estimates. 

      use variable
      implicit none

c     variables on the entry
      integer nd
      parameter (nd=nex+ncov)
      real*8 bval(nd),st_error(nd),st_errors(nd),st_errors1(nd)

c     dummy variables
      integer i,j, npa
      parameter (npa=(1+nex+ncov)*nex)
      real*8 fchat(nd,nd),fbhat(nd,nd),fbinv(nd,nd),fbsand(nd,nd)
      real*8 ftmp(nd,nd), paramatrix(nd,nd), paratemp1(nd,nd)
      real*8 st_erbe1(nd), st_erbe1s(nd)
      real*8 g_score(npa*ngrp,npa*ngrp)
      real*8 u2matrix(npa*ngrp,npa*ngrp)
      real*8 tempmatrix(npa*ngrp,npa*ngrp)
      real*8 st_erpsi(npa*ngrp,npa*ngrp)
      real*8 g_fval(npa*ngrp,nd),tempmatrix2(nd,npa*ngrp)

c     imsl subroutines used here
c     mrrrr,mxytf


      call score2rrc(bval,fchat)
      call rrchessian_cox(nd,bval,fbhat,nd)


c     inverse hessian matrix part 
      call dlinds(nd,fbhat,nd,fbinv,nd) 
      Do i=1,nd
         st_erbe1(i)=dsqrt(fbinv(i,i))
      End Do

c     pure sandwich variance estimator part
      call dmrrrr(nd,nd,fbinv,nd,nd,nd,fchat,nd,nd,nd,ftmp,nd)
      call dmrrrr(nd,nd,ftmp,nd,nd,nd,fbinv,nd,nd,nd,fbsand,nd)
      Do i=1,nd 
         st_erbe1s(i)=dsqrt(fbsand(i,i)) 
      End Do 


c     parameters' variability part
      call parascore_grad(ngrp, st_erpsi)
      call rrcscore_para(ngrp,bval,g_fval)

      
      call dmxtyf(npa*ngrp,nd,g_fval,npa*ngrp,npa*ngrp,npa*ngrp,
     *     st_erpsi,npa*ngrp,nd,npa*ngrp,tempmatrix2,nd)

      call dmrrrr(nd,npa*ngrp,tempmatrix2,nd,npa*ngrp,nd,g_fval,
     *     npa*ngrp,nd,nd,paramatrix,nd)
      

      call dmrrrr(nd,nd,fbinv,nd,nd,nd,paramatrix,nd,nd,nd,
     *paratemp1,nd)
      call dmrrrr(nd,nd,paratemp1,nd,nd,nd,fbinv,nd,nd,nd,
     *paramatrix,nd)


c     three kinds of variance estimators
c     st_errors1: pure sandwich variance estimator
c     st_errors:  sandwich variance estimator + variability of parameters
c              :  the most accurate one.
c     st_error:   non-sandwich variance estimator + variability of 
c             :   parameters. not so mathematically correctly but
acceptable
c             :   in this situation.

      do i=1,nd
         st_error(i)=dsqrt(fbinv(i,i)+paramatrix(i,i))
         st_errors(i)=dsqrt(fbsand(i,i)+paramatrix(i,i))
      enddo

      st_errors1=st_erbe1s

      print *, 'Classical variance=', st_erbe1
      print *, 'Pure sandwich variance=', st_errors1
      print *, 'non_sandwich variance + variability of estimator=', 
     *     st_error
      
      return
      end

ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss


      subroutine score2rrc(bval,fval)
c     This is the updated one with the second term suggested by David.
c     02/21/09, finished the multi-dimensional version.

      use variable  
      implicit none

c     variables on the entry
      integer nd
      parameter (nd=nex+ncov)
      real*8 bval(nd),fval(nd,nd)

c     dummy variables
      real*8 didl
      parameter (didl=1.0d-8)
      Integer j, jtemp, i, k, j1, indk
      real*8 rrf1(nmain,nageobs), rrf2(nmain,nageobs,nd)
      real*8 s0(nageobs), s1(nageobs,nex+ncov), s0a(nageobs), s0b
      real*8 ss2(nmain,nex+ncov), ftmp(nd,1), temp, fvaltmp(nd,nd)
      real*8 temp1(nd,1), temp2(nd,1)

      rrf1=0.0d0; rrf2=0.0d0
      temp1(:,1)=bval
      do j=1,nfa
         do i=ry_main(j),nmain
            do j1=1,j_main(i)
               if ((age_main(i,j1).eq.tf(j)).and.
     *              (zshat_depmain(i,j1,1).ne.-100.0d0)) then
                  if (ncov.ge.1) then
                     temp2(:,1)=(/zshat_depmain(i,j1,1:nex),
     *                    covar_main(i,j1,1:ncov)/)
                  else
                     temp2(1:nex,1)=zshat_depmain(i,j1,1:nex)
                  endif
                  call dmxtyf(nd,1,temp1,nd,nd,1,temp2,nd,1,1,temp,1)

                  rrf1(i,j)=dexp(temp)
                  rrf2(i,j,1:nex)=rrf1(i,j)*zshat_depmain(i,j1,1:nex)
                  if (ncov.ge.1) then
                     rrf2(i,j,(nex+1):nd)=rrf1(i,j)*
     *                    covar_main(i,j1,1:ncov)
                  endif
               endif
            enddo
         enddo
      enddo


      s0=0.0d0
      s1=0.0d0
      do j=1,nfa
         do i=ry_main(j), nmain
            if (rrf1(i,j).ne.0.0d0) then
               s0(j)=s0(j)+rrf1(i,j)
               do k=1,nd
                  s1(j,k)=s1(j,k)+rrf2(i,j,k)
               enddo
            endif
         enddo
      enddo

 
      ss2=0.0d0
      Do i = 1,nmain
         Do j= 1,nfa
            if (rrf1(i,j).ne.0.0d0) then
               s0a(j) = dmax1(s0(j),didl)
               if (t_main(i).ge.tf(j)) then
                  do k=1,nd
                     ss2(i,k)=ss2(i,k)+rrf1(i,j)* 
     *   (rrf2(i,j,k)/rrf1(i,j)-s1
					
