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(j,k)/s0a(j))/(dble(nmain)*s0a(j)) enddo endif endif End Do End Do fval=0.0d0 do j=1,ndist if (evcnt(j).ne.0.0d0) then do j1=1,nfa if (t_main(ilow(j)).eq.tf(j1)) indk=j1 enddo do i=ilow(j),iup(j) if (rrf1(i,indk).ne.0.0d0) then if (del_main(i).eq.1.0d0) then do k=1,nd ftmp(k,1)= * (rrf2(i,indk,k)/rrf1(i,indk)-ss2(i,k)-s1(indk,k)/s0a(indk)) enddo call dmxytf(nd,1,ftmp,nd,nd,1,ftmp,nd,nd,nd, * fvaltmp,nd) fval=fval+fvaltmp endif endif enddo endif enddo Do i = 1,nmain if (del_main(i).eq.0.0) then do k=1,nd ftmp(k,1)=-ss2(i,k) enddo call dmxytf(nd,1,ftmp,nd,nd,1,ftmp,nd,nd,nd,fvaltmp,nd) fval=fval+fvaltmp endif End Do return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine rrcscore_para(nd,bval,g_fval) c this is the subroutine to calculate the gradient of rrc score w.r.t c the regression parameters (\alpha_0,\alpha_1,\alpha_2). c this is a wrapper for the Adifor generated subroutine 'g_scorerrcpara'. use variable implicit none integer npa parameter (npa=(1+nex+ncov)*nex) integer i, j, k, nd integer g_p_, ldg_aa, ldg_fval real*8 aa(npa),g_aa(npa,npa) real*8 fval(nex+ncov), g_fval(npa*nd,nex+ncov), bval(nex+ncov) real*8 g_ftemp(npa,nex+ncov) g_p_=npa ldg_aa=g_p_ ldg_fval=g_p_ g_aa=0.0d0 do j=1,npa g_aa(j,j)=1.0d0 enddo g_fval=0.0d0 Do j=1,nd if (count4(j).ge.2) then do i=1,1+nex+ncov do k=1,nex aa((i-1)*nex+k)=acom(j,i,k) enddo enddo call g_scorerrcpara(g_p_,j, aa, g_aa, ldg_aa, fval, * g_ftemp, ldg_fval, bval) g_fval(1+npa*(j-1):npa*j,:)=g_ftemp endif Enddo end Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine g_scorerrcpara(g_p_, kk, aa, g_aa, ldg_aa, fval, g_fval *, ldg_fval, bval) C use variable implicit none C real*8 didl parameter (didl = 1.0d-8) integer nd parameter (nd = nex + ncov) integer kk, jtemp, i, indk, k, j, cnttmp, k1 real*8 aa((1 + nd) * nex), a(1 + nd, nex), atran(nex, 1 + nd) real*8 s0, s1(nd), fval(nd), s0a, bval(nd) real*8 temp(1 + nd, 1), temp1(nex, 1), temp2 C integer g_pmax_ parameter (g_pmax_ = (1+nd)*nex) integer g_i_, g_p_, ldg_fval, ldg_aa double precision d2_p, d6_b, d5_b, d4_v, d4_b, d1_p, d2_v, d3_v, * g_fval(ldg_fval, nd), g_a(g_pmax_, 1 + nd, nex) double precision g_aa(ldg_aa, (1 + nd) * nex), g_s0(g_pmax_), g_ *s1(g_pmax_, nd), g_temp2(g_pmax_), g_temp1(g_pmax_, nex, 1), g_s0a *(g_pmax_), g_atran(g_pmax_, nex, 1 + nd) save g_a, g_s0, g_s1, g_temp2, g_temp1, g_s0a, g_atran external g_matrixprod external g_matrixtranspose integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nd do g_i_ = 1, g_p_ g_fval(g_i_, i) = 0.0d0 enddo fval(i) = 0.0d0 C-------- enddo C do i = 1, 1 + nd do j = 1, nex do g_i_ = 1, g_p_ g_a(g_i_, i, j) = g_aa(g_i_, (i - 1) * nex + j) enddo a(i, j) = aa((i - 1) * nex + j) C-------- enddo enddo C call g_matrixtranspose1(g_p_, 1 + nd, nex, a, g_a, g_pmax_, atran *, g_atran, g_pmax_) C do k = 1, nfa do g_i_ = 1, g_p_ g_s0(g_i_) = 0.0d0 enddo s0 = 0.0d0 C-------- do k1 = 1, nd do g_i_ = 1, g_p_ g_s1(g_i_, k1) = 0.0d0 enddo s1(k1) = 0.0d0 C-------- enddo if ((tf(k) .ge. age_low(kk)) .and. (tf(k) .le. age_up(kk))) th *en jtemp = ry_main(k) C 5 if (jtemp .le. nmain) then do j = 1, j_main(jtemp) if (age_main(jtemp, j) .eq. tf(k)) then C temp(1, 1) = 1.0d0 do k1 = 1, nex temp(1 + k1, 1) = zs_depmain(jtemp, j, k1) enddo if (ncov .ge. 1) then do k1 = 1, ncov temp(1 + nex + k1, 1) = covar_main(jtemp, j, k1) enddo endif C call g_matrixprod1(g_p_, nex, 1 + nd, atran, g_atran, g *_pmax_, 1 + nd, 1, temp, nex, 1, temp1, g_temp1, g_pmax_) C do g_i_ = 1, g_p_ g_temp2(g_i_) = 0.0d0 enddo temp2 = 0.0d0 C-------- do k1 = 1, nex do g_i_ = 1, g_p_ g_temp2(g_i_) = bval(k1) * g_temp1(g_i_, k1, 1) + *g_temp2(g_i_) enddo temp2 = temp2 + bval(k1) * temp1(k1, 1) C-------- enddo if (ncov .ge. 1) then do k1 = 1, ncov do g_i_ = 1, g_p_ g_temp2(g_i_) = g_temp2(g_i_) enddo temp2 = temp2 + bval(nex + k1) * covar_main(jtemp, * j, k1) C-------- enddo endif C d3_v = exp(temp2) d1_p = d3_v do g_i_ = 1, g_p_ g_s0(g_i_) = d1_p * g_temp2(g_i_) + g_s0(g_i_) enddo s0 = s0 + d3_v C-------- do k1 = 1, nex d4_v = exp(temp2) d1_p = d4_v d6_b = temp1(k1, 1) * d1_p do g_i_ = 1, g_p_ g_s1(g_i_, k1) = d6_b * g_temp2(g_i_) + d4_v * g_t *emp1(g_i_, k1, 1) + g_s1(g_i_, k1) enddo s1(k1) = s1(k1) + temp1(k1, 1) * d4_v C-------- enddo if (ncov .ge. 1) then do k1 = 1, ncov d3_v = exp(temp2) d1_p = d3_v d5_b = covar_main(jtemp, j, k1) * d1_p do g_i_ = 1, g_p_ g_s1(g_i_, nex + k1) = d5_b * g_temp2(g_i_) + g_ *s1(g_i_, nex + k1) enddo s1(nex + k1) = s1(nex + k1) + covar_main(jtemp, j, * k1) * d3_v C-------- enddo endif C jtemp = jtemp + 1 goto 5 endif enddo endif C cnttmp = 0 do i = fsub(k, 1), fsub(k, 2) if (del_main(i) .eq. 1.0d0) then cnttmp = cnttmp + 1 do j = 1, j_main(i) if (age_main(i, j) .eq. tf(k)) then C temp(1, 1) = 1.0d0 do k1 = 1, nex temp(1 + k1, 1) = zs_depmain(i, j, k1) enddo if (ncov .ge. 1) then do k1 = 1, ncov temp(1 + nex + k1, 1) = covar_main(i, j, k1) enddo endif C call g_matrixprod1(g_p_, nex, 1 + nd, atran, g_atran, * g_pmax_, 1 + nd, 1, temp, nex, 1, temp1, g_temp1, g_pmax_) C do k1 = 1, nex do g_i_ = 1, g_p_ g_fval(g_i_, k1) = g_temp1(g_i_, k1, 1) + g_fval *(g_i_, k1) enddo fval(k1) = fval(k1) + temp1(k1, 1) C-------- enddo if (ncov .ge. 1) then do k1 = 1, ncov do g_i_ = 1, g_p_ g_fval(g_i_, nex + k1) = g_fval(g_i_, nex + k1 *) enddo fval(nex + k1) = fval(nex + k1) + covar_main(i, *j, k1) C-------- enddo endif C endif enddo endif enddo C d2_v = max (s0, didl) if (s0 .gt. didl) then d1_p = 1.0d0 d2_p = 0.0d0 else if (s0 .lt. didl) then d1_p = 0.0d0 d2_p = 1.0d0 else d1_p = 0.5d0 d2_p = 0.5d0 endif do g_i_ = 1, g_p_ g_s0a(g_i_) = d1_p * g_s0(g_i_) enddo s0a = d2_v C-------- do k1 = 1, nd d4_v = s1(k1) / s0a d4_b = -dble(cnttmp) d5_b = d4_b * (1.0d0 / s0a) d6_b = d4_b * (-d4_v / s0a) do g_i_ = 1, g_p_ g_fval(g_i_, k1) = d6_b * g_s0a(g_i_) + d5_b * g_s1(g_i_ *, k1) + g_fval(g_i_, k1) enddo fval(k1) = fval(k1) - dble(cnttmp) * d4_v C-------- enddo C endif C enddo C C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_matrixprod1(g_p_, nar, nac, a, g_a, ldg_a, nbr, nbc, *b, ncr, ncc, c, g_c, ldg_c) C C used to compute the real matrix(or vector) product c=a*b C use variable implicit none C C on the entry variables integer nar, nac, nbr, nbc, ncr, ncc real*8 a(nar, nac), b(nbr, nbc), c(ncr, ncc) C C dummy variables integer i, j, k C integer g_pmax_ parameter (g_pmax_ = (1+nex+ncov)*nex) integer g_i_, g_p_, ldg_c, ldg_a double precision g_c(ldg_c, ncr, ncc), g_a(ldg_a, nar, nac) integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif if ((nac .ne. nbr) .or. (nar .ne. ncr) .or. (nbc .ne. ncc)) then print *, 'Dimension does not match' stop endif C do i = 1, nar do j = 1, ncc do g_i_ = 1, g_p_ g_c(g_i_, i, j) = 0.0d0 enddo c(i, j) = 0.0d0 C-------- do k = 1, nac do g_i_ = 1, g_p_ g_c(g_i_, i, j) = b(k, j) * g_a(g_i_, i, k) + g_c(g_i_, *i, j) enddo c(i, j) = c(i, j) + a(i, k) * b(k, j) C-------- enddo enddo enddo C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_matrixtranspose1(g_p_, nar, nac, a, g_a, ldg_a, at, *g_at, ldg_at) C used to transpose a real matrix(or vector) C use variable implicit none C C on the entry variables integer nar, nac real*8 a(nar, nac), at(nac, nar) C C dummy variables integer i, j C integer g_pmax_ parameter (g_pmax_ = (1+nex+ncov)*nex) integer g_i_, g_p_, ldg_at, ldg_a double precision g_at(ldg_at, nac, nar), g_a(ldg_a, nar, nac) integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nac do j = 1, nar do g_i_ = 1, g_p_ g_at(g_i_, i, j) = g_a(g_i_, j, i) enddo at(i, j) = a(j, i) C-------- enddo enddo C return end C C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine parascore_grad(nd,st_erpsi) c this subroutine is used to calculate the variance of c \psi=(\alpha_0,\alpha_1, \alpha_2) c this way consists of stacking of matrices, which is time consuming. use variable implicit none c variables on the entry integer nd,npa parameter (npa=(1+nex+ncov)*nex) real*8 st_erpsi(npa*nd,npa*nd) c dummy variables integer i,j,k,k1,k2 real*8 aa(npa),ua(npa) integer g_p_, ldg_aa, ldg_ua real*8 g_aa(npa,npa), g_ua(npa,npa), g_uas(npa,npa) real*8 temp1(nex+ncov,1), temp2(nex,1) real*8 g_temp(npa,npa), g_tempj(nd,npa,npa) real*8 uall(nd,nval,npa), u2matrix(nd,nd,npa,npa) real*8 uwork1(1,npa), uwork2(1,npa), uwork(nd,nd,npa,npa) g_p_=npa ldg_aa=g_p_ ldg_ua=g_p_ g_aa=0.0d0 do j=1,npa g_aa(j,j)=1.0d0 enddo uall=0.0d0 Do j = 1,nd g_ua=0.0d0 do i=1,1+nex+ncov do k=1,nex aa((i-1)*nex+k)=acom(j,i,k) enddo enddo Do i= 1,count4(j) temp1(:,1)=xxcom(i,j,:) temp2(:,1)=yycom(i,j,:) call g_parascore(g_p_, temp1, temp2, * aa, g_aa, ldg_aa, ua, g_uas,ldg_ua) g_ua=g_ua+g_uas uall(j,i,:)=ua Enddo g_temp=-g_ua c Compute the inverse of a real symmetric positive definite matrix call dlinds(npa,g_temp,npa,g_tempj(j,:,:),npa) Enddo c This part calculates \sum U*U^t. Do j=1,nd Do i=j,nd u2matrix(j,i,:,:)=0.0d0 Do k1=1,max0(count4(i),count4(j)) Do k2=1,max0(count4(i),count4(j)) if ((id_val(k1,i).ne.0).and.(id_val(k2,j).ne.0).and. * (id_val(k1,i).eq.id_val(k2,j))) then uwork1(1,:)=uall(j,k2,:) uwork2(1,:)=uall(i,k1,:) call dmxtyf(1,npa,uwork1,1,1,npa,uwork2,1,npa,npa, * uwork(j,i,:,:) ,npa) u2matrix(j,i,:,:)=u2matrix(j,i,:,:)+uwork(j,i,:,:) endif Enddo Enddo Enddo Do i=1,j-1 c print *, 'i=', i, 'j=', j c read(*,*) call dtrnrr(npa,npa,u2matrix(i,j,:,:),npa,npa,npa, * u2matrix(j,i,:,:),npa) Enddo Enddo c do i=1,nd c do j =1,nd c print *, 'u2=',i,j,u2matrix(i,j,:,:) c read(*,*) c enddo c enddo c Multiplying the matrices together. Do i=1,nd Do j=i,nd call dmrrrr(npa,npa,g_tempj(i,:,:),npa,npa,npa, *u2matrix(i,j,:,:),npa,npa,npa,uwork(i,j,:,:),npa) call dmrrrr(npa,npa,uwork(i,j,:,:),npa,npa,npa, *g_tempj(j,:,:),npa,npa,npa,u2matrix(i,j,:,:),npa) Enddo Do j=1,i-1 call dtrnrr(npa,npa,u2matrix(j,i,:,:),npa,npa,npa, *u2matrix(i,j,:,:),npa) Enddo Enddo Do i=1,nd Do j=1,nd st_erpsi(1+npa*(i-1):npa*i,1+npa*(j-1):npa*j)= *u2matrix(i,j,:,:) Enddo Enddo return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine g_parascore(g_p_, x, y, aa, g_aa, ldg_aa, ua, g_ua, ldg *_ua) C use variable implicit none integer i, j, nd parameter (nd = nex + ncov) real*8 x(nd, 1), y(nex, 1), a(1 + nd, nex), atran(nex, 1 + nd) real*8 aa((1 + nd) * nex), ua((1 + nd) * nex) real*8 ua0(nex, 1), ua1(nex, nex), ua2(nex, 0:ncov) real*8 uatot(nex, 1 + nd) real*8 temp1(1 + nd, 1), temp2(nex, 1), temp3(nex, 1), temp4(1, *nex) real*8 temp5(0:ncov, 1) C integer g_pmax_ parameter (g_pmax_ = (1+nd)*nex) integer g_i_, g_p_, ldg_aa, ldg_ua double precision g_a(g_pmax_, 1 + nd, nex), g_aa(ldg_aa, (1 + nd *) * nex), g_ua2(g_pmax_, nex, 0:ncov), g_temp4(g_pmax_, 1, nex), g *_uatot(g_pmax_, nex, 1 + nd), g_ua0(g_pmax_, nex, 1), g_ua1(g_pmax *_, nex, nex), g_ua(ldg_ua, (1 + nd) * nex), g_atran(g_pmax_, nex, *1 + nd), g_temp1(g_pmax_, 1 + nd, 1) double precision g_temp2(g_pmax_, nex, 1), g_temp3(g_pmax_, nex, * 1) save g_a, g_ua2, g_temp4, g_uatot, g_ua0, g_ua1, g_atran, g_temp *1, g_temp2, g_temp3 external g_matrixsubtract external g_matrixprod external g_matrixtranspose integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, 1 + nd do j = 1, nex do g_i_ = 1, g_p_ g_a(g_i_, i, j) = g_aa(g_i_, (i - 1) * nex + j) enddo a(i, j) = aa((i - 1) * nex + j) C-------- enddo enddo C temp1(1, 1) = 1.0d0 do i = 1, nd temp1(1 + i, 1) = x(i, 1) enddo C call g_matrixtranspose(g_p_, 1 + nd, nex, a, g_a, g_pmax_, atran *, g_atran, g_pmax_) call g_matrixprod(g_p_, nex, 1 + nd, atran, g_atran, g_pmax_, 1 *+ nd, 1, temp1, g_temp1, g_pmax_, nex, 1, temp2, g_temp2, g_pmax_) call g_matrixsubtract(g_p_, nex, 1, y, temp2, g_temp2, g_pmax_, *ua0, g_ua0, g_pmax_) C do i = 1, nex temp3(i, 1) = x(i, 1) enddo call g_matrixtranspose(g_p_, nex, 1, ua0, g_ua0, g_pmax_, temp4, * g_temp4, g_pmax_) call g_matrixprod(g_p_, nex, 1, temp3, g_temp3, g_pmax_, 1, nex, * temp4, g_temp4, g_pmax_, nex, nex, ua1, g_ua1, g_pmax_) C if (ncov .ge. 1) then do i = 1, ncov temp5(i, 1) = x(nex + i, 1) enddo C do i = 1, ncov do j = 1, nex do g_i_ = 1, g_p_ g_ua2(g_i_, j, i) = temp5(i, 1) * g_temp4(g_i_, 1, j) enddo ua2(j, i) = temp5(i, 1) * temp4(1, j) C-------- enddo enddo endif C do i = 1, nex do g_i_ = 1, g_p_ g_uatot(g_i_, i, 1) = g_ua0(g_i_, i, 1) enddo uatot(i, 1) = ua0(i, 1) C-------- do j = 1, nex do g_i_ = 1, g_p_ g_uatot(g_i_, i, 1 + j) = g_ua1(g_i_, j, i) enddo uatot(i, 1 + j) = ua1(j, i) C-------- enddo if (ncov .ge. 1) then do j = 1, ncov do g_i_ = 1, g_p_ g_uatot(g_i_, i, 1 + nex + j) = g_ua2(g_i_, i, j) enddo uatot(i, 1 + nex + j) = ua2(i, j) C-------- enddo endif enddo C do i = 1, 1 + nd do j = 1, nex do g_i_ = 1, g_p_ g_ua(g_i_, (i - 1) * nex + j) = g_uatot(g_i_, j, i) enddo ua((i - 1) * nex + j) = uatot(j, i) C-------- enddo enddo C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_matrixprod(g_p_, nar, nac, a, g_a, ldg_a, nbr, nbc, b *, g_b, ldg_b, ncr, ncc, c, g_c, ldg_c) C C used to compute the real matrix(or vector) product c=a*b C use variable implicit none C C on the entry variables integer nar, nac, nbr, nbc, ncr, ncc real*8 a(nar, nac), b(nbr, nbc), c(ncr, ncc) C C dummy variables integer i, j, k C integer g_pmax_ parameter (g_pmax_ = (1+nex+ncov)*nex ) integer g_i_, g_p_, ldg_c, ldg_a, ldg_b double precision g_c(ldg_c, ncr, ncc), g_a(ldg_a, nar, nac), g_b *(ldg_b, nbr, nbc) integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif if ((nac .ne. nbr) .or. (nar .ne. ncr) .or. (nbc .ne. ncc)) then print *, 'Dimension does not match' stop endif C do i = 1, nar do j = 1, ncc do g_i_ = 1, g_p_ g_c(g_i_, i, j) = 0.0d0 enddo c(i, j) = 0.0d0 C-------- do k = 1, nac do g_i_ = 1, g_p_ g_c(g_i_, i, j) = a(i, k) * g_b(g_i_, k, j) + b(k, j) * *g_a(g_i_, i, k) + g_c(g_i_, i, j) enddo c(i, j) = c(i, j) + a(i, k) * b(k, j) C-------- enddo enddo enddo C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_matrixsubtract(g_p_, nr, nc, a, b, g_b, ldg_b, c, g_c *, ldg_c) C C used to compute the real matrix(or vector) subtraction c=a-b C use variable implicit none C C on the entry variables integer nr, nc real*8 a(nr, nc), b(nr, nc), c(nr, nc) C C dummy variables integer i, j, k C integer g_pmax_ parameter (g_pmax_ = (1+nex+ncov)*nex) integer g_i_, g_p_, ldg_c, ldg_b double precision g_c(ldg_c, nr, nc), g_b(ldg_b, nr, nc) integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nr do j = 1, nc do g_i_ = 1, g_p_ g_c(g_i_, i, j) = -g_b(g_i_, i, j) enddo c(i, j) = a(i, j) - b(i, j) C-------- enddo enddo C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_matrixtranspose(g_p_, nar, nac, a, g_a, ldg_a, at, g_ *at, ldg_at) C used to transpose a real matrix(or vector) C use variable implicit none C C on the entry variables integer nar, nac real*8 a(nar, nac), at(nac, nar) C C dummy variables integer i, j C integer g_pmax_ parameter (g_pmax_ = (1+nex+ncov)*nex) integer g_i_, g_p_, ldg_at, ldg_a double precision g_at(ldg_at, nac, nar), g_a(ldg_a, nar, nac) integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nac do j = 1, nar do g_i_ = 1, g_p_ g_at(g_i_, i, j) = g_a(g_i_, j, i) enddo at(i, j) = a(j, i) C-------- enddo enddo C return end C C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine rrccox_likelihood(nd,betval,fcn_ran) use variable implicit none c variables on the entry integer nd real*8 betval(nd),fcn_ran call rrc_cum_cox(nmain,nage_main,zs_common,covar_common, * t_main,age_main,del_main,betval,fcn_ran) fcn_ran=-fcn_ran c print *, '-2Log L=', 2.0d0*fcn_ran return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine rrc_cum_cox(nmax,ntm,z,cov_vec,t,t_age,del, * betval,fcn_ran) use variable implicit none c variables on the entry integer nmax, ntm real*8 z(nmax,ntm,nex), t(nmax), del(nmax), t_age(nmax,ntm) real*8 cov_vec(nmax,ntm,0:ncov) real*8 betval(nex+ncov), fcn_ran c dummy variables integer i,j,k,j1,k1,indk,kch real*8 eps parameter(eps=1.0d-8) real*8 rrf(nmax,ntm), s0(nageobs), temp do i=1,nmax do k=1,j_main(i) if (z(i,k,1).ne.-100.0d0) then temp=0.0d0 do j1=1,nex temp=temp+betval(j1)*z(i,k,j1) enddo if (ncov.ge.1) then do j1=1,ncov temp=temp+betval(nex+j1)*dble(cov_vec(i,k,j1)) enddo endif rrf(i,k)=dexp(temp) endif enddo enddo do k=1,nfa s0(k)=0.0d0 do i=1,nmax do j=1,j_main(i) if ((t_age(i,j).eq.tf(k)) * .and.(z(i,j,1).ne.-100.0d0) * ) then s0(k)=s0(k)+rrf(i,j) endif enddo enddo enddo fcn_ran=0.0d0 do i=1,ndist if (evcnt(i).ne.0.0d0) then do j1=1,nfa if (t(ilow(i)).eq.tf(j1)) indk=j1 enddo do j=ilow(i),iup(i) if ((del(j).eq.1.0d0) * .and.(z(j,j_main(j),1).ne.-100.0d0) * )then do k1=1,nex fcn_ran=fcn_ran+betval(k1)*z(j,j_main(j),k1) enddo if (ncov.ge.1.0d0) then do k1=1,ncov fcn_ran=fcn_ran+betval(nex+k1)* * dble(cov_vec(j,j_main(j),k1)) enddo endif endif enddo fcn_ran=fcn_ran-evcnt(i)*dlog(dmax1(s0(indk),eps)) endif enddo return end ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine rrcgradient_cox(nd,betval,g_fval) use variable implicit none c variables on the entry integer nd real*8 betval(nd),g_fval(nd) c dummy variables integer g_p_, ldg_betval,ldg_fcn_ran, i parameter(g_p_=nex+ncov,ldg_betval=nex+ncov,ldg_fcn_ran=nex+ncov) real*8 g_betval(ldg_betval,nd),g_fcn_ran(ldg_fcn_ran), fcn_ran g_betval=0.0d0 do i=1,g_p_ g_betval(i,i)=1.0d0 enddo call g_rrccox_likelihood(g_p_,nd,betval, g_betval, ldg_betval, fc *n_ran, g_fcn_ran, ldg_fcn_ran) g_fval=g_fcn_ran c print *, 'Gradient=', g_fval return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_rrccox_likelihood(g_p_, nd, betval, g_betval, ldg_bet *val, fcn_ran, g_fcn_ran, ldg_fcn_ran) C use variable implicit none C variables on the entry integer nd real*8 betval(nd), fcn_ran C integer g_pmax_ parameter (g_pmax_ = nex+ncov) integer g_i_, g_p_, ldg_fcn_ran, ldg_betval double precision g_fcn_ran(ldg_fcn_ran), g_betval(ldg_betval, nd *) external g_rrc_cum_cox integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif call g_rrc_cum_cox(g_p_, nmain, nage_main, zs_common, covar_comm *on, t_main, age_main, del_main, betval, g_betval, ldg_betval, fcn_ *ran, g_fcn_ran, ldg_fcn_ran) do g_i_ = 1, g_p_ g_fcn_ran(g_i_) = -g_fcn_ran(g_i_) enddo fcn_ran = -fcn_ran C-------- C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_rrc_cum_cox(g_p_, nmax, ntm, z, cov_vec, t, t_age, de *l, betval, g_betval, ldg_betval, fcn_ran, g_fcn_ran, ldg_fcn_ran) C use variable implicit none C variables on the entry integer nmax, ntm real*8 z(nmax, ntm, nex), t(nmax), del(nmax), t_age(nmax, ntm) real*8 cov_vec(nmax, ntm, 0:ncov) real*8 betval(nex + ncov), fcn_ran C C dummy variables integer i, j, k, j1, k1, indk, kch real*8 eps parameter (eps = 1.0d-8) real*8 rrf(nmain, nage_main), s0(nageobs), temp C C integer g_pmax_ parameter (g_pmax_ = nex+ncov) integer g_i_, g_p_, ldg_betval, ldg_fcn_ran double precision d5_b, d2_v, d2_p, d1_w, d1_p, d3_v, g_temp(g_pm *ax_), g_betval(ldg_betval, nex + ncov), g_rrf(g_pmax_, nmain, nage *_main), g_s0(g_pmax_, nageobs) double precision g_fcn_ran(ldg_fcn_ran), g_d1_w(g_pmax_) save g_temp, g_rrf, g_s0, g_d1_w integer g_ehfid data g_ehfid /0/ C C if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nmax do k = 1, j_main(i) if (z(i, k, 1) .ne. (-100.0d0)) then do g_i_ = 1, g_p_ g_temp(g_i_) = 0.0d0 enddo temp = 0.0d0 C-------- do j1 = 1, nex do g_i_ = 1, g_p_ g_temp(g_i_) = z(i, k, j1) * g_betval(g_i_, j1) + g_te *mp(g_i_) enddo temp = temp + betval(j1) * z(i, k, j1) C-------- enddo if (ncov .ge. 1) then do j1 = 1, ncov do g_i_ = 1, g_p_ g_temp(g_i_) = cov_vec(i, k, j1) * g_betval(g_i_, ne *x + j1) + g_temp(g_i_) enddo temp = temp + betval(nex + j1) * cov_vec(i, k, j1) C-------- enddo endif d2_v = exp(temp) d1_p = d2_v do g_i_ = 1, g_p_ g_rrf(g_i_, i, k) = d1_p * g_temp(g_i_) enddo rrf(i, k) = d2_v C-------- endif enddo enddo C do k = 1, nfa do g_i_ = 1, g_p_ g_s0(g_i_, k) = 0.0d0 enddo s0(k) = 0.0d0 C-------- do i = 1, nmax do j = 1, j_main(i) if ((t_age(i, j) .eq. tf(k)) .and. (z(i, j, 1) .ne. (-100. *0d0))) then do g_i_ = 1, g_p_ g_s0(g_i_, k) = g_rrf(g_i_, i, j) + g_s0(g_i_, k) enddo s0(k) = s0(k) + rrf(i, j) C-------- endif enddo enddo enddo C C do g_i_ = 1, g_p_ g_fcn_ran(g_i_) = 0.0d0 enddo fcn_ran = 0.0d0 C-------- do i = 1, ndist C if (evcnt(i) .ne. 0.0d0) then do j1 = 1, nfa if (t(ilow(i)) .eq. tf(j1)) then indk = j1 endif enddo C do j = ilow(i), iup(i) if ((del(j) .eq. 1.0d0) .and. (z(j, j_main(j), 1) .ne. (-1 *00.0d0))) then do k1 = 1, nex do g_i_ = 1, g_p_ g_fcn_ran(g_i_) = z(j, j_main(j), k1) * g_betval(g_i *_, k1) + g_fcn_ran(g_i_) enddo fcn_ran = fcn_ran + betval(k1) * z(j, j_main(j), k1) C-------- enddo C if (ncov .ge. 1.0d0) then do k1 = 1, ncov do g_i_ = 1, g_p_ g_fcn_ran(g_i_) = cov_vec(j, j_main(j), k1) * g_be *tval(g_i_, nex + k1) + g_fcn_ran(g_i_) enddo fcn_ran = fcn_ran + betval(nex + k1) * cov_vec(j, j_ *main(j), k1) C-------- enddo endif endif enddo d2_v = max (s0(indk), eps) if (s0(indk) .gt. eps) then d1_p = 1.0d0 d2_p = 0.0d0 else if (s0(indk) .lt. eps) then d1_p = 0.0d0 d2_p = 1.0d0 else d1_p = 0.5d0 d2_p = 0.5d0 endif do g_i_ = 1, g_p_ g_d1_w(g_i_) = d1_p * g_s0(g_i_, indk) enddo d1_w = d2_v d3_v = log(d1_w) d1_p = 1.0d0 / d1_w d5_b = -evcnt(i) * d1_p do g_i_ = 1, g_p_ g_fcn_ran(g_i_) = d5_b * g_d1_w(g_i_) + g_fcn_ran(g_i_) enddo fcn_ran = fcn_ran - evcnt(i) * d3_v C-------- C endif C enddo C return end C C Ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine rrchessian_cox(nd,betval,hess_fval,ldh) use variable implicit none c variables on the entry integer nd,ldh real*8 betval(nd), hess_fval(ldh,nd) c dummy variables integer g__p_, ldg__betval, ldg_g_fval,i parameter(g__p_=nex+ncov,ldg__betval=nex+ncov,ldg_g_fval=nex+ncov) real*8 g__betval(ldg__betval,nd), fcn_ran, g_fval(nd) real*8 g_g_fval(ldg_g_fval, nd) g__betval=0.0d0 do i=1,g__p_ g__betval(i,i)=1.0d0 enddo call g_rrcgradient_cox(g__p_,nd,betval, g__betval, ldg__betval *, g_fval, g_g_fval, ldg_g_fval) hess_fval=g_g_fval return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_rrcgradient_cox(g__p_, nd, betval, g__betval, ldg__be *tval, g_fval, g_g_fval, ldg_g_fval) C use variable implicit none C C variables on the entry integer nd real*8 betval(nd), g_fval(nd) C C dummy variables integer g_p_, ldg_betval, ldg_fcn_ran, i, j parameter (g_p_ =nex+ncov,ldg_betval =nex+ncov) parameter (ldg_fcn_ran =nex+ncov) real*8 g_betval(ldg_betval, ldg_betval), g_fcn_ran(ldg_fcn_ran) real*8 fcn_ran C C integer g_pmax_ parameter (g_pmax_ = nex+ncov) integer g_i_, g__p_, ldg_g_fval, ldg__betval double precision g_g_fval(ldg_g_fval, nd), g_g_fcn_ran(g_pmax_, *ldg_fcn_ran), g__betval(ldg__betval, nd) save g_g_fcn_ran external g_g_rrccox_likelihood integer g_ehfid data g_ehfid /0/ C C if (g__p_ .gt. g_pmax_) then print *, 'Parameter g__p_ is greater than g_pmax_' stop endif do i = 1, g_p_ do j = 1, g_p_ if (i .ne. j) then g_betval(i, j) = 0.0d0 else g_betval(i, j) = 1.0d0 endif enddo enddo C call g_g_rrccox_likelihood(g__p_, g_p_, nd, betval, g__betval, l *dg__betval, g_betval, ldg_betval, fcn_ran, g_fcn_ran, g_g_fcn_ran, * g_pmax_, ldg_fcn_ran) C do i = 1, g_p_ do g_i_ = 1, g__p_ g_g_fval(g_i_, i) = g_g_fcn_ran(g_i_, i) enddo g_fval(i) = g_fcn_ran(i) C-------- enddo C return end C C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_g_rrccox_likelihood(g__p_, g_p_, nd, betval, g__betva *l, ldg__betval, g_betval, ldg_betval, fcn_ran, g_fcn_ran, g_g_fcn_ *ran, ldg_g_fcn_ran, ldg_fcn_ran) C use variable implicit none C C variables on the entry integer nd real*8 betval(nd), fcn_ran C integer g_pmax_ parameter (g_pmax_ = nex+ncov) integer g_i_, g_p_, ldg_fcn_ran, ldg_betval double precision g_fcn_ran(ldg_fcn_ran), g_betval(ldg_betval, nd *) integer g_ehfid integer g__pmax_ parameter (g__pmax_ = nex+ncov) integer g__i_, g__p_, ldg_g_fcn_ran, ldg__betval double precision g_g_fcn_ran(ldg_g_fcn_ran, ldg_fcn_ran), g__bet *val(ldg__betval, nd) external g_g_rrc_cum_cox data g_ehfid /0/ C C C C C if (g__p_ .gt. g__pmax_) then print *, 'Parameter g__p_ is greater than g__pmax_' stop endif if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif call g_g_rrc_cum_cox(g__p_, g_p_, nmain, nage_main, zs_common, c *ovar_common, t_main, age_main, del_main, betval, g__betval, ldg__b *etval, g_betval, ldg_betval, fcn_ran, g_fcn_ran, g_g_fcn_ran, ldg_ *g_fcn_ran, ldg_fcn_ran) do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_fcn_ran(g__i_, g_i_) = -g_g_fcn_ran(g__i_, g_i_) enddo g_fcn_ran(g_i_) = -g_fcn_ran(g_i_) C-------- enddo fcn_ran = -fcn_ran C-------- C return end C Ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss C subroutine g_g_rrc_cum_cox(g__p_, g_p_, nmax, ntm, z, cov_vec, t, *t_age, del, betval, g__betval, ldg__betval, g_betval, ldg_betval, *fcn_ran, g_fcn_ran, g_g_fcn_ran, ldg_g_fcn_ran, ldg_fcn_ran) C use variable implicit none C C C variables on the entry integer nmax, ntm real*8 z(nmax, ntm, nex), t(nmax), del(nmax), t_age(nmax, ntm) real*8 cov_vec(nmax, ntm, 0:ncov) real*8 betval(nex + ncov), fcn_ran C C dummy variables integer i, j, k, j1, k1, indk, kch real*8 eps parameter (eps = 1.0d-8) real*8 rrf(nmain, nage_main), s0(nageobs), temp C C integer g_pmax_ parameter (g_pmax_ = nex+ncov) integer g_i_, g_p_, ldg_betval, ldg_fcn_ran double precision d5_b, d2_v, d2_p, d1_w, d1_p, d3_v, g_temp(g_pm *ax_), g_betval(ldg_betval, nex + ncov), g_rrf(g_pmax_, nmain, nage *_main), g_s0(g_pmax_, nageobs) double precision g_fcn_ran(ldg_fcn_ran), g_d1_w(g_pmax_) save g_temp, g_rrf, g_s0, g_d1_w integer g_ehfid integer g__pmax_ parameter (g__pmax_ = nex+ncov) integer g__i_, g__p_, ldg__betval, ldg_g_fcn_ran double precision d2_b, d4_v, d4_p, d3_p, g__temp(g__pmax_), g__b *etval(ldg__betval, nex + ncov), g_d2_v(g__pmax_), g_d1_p(g__pmax_) *, g_g_rrf(g__pmax_, g_pmax_, nmain, nage_main), g__rrf(g__pmax_, n *main, nage_main) double precision g_g_s0(g__pmax_, g_pmax_, nageobs), g__s0(g__pm *ax_, nageobs), g_g_fcn_ran(ldg_g_fcn_ran, ldg_fcn_ran), g_g_d1_w(g *__pmax_, g_pmax_), g__d1_w(g__pmax_), g_d5_b(g__pmax_) save g__temp, g_d2_v, g_d1_p, g_g_rrf, g__rrf, g_g_s0, g__s0, g_ *g_d1_w, g__d1_w, g_d5_b data g_ehfid /0/ C C C integer g_ehfid data g_ehfid /0/ C C if (g__p_ .gt. g__pmax_) then print *, 'Parameter g__p_ is greater than g__pmax_' stop endif if (g_p_ .gt. g_pmax_) then print *, 'Parameter g_p_ is greater than g_pmax_' stop endif do i = 1, nmax do k = 1, j_main(i) if (z(i, k, 1) .ne. (-100.0d0)) then do g_i_ = 1, g_p_ g_temp(g_i_) = 0.0d0 enddo do g__i_ = 1, g__p_ g__temp(g__i_) = 0.0d0 enddo temp = 0.0d0 C-------- C-------- do j1 = 1, nex do g_i_ = 1, g_p_ g_temp(g_i_) = z(i, k, j1) * g_betval(g_i_, j1) + g_te *mp(g_i_) enddo do g__i_ = 1, g__p_ g__temp(g__i_) = z(i, k, j1) * g__betval(g__i_, j1) + *g__temp(g__i_) enddo temp = temp + betval(j1) * z(i, k, j1) C-------- C-------- enddo if (ncov .ge. 1) then do j1 = 1, ncov do g_i_ = 1, g_p_ g_temp(g_i_) = cov_vec(i, k, j1) * g_betval(g_i_, ne *x + j1) + g_temp(g_i_) enddo do g__i_ = 1, g__p_ g__temp(g__i_) = cov_vec(i, k, j1) * g__betval(g__i_ *, nex + j1) + g__temp(g__i_) enddo temp = temp + betval(nex + j1) * cov_vec(i, k, j1) C-------- C-------- enddo endif d4_v = exp(temp) d3_p = d4_v do g__i_ = 1, g__p_ g_d2_v(g__i_) = d3_p * g__temp(g__i_) enddo d2_v = d4_v C-------- do g__i_ = 1, g__p_ g_d1_p(g__i_) = g_d2_v(g__i_) enddo d1_p = d2_v C-------- do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_rrf(g__i_, g_i_, i, k) = g_temp(g_i_) * g_d1_p(g__ *i_) enddo g_rrf(g_i_, i, k) = d1_p * g_temp(g_i_) C-------- enddo do g__i_ = 1, g__p_ g__rrf(g__i_, i, k) = g_d2_v(g__i_) enddo rrf(i, k) = d2_v C-------- C-------- endif enddo enddo C do k = 1, nfa do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_s0(g__i_, g_i_, k) = 0.0d0 enddo g_s0(g_i_, k) = 0.0d0 C-------- enddo do g__i_ = 1, g__p_ g__s0(g__i_, k) = 0.0d0 enddo s0(k) = 0.0d0 C-------- C-------- do i = 1, nmax do j = 1, j_main(i) if ((t_age(i, j) .eq. tf(k)) .and. (z(i, j, 1) .ne. (-100. *0d0))) then do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_s0(g__i_, g_i_, k) = g_g_s0(g__i_, g_i_, k) + g_ *g_rrf(g__i_, g_i_, i, j) enddo g_s0(g_i_, k) = g_rrf(g_i_, i, j) + g_s0(g_i_, k) C-------- enddo do g__i_ = 1, g__p_ g__s0(g__i_, k) = g__rrf(g__i_, i, j) + g__s0(g__i_, k *) enddo s0(k) = s0(k) + rrf(i, j) C-------- C-------- endif enddo enddo enddo C C do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_fcn_ran(g__i_, g_i_) = 0.0d0 enddo g_fcn_ran(g_i_) = 0.0d0 C-------- enddo fcn_ran = 0.0d0 C-------- do i = 1, ndist C if (evcnt(i) .ne. 0.0d0) then do j1 = 1, nfa if (t(ilow(i)) .eq. tf(j1)) then indk = j1 endif enddo C do j = ilow(i), iup(i) if ((del(j) .eq. 1.0d0) .and. (z(j, j_main(j), 1) .ne. (-1 *00.0d0))) then do k1 = 1, nex do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_fcn_ran(g__i_, g_i_) = g_g_fcn_ran(g__i_, g_i_ *) enddo g_fcn_ran(g_i_) = z(j, j_main(j), k1) * g_betval(g_i *_, k1) + g_fcn_ran(g_i_) C-------- enddo fcn_ran = fcn_ran + betval(k1) * z(j, j_main(j), k1) C-------- enddo C if (ncov .ge. 1.0d0) then do k1 = 1, ncov do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_fcn_ran(g__i_, g_i_) = g_g_fcn_ran(g__i_, g_ *i_) enddo g_fcn_ran(g_i_) = cov_vec(j, j_main(j), k1) * g_be *tval(g_i_, nex + k1) + g_fcn_ran(g_i_) C-------- enddo fcn_ran = fcn_ran + betval(nex + k1) * cov_vec(j, j_ *main(j), k1) C-------- enddo endif endif enddo d4_v = max (s0(indk), eps) if (s0(indk) .gt. eps) then d3_p = 1.0d0 d4_p = 0.0d0 else if (s0(indk) .lt. eps) then d3_p = 0.0d0 d4_p = 1.0d0 else d3_p = 0.5d0 d4_p = 0.5d0 endif do g__i_ = 1, g__p_ g_d2_v(g__i_) = d3_p * g__s0(g__i_, indk) enddo d2_v = d4_v C-------- if (s0(indk) .gt. eps) then do g__i_ = 1, g__p_ g_d1_p(g__i_) = 0.0d0 enddo d1_p = 1.0d0 C-------- d2_p = 0.0d0 else if (s0(indk) .lt. eps) then do g__i_ = 1, g__p_ g_d1_p(g__i_) = 0.0d0 enddo d1_p = 0.0d0 C-------- d2_p = 1.0d0 else do g__i_ = 1, g__p_ g_d1_p(g__i_) = 0.0d0 enddo d1_p = 0.5d0 C-------- d2_p = 0.5d0 endif endif do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_d1_w(g__i_, g_i_) = d1_p * g_g_s0(g__i_, g_i_, indk) * + g_s0(g_i_, indk) * g_d1_p(g__i_) enddo g_d1_w(g_i_) = d1_p * g_s0(g_i_, indk) C-------- enddo do g__i_ = 1, g__p_ g__d1_w(g__i_) = g_d2_v(g__i_) enddo d1_w = d2_v C-------- d3_v = log(d1_w) d4_v = 1.0d0 / d1_w d2_b = -d4_v / d1_w do g__i_ = 1, g__p_ g_d1_p(g__i_) = d2_b * g__d1_w(g__i_) enddo d1_p = d4_v C-------- do g__i_ = 1, g__p_ g_d5_b(g__i_) = -evcnt(i) * g_d1_p(g__i_) enddo d5_b = -evcnt(i) * d1_p C-------- do g_i_ = 1, g_p_ do g__i_ = 1, g__p_ g_g_fcn_ran(g__i_, g_i_) = g_g_fcn_ran(g__i_, g_i_) + d5 *_b * g_g_d1_w(g__i_, g_i_) + g_d1_w(g_i_) * g_d5_b(g__i_) enddo g_fcn_ran(g_i_) = d5_b * g_d1_w(g_i_) + g_fcn_ran(g_i_) C-------- enddo fcn_ran = fcn_ran - evcnt(i) * d3_v C-------- C endif C enddo C return end C C Ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine mvlm2(p,q,n,y,xin,gamma) c c this routine computes the parameter estimates and their variance-covariance c matrix for the multivariate linear model of p-vector y on (p+q)-vector xin c the slopes are in (p+q)xp matrix g; aprim has the intercepts c c n = number of observations in the dataset c c p = dim(Y) ! which are the true covariates c c q = dim(X) ! which are the error measured covariates and other explanatory variable c integer p,q,q1 real*8 y(n,p),xin(n,q),g(q,p),aprime(p),x(n,q+1), 1 xtx(q+1,q+1),xix(q+1,n),gamma(q+1,p),vargam(p*q,p*q),e(n,n), 2 proj(n,n),yp(p,n),s(p,p),xtx2(q,q),temp(q,q), 3 isubn(n,n) c c add 1 to the dimension of the problem to account for intercepts c c augment x matrix with a vector of 1's c q1=q+1 do 10 i=1,n do 11 j=1,q x(i,j+1)=xin(i,j) 11 continue x(i,1)=1.d0 10 continue c c get means of x,y c c do 83 j=1,q+1 c xb(j)=0.d0 c83 continue c do 82 i=1,n c do 82 j=1,q+1 c xb(j)=xb(j)+x(i,j) c82 continue c do 87 j=1,q+1 c xb(j)=xb(j)/float(n) c87 continue c write(*,*) ' means of x ',xb c c compute gammas c c first, get x'x c call dmxtxf(n,q1,x,n,q1,xtx,q1) c c get inverse (x'x)-1 c call dlinrg(q1,xtx,q1,xtx,q1) c c (x'x)-1 x' c call dmxytf(q1,q1,xtx,q1,n,q1,x,n,q1,n,xix,q1) c c (x'x)-1 x' y c call dmrrrr(q1,n,xix,q1,n,p,y,n,q1,p,gamma,q1) c return end Ccccccsssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine data_summary c this subroutine is used to summarize the data statistics. use variable implicit none integer i,nd parameter (nd=nex+ncov) real*8 naivci_left(nd), naivci_right(nd), or_naiv(nd) real*8 rrcci_left(nd), rrcci_right(nd), or_rrc(nd), pval(nd) character*7 names(nd) open(unit=outfile1,file='Output.txt',status='replace') write(outfile1,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' write(outfile1,*) '' write(outfile1,*) 'The number of risk sets:', ngrp write(outfile1,*) '# of subjects in the main study:', nmain write(outfile1,*) '# of person-year observations in the main', * ' study:', mainobs write(outfile1,*) '# of prostate cancer cases in the main study:', * nint(sum(del_main)) write(outfile1,*) '# of subjects in the validation study:',nval write(outfile1,*) '# of person-year observations in the ', * 'validation study:', valobs write(outfile1,*) '# of prostate cancer cases in the validation', * ' study:',nint(sum(del_val)) write(outfile1,*) '# of unique failure time:', nfa if ((nex.eq.1).and.(ncov.eq.1)) then write(outfile1,*) 'The Vit-E is a error-free covariate,' write(outfile1,*) 'The total calcium is measured with error.' else if ((nex.eq.2).and.(ncov.eq.0)) then write(outfile1,*) 'Both the Vit-E and calcium are measured *with error.' endif write(outfile1,*) '' names(1)='Calcium' if ((nex.gt.1).or.(ncov.ge.1)) then names(nex+ncov)='Tot VE' endif write(outfile1,*) '---------------------------------------------' do i=1,nd call CI_95_alco(betaes_naiv(i), serror_naiv(i),or_naiv(i), * naivci_left(i), naivci_right(i), pval(i)) enddo write(outfile1,*) 'Naive cox result:' write(outfile1,90) 'Variable','Para_es','S.E.','p_val','H.R.', * '95% C.I.' do i=1,nd write(outfile1,100) names(i), betaes_naiv(i),serror_naiv(i), * pval(i),or_naiv(i),'[',naivci_left(i),',',naivci_right(i),']' enddo write(outfile1,*) '---------------------------------------------' do i=1,nd call CI_95_alco(betaes_rrc(i), serror_rrc(i), or_rrc(i), * rrcci_left(i),rrcci_right(i), pval(i)) enddo write(outfile1,*) 'RRC result:' write(outfile1,90) 'Variable','Para_es','S.E.','p_val','H.R.', * '95% C.I.' do i=1,nd write(outfile1,100) names(i), betaes_rrc(i),serror_rrc(i), * pval(i),or_rrc(i),'[',rrcci_left(i),',',rrcci_right(i),']' enddo write(outfile1,*) '---------------------------------------------' write(outfile1,*) '' write(outfile1,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' 90 format(1X,A8,2X,A7,1X,A6,2X,A5,2X,A5,2X,A14) 100 format(1X,A8,2X,F6.3,2X,F6.3,2X,F5.3,2X,F5.3,2X, A, F6.2, A, F6.2, * A) return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine CI_95_alco(esti, serror, oddr, in_left, in_right, pval) c this subroutine is used to calculate the coverage probability c and the relative percentage bias. implicit none c variables on the entry real*8 esti, serror,oddr,in_left, in_right, pval c dummy variables real*8 cr, scale, dchidf cr=1.960d0 scale=1.0d0 oddr=dexp(esti*scale) in_left=dexp((esti-cr*serror)*scale) in_right=dexp((esti+cr*serror)*scale) pval=1.0d0-dchidf(esti**2/serror**2,1.0d0) return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss subroutine cal_cov(ldx,x1,x2,x1mean,x2mean,corr, nmiss) use variable implicit none c variables for IMSL subroutine corvc c------------------------------------------------------------- integer i integer ido,nrow,nvar,ldx,ifrq,iwt,mopt,icopt,ldcov integer ldincd,nobs,nmiss parameter (nvar=2,ldcov=2,ldincd=1) real*8 x(ldx,nvar),cov(ldcov,nvar),sumwt,xmean(nvar) integer incd(nvar,nvar) c variables on the entry real*8 x1(ldx), x2(ldx), x1mean, x2mean, corr external dcorvc c------------------------------------------------------------- ido=0 nrow=ldx ifrq=0 iwt=0 mopt=3 c COV contains the correlation matrix, except c for the diagonal elements, which are the standard deviations. icopt=3 do i=1,ldx x(i,1)=x1(i) x(i,2)=x2(i) enddo call dcorvc(ido,nrow,nvar,x,ldx,ifrq,iwt,mopt,icopt,xmean,cov, * ldcov,incd,ldincd,nobs,nmiss,sumwt) x1mean=xmean(1) x2mean=xmean(2) corr=cov(1,2) return end ccccccssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss