C-------------------------------------------------------------------------- C PROGRAM: a5055.for C DYNAMIC MODEL: C dT/dt = lambda-rho*T-[1-gamma(t)]*k*T*V C dTp/dt = [1-gamma(t)]*k*T*V - delta*Tp C dV/dt = N*delta*Tp - c*V C FORM/INPUT(S): hyper-parameters: alpha,beta,eta,big_GAMMA, big_OMAGA,nu C C12h, IC50(t), A(t) gamma(t), C Data--Viral loads, IC50, PK and adherence C OUTPUT(S): estimate parameters: phi,c, delta,lambda,rho,N,k for both po[ulation and C individuals C PROTOCOL: 1. THE Runge-Kutta-Verner METHOD OF DIVPRK IS USED C TO SOLVE DIFFERENTIAL EQUATIONS C 2. TO ESTIMATE PARAMETERS USING BAYESIAN METHOD VIA MCMC C Original program: 07/05/02 Y.Huang C Modification history: 07/05/05 Y.Huang C---------------------------------------------------------------------------- C **************************************************** C *(1) MAIN PROGRAM C **************************************************** program A5055 implicit double precision (a-h, o-z) integer i,j,k,n_sub,n_obs,np,lda,ldainv,nra,nca integer nx,ny,ipath,nr,iseed parameter(n_sub=42,n_obs=11,n_par=7, nn=328, + n_tot=1+n_par+(n_par*(1+n_par)/2)) parameter(n_obs10=10,n_obs9=9,n_obs8=8) parameter(n_obs7=7,n_obs6=6,n_obs5=5,n_obs4=4,n_obs3=3) parameter(nmcmc=100000,nsmcmc=20000,nk=4, + length=(nmcmc-nsmcmc)/nk) dimension m(n_sub) dimension x(n_sub, n_obs),z(n_sub,n_obs),tp(n_obs+1,n_sub) dimension a_mu(n_par),smatinv(n_par,n_par),theta(n_sub,n_par) dimension eta(n_par),a_laminv(n_par,n_par),omeginv(n_par,n_par) dimension c(n_par),cc(n_par),cct(n_par),tvc(n_par) dimension d(n_par,n_par),e(n_par,n_par),estexp(n_par) dimension a_mulower(n_par),a_muupper(n_par),a_muexp(n_par) dimension a_mumedian(n_par),sd(n_par),sdlog(n_par),vinlog(n_par) dimension t_mu(length,n_par),tmu1(length),tmu2(length) dimension dinv(n_par,n_par),einv(n_par,n_par), AD(2) dimension thetap(n_par),propcov(n_par,n_par),ss(n_sub) dimension tmp(n_tot),tmpexp(n_par),exptmp(n_tot),est(n_tot) dimension Aidv_Cmin(n_sub),Aidv_IC1(n_sub),Aidv_IC2(n_sub), + Ridv(n_sub) dimension Artv_Cmin(n_sub),Artv_IC1(n_sub),Artv_IC2(n_sub), + Rrtv(n_sub) dimension theta_sum(n_sub,n_par) dimension thetaexp(n_sub,n_par) dimension theta_exp(n_sub,n_par) dimension est_exp(n_sub,n_par) dimension t_theta(n_sub,length,n_par) dimension atheta(n_sub,length) dimension a_theta(length), b_theta(length) dimension sd_s(n_sub,n_par) dimension sdlog_s(n_sub,n_par) dimension vinlog_s(n_sub,n_par) dimension a_me(n_sub,n_par) dimension a_l(n_sub,n_par) dimension a_u(n_sub,n_par) EXTERNAL RNUND,DSCAL,DADD,DRNNOA external gamma,wish3,mvnor3,dlinrg,dmurrv,dblinf,drnun,rnset, + dsvrgn,drnunf C The values of hyper-parameters parameter (alpha=4.5,beta=9.0,a_nu=8.0) data eta/1.5,1.1,-1.0,4.6,-2.5,6.9,-11.0/ data a_laminv/0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 400, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 400, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 400, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 400, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 400, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1000/ data omeginv/0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5/ C Data data m/n_obs8,n_obs9,n_obs4,n_obs9,n_obs6,n_obs8,n_obs3, + n_obs8,n_obs9,n_obs8,n_obs9,n_obs9,n_obs9,n_obs9, + n_obs9,n_obs9,n_obs6,n_obs8,n_obs, n_obs10,n_obs10, + n_obs5,n_obs8,n_obs9,n_obs7,n_obs9,n_obs6,n_obs9, + n_obs6,n_obs3,n_obs6,n_obs8,n_obs9,n_obs7,n_obs9, + n_obs7,n_obs8,n_obs10,n_obs5,n_obs10,n_obs6,n_obs10/ data propcov/0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.01,0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.01,0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.01,0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.01,0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01/ C PARAMETER (MXPARM=50, N=3) PARAMETER(MABSE=1,MNS=3000000) DIMENSION PARAM(MXPARM), Y(N) DIMENSION V(n_sub, n_obs+1), A_V0(n_sub),A_PHI(n_sub),A_C0(n_sub) DIMENSION A_NUM(n_sub),A_DELTA(n_sub),A_LAMBDA(n_sub),A_AK(n_sub) DIMENSION A_RHO(n_sub) EXTERNAL DIVPRK,SSET,UMACH,FCN COMMON /COM/PHI,C0,DELTA,ALAMBDA,RHO,ANUM,AK COMMON /COM1/Cmin_idv, C1_50idv, C2_50idv, AR50idv COMMON /COM2/Cmin_rtv, C1_50rtv, C2_50rtv, AR50rtv COMMON /COM3/ad_idv, ad_rtv C PK and IC50 data for two drugs C --- Ridv is individual failure time point data Ridv/ 105.0, 0.0, 0.0, 88.0, 0.0, 56.0, 0.0, + 0.0, 0.0, 105.0, 0.0, 56.0,0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 77.0, 56.0, 0.0, + 28.0, 0.0, 78.0, 0.0, 28.0, 0.0, 0.0, + 71.0, 0.0, 28.0, 0.0, 0.0, 0.0, 86.0, + 0.0, 0.0, 54.0, 0.0, 0.0, 0.0/ C --- Aidv_Cmin is c12h data of IDV drug--- data Aidv_Cmin/1509.0,1828.4,539.9,529.1,80.4,614.0,770.7, + 658.9,633.8,2674.5,422.2,828.9,871.9,105.5, + 601.2,1775.2,1037.3,1909.6,744.6,2918.8,648.3,305.4, + 655.5,566.8,1451.4,113.5,818.5,343.7,1497.1,548.7, + 664.7,164.8, 854.3,1561.0,144.3,793.3,221.6, + 739.8,679.1,218.5,620.8,332.9/ C --- IC50 change at baseline and failure time for IDV drug--- data Aidv_IC1/2.4,9.6,7.9,4.2,2.5,2.8,6.0,2.8,7.5,2.9,4.3, + 2.3,9.9,3.3,7.9,2.5,7.1,11.8,3.1,2.6,2.0,2.8, + 2.58,3.31,2.76,1.96,3.44,9.63,3.87,4.42,18.53,5.46, + 4.91,3.31,8.71,4.97,5.95,10.12,2.88,10.0,2.45,8.10/ data Aidv_IC2/2.4,9.6,7.9,5.6,2.5,3.3,6.0,2.8,7.5,3.1,4.3, + 1.3,9.9,3.3,7.9,2.5,7.1,11.8,3.1,1.7,2.4,2.8, + 3.68,3.31,2.58,1.96,1.90, + 9.63,3.87,4.05,18.53,3.80, + 4.91,3.31,8.71,11.96,5.95, + 10.12,1.66,10.0,2.45,8.10/ C --- Ridv is individual failure time point data Rrtv/105.0, 0.0, 0.0, 88.0, 0.0, 56.0, 0.0, + 0.0, 0.0, 105.0, 0.0, 56.0,0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 77.0, 56.0, 0.0, + 28.0, 0.0, 78.0, 0.0, 28.0, 0.0, 0.0, + 71.0, 0.0, 28.0, 0.0, 0.0, 0.0, 86.0, + 0.0, 0.0, 54.0, 0.0, 0.0, 0.0/ C --- Artv_Cmin is c12h data of RTV drug--- data Artv_Cmin/1160.3,2964.9,854.3,806.2,666.8,399.4,670.0, + 649.2,841.8,2515.9,761.4,2124.7,2714.1,581.3, + 448.5,4627.3,1462.7,1067.0,499.5,1727.2,971.2,189.5, + 3264.2,2511.1,6107.5,968.6,2539.6, + 785.9,6043.6,2322.7,2359.3,655.4, + 3681.5,5521.0,258.4,2978.6,556.0, + 2616.0,3188.6,2547.0,3017.7,1491.9/ C --- IC50 change at baseline and failure time for RTV drug--- data Artv_IC1/9.8,13.9,11.8,17.7,9.6,11.2,28.2,11.8,16.0,11.9,7.5, + 8.0,8.7,7.7,37.7,97.1,44.5,222.4,7.7,20.0,3.5,14.6, + 8.85,6.69,5.04,4.82,12.09,25.18,13.74,6.91,75.68, + 13.02,6.26,7.70,26.62,22.66,20.58,35.54,5.68,49.06,7.70,17.19/ data Artv_IC2/9.8,13.9,11.8,37.5,9.6,15.8,28.2,11.8,16.0,14.7,7.5, + 5.8,8.7,7.7,37.7,97.1,44.5,222.4,7.7,10.0,7.1,14.6, + 12.81,6.69,9.28,4.82,6.12,25.18,13.74,13.60,75.68,16.26, + 6.26,7.70,26.62,55.04,20.58,35.54,5.54,49.06,7.70,17.19/ open (unit=9,file="a5055.out") open (unit=4,file="a5055.dat") open (unit=5,file="a5055.par.dat") open (unit=6,file="a5055.mu.dat") open (unit=7,file="a5055.sample.dat") open (unit=2,file="V42.dat") open (unit=3,file="day42.dat") C Read rna and time data DO 2 i=1,n_sub READ(2,*) (x(i,j),j=1,n_obs) 2 CONTINUE DO 3 i=1,n_obs+1 READ(3,*) (tp(i,j),j=1,n_sub) 3 CONTINUE ISEED=12344 CALL RNSET(ISEED) CALL UMACH(2,NOUT) c Setup firt value do 4 i=1,n_tot tmp(i)=0.0 tmpexp(i)=0.0 4 continue do 41 i=1,n_par do 42 l=1,n_sub theta_sum(l,i)=0.0 thetaexp(l,i)=0.0 42 continue 41 continue do 7 i=1,n_sub do 6 j=1,m(i) z(i,j)=log10(x(i,j)) 6 continue 7 continue C Specifications for optiona in ODE TOL=0.01 CALL SSET(MXPARM,0.0,PARAM,1) PARAM(10)=MABSE PARAM(4)=MNS C initialize DO 61 k=1,n_sub A_V0(k)=x(k,1) 61 CONTINUE c data A_PHI/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, c + 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/ c data A_PHI/12,12,12,12,12,12,12,12,12,12,12,12,12,12,12, c + 12,12,12,12,12,12,12,12,12,12,12,12,12,12,12, c + 12,12,12,12,12,12,12,12,12,12,12,12/ data A_PHI/4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5, + 4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5, + 4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5, + 4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5/ data A_C0/3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3, + 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3/ data A_DELTA/0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37, + 0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37, + 0.37,0.37,0.37,0.37,0.37,3.7, + 0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37, + 0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37, + 0.37,0.37,0.37,0.37/ data A_LAMBDA/100,100,100,100,100,100,100,100, + 100,100,100,100,100,100,100,100, + 100,100,100,100,100,100, + 100,100,100,100,100,100,100,100,100,100, + 100,100,100,100,100,100,100,100,100,100/ data A_RHO/0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08, + 0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08, + 0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08, + 0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08, + 0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08, + 0.08,0.08/ data A_AK/0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001,0.00001,0.00001,0.00001, + 0.00001,0.00001/ do 8 i=1,n_sub theta(i,1)=log(A_PHI(i)) theta(i,2)=log(A_C0(i)) theta(i,3)=log(A_DELTA(i)) theta(i,4)=log(A_LAMBDA(i)) theta(i,5)=log(A_RHO(i)) theta(i,6)=log((A_V0(i)+A_RHO(i)/A_AK(i))* + (A_C0(i)/A_LAMBDA(i))) theta(i,7)=log(A_AK(i)) 8 continue siginv=1.0 do 10 i=1,n_par do 9 j=1,n_par smatinv(i,j)=0.0 9 continue 10 continue smatinv(1,1)=5.0 smatinv(2,2)=5.0 smatinv(3,3)=5.0 smatinv(4,4)=0.5 smatinv(5,5)=1.0 smatinv(6,6)=1.0 smatinv(7,7)=0.5 a_mu(1)=0.0 a_mu(2)=1.0 a_mu(3)=-1.0 a_mu(4)=6.5 a_mu(5)=-3.5 a_mu(6)=4.5 a_mu(7)=-12 acc=0.0 n1=0 n2=0 n3=0 Vlimit=25.0 C Start MCMC do 999 imcmc=1,nmcmc a0=alpha+nn/2 99 continue acc=acc-n3 n3=0 do 13 i=1,n_sub ss(i)=0.0 PHI=exp(theta(i,1)) C0=exp(theta(i,2)) DELTA=exp(theta(i,3)) ALAMBDA=exp(theta(i,4)) RHO=exp(theta(i,5)) AK=exp(theta(i,7)) ANUM=exp(theta(i,6)) c print *, 'i and First est par is:', i c write(*,'(7E12.5)') PHI,C0,DELTA,ALAMBDA, RHO,ANUM,AK Cmin_idv=Aidv_Cmin(i) C1_50idv=Aidv_IC1(i) C2_50idv=Aidv_IC2(i) AR50idv=Ridv(i) Cmin_rtv=Artv_Cmin(i) C1_50rtv=Artv_IC1(i) C2_50rtv=Artv_IC2(i) AR50rtv=Rrtv(i) V0=ALAMBDA*ANUM/C0-RHO/AK IF (V0 .LE. 0.0) THEN n1=n1+1 END IF T0=C0/(AK*ANUM) Tp0=C0*V0/(DELTA*ANUM) T=tp(1,i) Y(1)=T0 Y(2)=Tp0 Y(3)=V0 V(i,1)=log10(Y(3)) istep=1 ss(i)=ss(i)+(z(i,1)-V(i,1))**2 IDO=1 12 CONTINUE istep=istep+1 j=istep TEND=tp(j,i) C ****** CALL SUBROUTINE FOR ADHERENCE***** CALL ADHEREN(i, TEND, AD) ad_idv=AD(1) ad_rtv=AD(2) C ****** CALL IMSL SUBROUTINE to solve ODE CALL DIVPRK(IDO,N,FCN,T,TEND,TOL,PARAM,Y) IF (Y(3) .LT. Vlimit) THEN Y(3)=Vlimit END IF V(i,j)=log10(Y(3)) IF (tp(j,i) .LE. tp(m(i),i)) THEN IF (tp(j,i) .EQ. tp(m(i),i)) IDO=3 ss(i)=ss(i)+(z(i,j)-V(i,j))**2 GO TO 12 END IF 13 continue sum_ssi=0 do 133 i=1,n_sub sum_ssi=sum_ssi+ss(i) 133 continue b=1.0/(1.0/beta+0.5*sum_ssi) call gamma(siginv,a0,b) do 15 i=1,n_par do 14 j=1,n_par dinv(i,j)=dble(n_sub)*smatinv(i,j)+a_laminv(i,j) 14 continue 15 continue np=n_par lda=n_par ldainv=n_par call dlinrg(np,dinv,lda,d,ldainv) nra=n_par nca=n_par lda=n_par nx=n_par ny=n_par ipath=2 call dmurrv(nra,nca,a_laminv,lda,nx,eta,ipath,ny,cc) do 18 i=1,n_sub do 16 j=1,n_par tvc(j)=theta(i,j) 16 continue call dmurrv(nra,nca,smatinv,lda,nx,tvc,ipath,ny,cct) do 17 j=1,n_par cc(j)=cc(j)+cct(j) 17 continue 18 continue call dmurrv(nra,nca,d,lda,nx,cc,ipath,ny,c) call mvnor3(a_mu,c,d) f=dble(n_sub)+a_nu do 21 i=1,n_par do 20 j=1,n_par einv(i,j)=omeginv(i,j) do 19 k=1,n_sub einv(i,j)=einv(i,j)+(theta(k,i)-a_mu(i))* + (theta(k,j)-a_mu(j)) 19 continue 20 continue 21 continue call dlinrg(np,einv,lda,e,ldainv) call wish3(smatinv,e,f) do 510 i=1,n_sub do 520 j=1,n_par tvc(j)=theta(i,j) 520 continue call mvnor3(thetap,tvc,propcov) prob=0.5*siginv*ss(i) do 530 j=1,n_par tvc(j)=theta(i,j)-a_mu(j) 530 continue prob=prob+0.5*dblinf(nra,nca,smatinv,lda,tvc,tvc) do 540 j=1,n_par tvc(j)=thetap(j)-a_mu(j) 540 continue prob=prob-0.5*dblinf(nra,nca,smatinv,lda,tvc,tvc) C A_PHI(i)=exp(thetap(1)) A_C0(i)=exp(thetap(2)) A_DELTA(i)=exp(thetap(3)) A_LAMBDA(i)=exp(thetap(4)) A_RHO(i)=exp(thetap(5)) A_NUM(i)=exp(thetap(6)) A_AK(i)=exp(thetap(7)) PHI=A_PHI(i) C0=A_C0(i) DELTA=A_DELTA(i) ALAMBDA=A_LAMBDA(i) RHO=A_RHO(i) ANUM=A_NUM(i) AK=A_AK(i) c print *, 'i and Second est par is', i c write(*,'(7E12.5)') PHI,C0,DELTA,ALAMBDA, RHO,ANUM,AK Cmin_idv=Aidv_Cmin(i) C1_50idv=Aidv_IC1(i) C2_50idv=Aidv_IC2(i) AR50idv=Ridv(i) Cmin_rtv=Artv_Cmin(i) C1_50rtv=Artv_IC1(i) C2_50rtv=Artv_IC2(i) AR50rtv=Rrtv(i) V0=ALAMBDA*ANUM/C0-RHO/AK IF (V0 .LE. 0.0) THEN C V0=A_V0(i) n2=n2+1 END IF IF (V0 .LE. 1.0) GO TO 99 T0=C0/(AK*ANUM) Tp0=C0*V0/(DELTA*ANUM) T=tp(1,i) Y(1)=T0 Y(2)=Tp0 Y(3)=V0 V(i,1)=log10(Y(3)) istep=1 prob=prob-0.5*siginv*(z(i,1)-V(i,1))**2 IDO=1 550 CONTINUE istep=istep+1 j=istep TEND=tp(j,i) C ****** CALL SUBROUTINE FOR ADHERENCE***** CALL ADHEREN(i, TEND, AD) ad_idv=AD(1) ad_rtv=AD(2) C ****** CALL IMSL SUBROUTINE to solve ODE CALL DIVPRK(IDO,N,FCN,T,TEND,TOL,PARAM,Y) IF (Y(3) .LT. Vlimit) THEN Y(3)=Vlimit END IF V(i,j)=log10(Y(3)) IF (tp(j,i) .LE. tp(m(i),i)) THEN IF (tp(j,i) .EQ. tp(m(i),i)) IDO=3 prob=prob-0.5*siginv*(z(i,j)-V(i,j))**2 GO TO 550 END IF prob=exp(prob) nr=1 call drnun(nr,uni) if (uni .le. prob) then do 560 j=1,n_par theta(i,j)=thetap(j) 560 continue acc=acc+1.0 n3=n3+1 endif 510 continue n3=0 if (imcmc .gt. nsmcmc .and. mod(imcmc,nk) .eq. 0) then do 721 i=1,n_par tmp(i)=tmp(i)+a_mu(i) tmpexp(i)=tmpexp(i)+exp(a_mu(i)) a_muexp(i)=exp(a_mu(i)) t_mu((imcmc-nsmcmc)/nk,i)=exp(a_mu(i)) 721 continue tmp(8)=tmp(8)+smatinv(1,1) tmp(9)=tmp(9)+smatinv(1,2) tmp(10)=tmp(10)+smatinv(1,3) tmp(11)=tmp(11)+smatinv(1,4) tmp(12)=tmp(12)+smatinv(1,5) tmp(13)=tmp(13)+smatinv(1,6) tmp(14)=tmp(14)+smatinv(1,7) tmp(15)=tmp(15)+smatinv(2,2) tmp(16)=tmp(16)+smatinv(2,3) tmp(17)=tmp(17)+smatinv(2,4) tmp(18)=tmp(18)+smatinv(2,5) tmp(19)=tmp(19)+smatinv(2,6) tmp(20)=tmp(20)+smatinv(2,7) tmp(21)=tmp(21)+smatinv(3,3) tmp(22)=tmp(22)+smatinv(3,4) tmp(23)=tmp(23)+smatinv(3,5) tmp(24)=tmp(24)+smatinv(3,6) tmp(25)=tmp(25)+smatinv(3,7) tmp(26)=tmp(26)+smatinv(4,4) tmp(27)=tmp(27)+smatinv(4,5) tmp(28)=tmp(28)+smatinv(4,6) tmp(29)=tmp(29)+smatinv(4,7) tmp(30)=tmp(30)+smatinv(5,5) tmp(31)=tmp(31)+smatinv(5,6) tmp(32)=tmp(32)+smatinv(5,7) tmp(33)=tmp(33)+smatinv(6,6) tmp(34)=tmp(34)+smatinv(6,7) tmp(35)=tmp(35)+smatinv(7,7) tmp(36)=tmp(36)+siginv write(6, 881) (a_mu(k),k=1,n_par) do 741 i=1,n_par ndiff=(imcmc-nsmcmc)/nk est(i)=tmp(i)/ndiff 741 continue write(7, 882) imcmc-nsmcmc, (est(k), k=1,n_par) C *******ADD3******* do 7211 j=1,n_par do 7212 l=1,n_sub theta_sum(l,j)=theta_sum(l,j)+theta(l,j) thetaexp(l,j)=thetaexp(l,j)+exp(theta(l,j)) t_theta(l,(imcmc-nsmcmc)/nk,j)=exp(theta(l,j)) 7212 continue 7211 continue endif c if (mod(imcmc,1000) .eq. 0) print *,imcmc 999 continue n_lower1=length*25/1000 n_lower2=n_lower1+1 n_upper1=length*975/1000 n_upper2=n_upper1+1 do 805 j=1,n_tot tmp(j)=tmp(j)/length 805 continue do 8050 j=1,n_par tmpexp(j)=tmpexp(j)/length exptmp(j)=exp(tmp(j)) 8050 continue do 806 i=1,n_par sd(i)=0.0 sdlog(i)=0.0 do 807 k=1,length tmu1(k)=t_mu(k,i) sd(i)=sd(i)+(t_mu(k,i)-tmpexp(i))**2 sdlog(i)=sdlog(i)+(log(t_mu(k,i))-tmp(i))**2 807 continue a_mumedian(i)=t_mu(int(length/2),i) sd(i)=sqrt(sd(i)/(length-1)) vinlog(i)=1/(sdlog(i)/(length-1)) sdlog(i)=sqrt(sdlog(i)/(length-1)) call dsvrgn(length,tmu1,tmu2) a_mulower(i)=0.5*(tmu2(n_lower1)+tmu2(n_lower2)) a_muupper(i)=0.5*(tmu2(n_upper1)+tmu2(n_upper2)) estexp(i)=exp(tmp(i)+0.5*sdlog(i)**2) 806 continue C ******ADD4******* do 8051 j=1,n_par do 8052 l=1,n_sub theta_sum(l,j)=theta_sum(l,j)/length theta_exp(l,j)=exp(theta_sum(l,j)) thetaexp(l,j)=thetaexp(l,j)/length 8052 continue 8051 continue do 8061 j=1,n_par do 8062 l=1,n_sub sd_s(l,j)=0.0 sdlog_s(l,j)=0.0 do 8071 k=1,length atheta(l,k)=t_theta(l,k,j) sd_s(l,j)=sd_s(l,j)+(atheta(l,k)-thetaexp(l,j))**2 sdlog_s(l,j)=sdlog_s(l,j)+(log(atheta(l,k))- + theta_sum(l,j))**2 8071 continue vinlog_s(l,j)=1/(sdlog_s(l,j)/(length-1)) a_me(l,j)=atheta(l,int(length/2)) sd_s(l,j)=sqrt(sd_s(l,j)/(length-1)) sdlog_s(l,j)=sqrt(sdlog_s(l,j)/(length-1)) do 8063 k=1,length a_theta(k)=atheta(l,k) 8063 continue call dsvrgn(length,a_theta,b_theta) a_l(l,j)=0.5*(b_theta(n_lower1)+b_theta(n_lower2)) a_u(l,j)=0.5*(b_theta(n_upper1)+b_theta(n_upper2)) est_exp(l,j)=exp(theta_sum(l,j)+0.5*sdlog_s(l,j)**2) 8062 continue 8061 continue acc=acc/(nmcmc*n_sub) write (9,870) nmcmc,nsmcmc,length,nk,ISEED write (9,883) acc, n1,n2,ISEED, TOL write (9,871) (exptmp(k),k=1,7) write (9,8710) (tmpexp(k),k=1,7) write (9,8711) estexp write (9,872) a_mulower write (9,873) a_muupper write (9,874) a_mumedian write (9,875) sd write (9,8751) sdlog write (9,8752) vinlog write (9,876) tmp do 100 l=1,n_sub write (9,884) l write (9,871) (theta_exp(l,j),j=1,n_par) write (9,8710) (thetaexp(l,j),j=1,n_par) write (9,8711) (est_exp(l,j),j=1,n_par) write (9,872) (a_l(l,j),j=1,n_par) write (9,873) (a_u(l,j),j=1,n_par) write (9,874) (a_me(l,j),j=1,n_par) write (9,875) (sd_s(l,j),j=1,n_par) write (9,8751) (sdlog_s(l,j),j=1,n_par) write (9,8752) (vinlog_s(l,j),j=1,n_par) write (9,876) (theta_sum(l,j),j=1,n_par) write (4,770) (est_exp(l,j),j=1,n_par) write (5,771) (theta_exp(l,j),j=1,n_par) 100 continue close (unit=9) close (unit=2) close (unit=3) close (unit=4) close (unit=5) close (unit=6) close (unit=7) 870 format (1X, i8, 1X,i6,1X,i8,1X,i3,1X,i7) 871 format('exptmp=', 5f10.3,2x,f10.3,1x,f11.7) 8710 format('tmpexp=', 5f10.3,2x,f10.3,1x,f11.7) 8711 format('estexp=', 5f10.3,2x,f10.3,1x,f11.7) 872 format('mulow=', 1x,5f10.3,2x,f10.3,1x,f11.7) 873 format('muupp=', 1x,5f10.3,2x,f10.3,1x,f11.7) 874 format('median=', 5f10.3,2x,f10.3,1x,f11.7) 875 format('sd=', 4x,5f10.3,2x,f10.3,1x,f11.7) 8751 format('sdlog=', 1x,5f10.3,2x,f10.3,1x,f11.7) 8752 format('vinlog=', 5f10.3,2x,f10.3,1x,f11.3) 876 format('tmp=', 3x,5f10.3,2x,f10.3,1x,f11.4) 770 format(5f10.3,2x,f10.3,1x,f11.7) 771 format(5f10.3,2x,f10.3,1x,f11.7) 772 format(5f10.3,2x,f10.3,1x,f11.7) 773 format(5f10.3,2x,f10.3,1x,f11.3) 774 format(5f10.3,2x,f10.3,1x,f11.6) 881 format (7(1X,F10.4)) 882 format(i6,1x,7(1X,F10.4)) 883 format(f7.4,1X,i6,1X,i6,1X,i7,1x, F8.5) 884 format('nosub=', i3) 990 format('1st=', f10.2,1x,f10.4) 991 format('2nd=', f10.2,1x,f10.4,f10.5,1x,f14.9) stop end C **************************************************** C * SUBROUTINES AND FUNCTIONS C **************************************************** C (1) To generate gamma samples subroutine gamma (x,alpha,beta) implicit double precision (a-h, o-z) external drngam nr=1 call drngam(nr,alpha,x) x=x*beta return end C (2) To generate Wishart samples subroutine wish3(v,omega,a_nu) implicit double precision (a-h, o-z) integer n,ldr,irank,nr,nra,nca,lda,nrb,ncb,ldb,nrc,ncc,ldc,nb parameter(n_par=7) dimension v(n_par,n_par),omega(n_par,n_par) dimension chol(n_par,n_par),awish(n_par,n_par),hwish(n_par,n_par) intrinsic sqrt external drnnor,drngam,dchfac,dmrrrr,dmxtxf,dmach nr=1 n=n_par lda=n_par ldr=n_par tol=100.0*dmach(4) call dchfac(n,omega,lda,tol,irank,chol,ldr) call drnnor(nr,awish(1,2)) call drnnor(nr,awish(1,3)) call drnnor(nr,awish(1,4)) call drnnor(nr,awish(1,5)) call drnnor(nr,awish(1,6)) call drnnor(nr,awish(1,7)) call drnnor(nr,awish(2,3)) call drnnor(nr,awish(2,4)) call drnnor(nr,awish(2,5)) call drnnor(nr,awish(2,6)) call drnnor(nr,awish(2,7)) call drnnor(nr,awish(3,4)) call drnnor(nr,awish(3,5)) call drnnor(nr,awish(3,6)) call drnnor(nr,awish(3,7)) call drnnor(nr,awish(4,5)) call drnnor(nr,awish(4,6)) call drnnor(nr,awish(4,7)) call drnnor(nr,awish(5,6)) call drnnor(nr,awish(5,7)) call drnnor(nr,awish(6,7)) shape=a_nu/2.0 call drngam(nr,shape,awish(1,1)) awish(1,1)=sqrt(2.0*awish(1,1)) shape=(a_nu-1.0)/2.0 call drngam(nr,shape,awish(2,2)) awish(2,2)=sqrt(2.0*awish(2,2)) shape=(a_nu-2.0)/2.0 call drngam(nr,shape,awish(3,3)) awish(3,3)=sqrt(2.0*awish(3,3)) shape=(a_nu-3.0)/2.0 call drngam(nr,shape,awish(4,4)) awish(4,4)=sqrt(2.0*awish(4,4)) shape=(a_nu-4.0)/2.0 call drngam(nr,shape,awish(5,5)) awish(5,5)=sqrt(2.0*awish(5,5)) shape=(a_nu-5.0)/2.0 call drngam(nr,shape,awish(6,6)) awish(6,6)=sqrt(2.0*awish(6,6)) shape=(a_nu-6.0)/2.0 call drngam(nr,shape,awish(7,7)) awish(7,7)=sqrt(2.0*awish(7,7)) awish(2,1)=0.0 awish(3,1)=0.0 awish(3,2)=0.0 awish(4,1)=0.0 awish(4,2)=0.0 awish(4,3)=0.0 awish(5,1)=0.0 awish(5,2)=0.0 awish(5,3)=0.0 awish(5,4)=0.0 awish(6,1)=0.0 awish(6,2)=0.0 awish(6,3)=0.0 awish(6,4)=0.0 awish(6,5)=0.0 awish(7,1)=0.0 awish(7,2)=0.0 awish(7,3)=0.0 awish(7,4)=0.0 awish(7,5)=0.0 awish(7,6)=0.0 nra=n_par nca=n_par lda=n_par nrb=n_par ncb=n_par ldb=n_par nrc=n_par ncc=n_par ldc=n_par call dmrrrr(nra,nca,awish,lda,nrb,ncb,chol,ldb,nrc,ncc,hwish,ldc) nb=n_par call dmxtxf(nra,nca,hwish,lda,nb,v,ldb) return end C (3) To generate multivariate normal samples subroutine mvnor3(x,a_mu,sigma) implicit double precision (a-h, o-z) integer i,irank,k,ldr,ldrsig,nr,lda parameter(n_par=7) dimension x(n_par),a_mu(n_par),sigma(n_par,n_par), + rsig(n_par,n_par) external dchfac,drnmvn,dmach k=n_par lda=n_par ldrsig=n_par tol=100.0*dmach(4) call dchfac(k,sigma,lda,tol,irank,rsig,ldrsig) nr=1 ldr=1 call drnmvn(nr,k,rsig,ldrsig,x,ldr) do 10 i=1,n_par x(i)=x(i)+a_mu(i) 10 continue return end C To solve differential equations C (4) SUBROUTINE FCN(N, T, Y, YPRIME) IMPLICIT DOUBLE PRECISION (A-H, O-Z) INTEGER N DIMENSION Y(N), YPRIME(N) COMMON /COM/PHI,C0,DELTA,ALAMBDA,RHO,ANUM,AK COMMON /COM1/Cmin_idv, C1_50idv, C2_50idv, AR50idv COMMON /COM2/Cmin_rtv, C1_50rtv, C2_50rtv, AR50rtv COMMON /COM3/ad_idv, ad_rtv C50idv=C50FUN(AR50idv,C1_50idv,C2_50idv,T) C50rtv=C50FUN(AR50rtv,C1_50rtv,C2_50rtv,T) adh1=ad_idv adh2=ad_rtv DNOM=Cmin_idv*adh1/C50idv+Cmin_rtv*adh2/C50rtv DEN=PHI+DNOM IF (T .EQ. 0.0) THEN GAMMA=0.0 ELSE GAMMA=DNOM/DEN END IF YPRIME(1)=ALAMBDA-RHO*Y(1)-(1-GAMMA)*AK*Y(1)*Y(3) YPRIME(2)=(1-GAMMA)*AK*Y(1)*Y(3)-DELTA*Y(2) YPRIME(3)=ANUM*DELTA*Y(2)-C0*Y(3) RETURN END C (5) Adherence SUBROUTINE ADHEREN(i, TEND, AD) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION AD(2) C **(1)** IF (i .EQ. 1) THEN IF (TEND .LE. 56.0) THEN ad_idv=0.79 ad_rtv=0.77 ELSE IF (TEND .LE. 149.0) THEN ad_idv=0.65 ad_rtv=0.54 ELSE IF (TEND .LE. 169.0) THEN ad_idv=0.83 ad_rtv=1.00 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(2)** ELSE IF (i .EQ. 2) THEN IF (TEND .LE. 13.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 26.0) THEN ad_idv=0.88 ad_rtv=0.77 ELSE IF (TEND .LE. 111.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 167.0) THEN ad_idv=0.98 ad_rtv=1.00 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(3)** ELSE IF (i .EQ. 3) THEN ad_idv=1.00 ad_rtv=1.00 C **(4)** ELSE IF (i .EQ. 4) THEN ad_idv=1.00 ad_rtv=1.00 C **(5)** ELSE IF (i .EQ. 5) THEN IF (TEND .LE. 34.0) THEN ad_idv=0.91 ad_rtv=0.93 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(6)** ELSE IF (i .EQ. 6) THEN IF (TEND .LE. 140.0) THEN ad_idv=0.96 ad_rtv=0.85 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(7)** ELSE IF (i .EQ. 7) THEN ad_idv=1.00 ad_rtv=1.00 C **(8)** ELSE IF (i .EQ. 8) THEN IF (TEND .LE. 64.0) THEN ad_idv=0.95 ad_rtv=0.95 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(9)** ELSE IF (i .EQ. 9) THEN IF (TEND .LE. 29.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 56.0) THEN ad_idv=1.00 ad_rtv=0.89 ELSE IF (TEND .LE. 86.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 114.0) THEN ad_idv=1.00 ad_rtv=0.89 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(10)** ELSE IF (i .EQ. 10) THEN ad_idv=1.00 ad_rtv=1.00 C **(11)** ELSE IF (i .EQ. 11) THEN IF (TEND .LE. 14.0) THEN ad_idv=0.97 ad_rtv=1.00 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(12)** ELSE IF (i .EQ. 12) THEN ad_idv=1.00 ad_rtv=1.00 C **(13)** ELSE IF (i .EQ. 13) THEN IF (TEND .LE. 112.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 141.0) THEN ad_idv=0.98 ad_rtv=0.98 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(14)** ELSE IF (i .EQ. 14) THEN IF (TEND .LE. 57.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 108.0) THEN ad_idv=1.00 ad_rtv=0.94 ELSE IF (TEND .LE. 128.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 163.0) THEN ad_idv=0.99 ad_rtv=0.99 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(15)** ELSE IF (i .EQ. 15) THEN IF (TEND .LE. 112.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 144.0) THEN ad_idv=0.98 ad_rtv=0.98 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(16)** ELSE IF (i .EQ. 16) THEN IF (TEND .LE. 28.0) THEN ad_idv=0.82 ad_rtv=0.84 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(17)** ELSE IF (i .EQ. 17) THEN IF (TEND .LE. 34.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 63.0) THEN ad_idv=0.98 ad_rtv=0.98 ELSE IF (TEND .LE. 91.0) THEN ad_idv=0.96 ad_rtv=0.91 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(18)** ELSE IF (i .EQ. 18) THEN ad_idv=1.00 ad_rtv=1.00 C **(19)** ELSE IF (i .EQ. 19) THEN IF (TEND .LE. 13.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 28.0) THEN ad_idv=0.33 ad_rtv=0.33 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(20)** ELSE IF (i .EQ. 20) THEN ad_idv=1.00 ad_rtv=1.00 C **(21)** ELSE IF (i .EQ. 21) THEN ad_idv=1.00 ad_rtv=1.00 c **(22)** ELSE IF (i .EQ. 22) THEN ad_idv=1.00 ad_rtv=1.00 C **(23)** ELSE IF (i .EQ. 23) THEN ad_idv=1.00 ad_rtv=1.00 C **(24)** ELSE IF (i .EQ. 24) THEN ad_idv=1.00 ad_rtv=1.00 C **(25)** ELSE IF (i .EQ. 25) THEN IF (TEND .LE. 13.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 56.0) THEN ad_idv=0.99 ad_rtv=0.99 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(26)** ELSE IF (i .EQ. 26) THEN IF (TEND .LE. 14.0) THEN ad_idv=0.93 ad_rtv=1.00 ELSE IF (TEND .LE. 28.0) THEN ad_idv=0.96 ad_rtv=0.96 ELSE IF (TEND .LE. 54.0) THEN ad_idv=1.00 ad_rtv=0.98 ELSE IF (TEND .LE. 82.0) THEN ad_idv=0.98 ad_rtv=0.99 ELSE IF (TEND .LE. 112.0) THEN ad_idv=0.83 ad_rtv=0.83 ELSE IF (TEND .LE. 146.0) THEN ad_idv=0.72 ad_rtv=0.68 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(27)** ELSE IF (i .EQ. 27) THEN IF (TEND .LE. 14.0) THEN ad_idv=0.34 ad_rtv=0.66 ELSE IF (TEND .LE. 28.0) THEN ad_idv=0.75 ad_rtv=0.77 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(28)** ELSE IF (i .EQ. 28) THEN IF (TEND .LE. 82.0) THEN ad_idv=0.87 ad_rtv=0.96 ELSE IF (TEND .LE. 112.0) THEN ad_idv=0.93 ad_rtv=0.95 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(29)** ELSE IF (i .EQ. 29) THEN IF (TEND .LE. 20.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 99.0) THEN ad_idv=0.97 ad_rtv=0.99 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(30)** ELSE IF (i .EQ. 30) THEN ad_idv=1.00 ad_rtv=1.00 C **(31)** ELSE IF (i .EQ. 31) THEN ad_idv=1.00 ad_rtv=1.00 C **(32)** ELSE IF (i .EQ. 32) THEN ad_idv=1.00 ad_rtv=1.00 C **(33)** ELSE IF (i .EQ. 33) THEN IF (TEND .LE. 112.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 140.0) THEN ad_idv=0.98 ad_rtv=0.96 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(34)** ELSE IF (i .EQ. 34) THEN IF (TEND .LE. 20.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 61.0) THEN ad_idv=0.96 ad_rtv=0.96 ELSE IF (TEND .LE. 90.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 119.0) THEN ad_idv=0.98 ad_rtv=0.98 ELSE IF (TEND .LE. 166.0) THEN ad_idv=0.98 ad_rtv=0.96 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(35)** ELSE IF (i .EQ. 35) THEN IF (TEND .LE. 26.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 53.0) THEN ad_idv=0.99 ad_rtv=1.00 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(36)** ELSE IF (i .EQ. 36) THEN ad_idv=1.00 ad_rtv=1.00 C **(37)** ELSE IF (i .EQ. 37) THEN ad_idv=1.00 ad_rtv=1.00 C **(38)** ELSE IF (i .EQ. 38) THEN IF (TEND .LE. 56.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 84.0) THEN ad_idv=0.93 ad_rtv=0.94 ELSE IF (TEND .LE. 111.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 140.0) THEN ad_idv=0.98 ad_rtv=0.98 ELSE ad_idv=1.00 ad_rtv=1.00 END IF C **(39)** ELSE IF (i .EQ. 39) THEN ad_idv=1.00 ad_rtv=1.00 C **(40)** ELSE IF (i .EQ. 40) THEN ad_idv=1.00 ad_rtv=1.00 C **(41)** ELSE IF (i .EQ. 41) THEN ad_idv=1.00 ad_rtv=1.00 C **(42)** ELSE IF (i .EQ. 42) THEN IF (TEND .LE. 15.0) THEN ad_idv=0.97 ad_rtv=0.90 ELSE IF (TEND .LE. 25.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 56.0) THEN ad_idv=0.97 ad_rtv=0.97 ELSE IF (TEND .LE. 84.0) THEN ad_idv=1.00 ad_rtv=1.00 ELSE IF (TEND .LE. 140.0) THEN ad_idv=1.00 ad_rtv=0.99 ELSE ad_idv=1.00 ad_rtv=1.00 END IF END IF AD(1)=ad_idv AD(2)=ad_rtv RETURN END C (6) Function IC50(t) FUNCTION C50FUN(D,C1,C2,T) IMPLICIT DOUBLE PRECISION (A-H, O-Z) IF (D .EQ. 0.0) THEN C50=C1 ELSE IF (T .LE. D) THEN C50=C1+T*(C2-C1)/D ELSE C50=C2 END IF C50FUN=C50 RETURN END