Model for Simulation

;;;Model for simulation

$SIZES NO=1000 LIM6=1000 ;Re-set predefault value

$PROB control stream ;Model for simulation

$DATA screening_V4_sde5.csv IGNORE=@

$INPUT ID DROP TIME DV PCNT DROP DROPDROPDROPDROPDROPDROP

PDV DROP DROPDROPDROPDROPDROPMDUR EVID CMT SDE AMT=DROP

$PRED

ALPHA=THETA(2)

LAMB =(THETA(1)+ALPHA*TIME)*EXP(ETA(1))

;;Simulation code for Poisson model

IF(ICALL.EQ.4.AND.EVID.EQ.0) THEN

T=0

N=0

DO WHILE (T.LE.1)

CALL RANDOM(2,R)

T=T-LOG(1-R)/LAMB

IF(T.LT.1) N=N+1

END DO

DV=N

END IF

IF(ICALL.EQ.4.AND.EVID.NE.0) DV=0

; simulation end

$THETA (0,50) ; Lambda

$THETA (-0.1) ; ALPHA

$OMEGA

0.09 ;

;$ESTIMATION MAXEVAL=9999 METHOD=COND LAPLACE -2LL MSFO=msf3100

;$COV PRINT=E

$SIM (123091) (123091 UNIFORM) ONLYSIM NOPRED

; Estimate using constant model

$SIZES NO=1000 LIM6=1000 ;Re-set pre-default value

$PROB control stream ;Constant model

$DATA screening_V4_sde5.csv IGNORE=@ IGNORE(EVID.GT.0)

$INPUT ID DROP TIME DV PCNT DROP DROPDROPDROPDROPDROPDROP

PDV DROP DROPDROPDROPDROPDROPMDUR EVID CMT SDE AMT=DROP

$PRED

LAMB=(THETA(1)+TIME*THETA(2))*EXP(ETA(1))

LDVFAC=0

IF (DV.GT.0) LDVFAC=DV*LOG(DV)- DV+LOG(DV*(1+4*DV*(1+2*DV)))/6+LOG(3.1415)/2

LYY=-LAMB+DV*LOG(LAMB)-LDVFAC

Y = -2*LYY

$THETA (0,50) ; Lambda pdv == 1

$THETA (-0.05) ; SLOPE

$OMEGA

0.09 ;

$ESTIMATION MAXEVAL=9999 METHOD=COND LAPLACE -2LL MSFO=msf3101

$COV PRINT=E

; Estimate using Slope model

$SIZES NO=1000 LIM6=1000 ;Re-set pre-default value

$PROB control stream ;Slope model

$DATA screening_V4_sde5.csv IGNORE=@ IGNORE(EVID.GT.0)

$INPUT ID DROP TIME DV PCNT DROP DROPDROPDROPDROPDROPDROP

PDV DROP DROPDROPDROPDROPDROPMDUR EVID CMT SDE AMT=DROP

$PRED

LAMB=(THETA(1)+TIME*THETA(2))*EXP(ETA(1))

LDVFAC=0

IF (DV.GT.0) LDVFAC=DV*LOG(DV)- DV+LOG(DV*(1+4*DV*(1+2*DV)))/6+LOG(3.1415)/2

LYY=-LAMB+DV*LOG(LAMB)-LDVFAC

Y = -2*LYY

$THETA (0,50) ; Lambda pdv == 1

$THETA (0) FIX ; SLOPE

$OMEGA

0.09 ;

$ESTIMATION MAXEVAL=9999 METHOD=COND LAPLACE -2LL MSFO=msf3102

$COV PRINT=E

; Estimate using dIOV model

$SIZES NO=1000 LIM6=1000 LVR=40 LVR2=40 ;Re-set pre-default value

$PROB control stream ;dIOV model

$DATA screening_V4_sde5.csv IGNORE=@ IGNORE(EVID.GT.0)

$INPUT ID DROP TIME DV PCNT DROP DROPDROPDROPDROPDROPDROP

PDV DROP DROPDROPDROPDROPDROPMDUR EVID CMT SDE AMT=DROP

$PRED

OCCL1=THETA(2) ;occasion length

GAM1=THETA(3) ;shape factor

SW1=OCCL1/2. ;width of surge

OCC1=ETA(2)/(((OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1) ;1st occasion

OCC2=ETA(3)/(((2*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1) ;2nd occasion

OCC3=ETA(4)/(((3*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1) ;.

OCC4=ETA(5)/(((4*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1) ;.

OCC5=ETA(6)/(((5*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1) ;.

OCC6=ETA(7)/(((6*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC7=ETA(8)/(((7*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC8=ETA(9)/(((8*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC9=ETA(10)/(((9*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC10=ETA(11)/(((10*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC11=ETA(12)/(((11*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC12=ETA(13)/(((12*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC13=ETA(14)/(((13*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC14=ETA(15)/(((14*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC15=ETA(16)/(((15*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC16=ETA(17)/(((16*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC17=ETA(18)/(((17*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC18=ETA(19)/(((18*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC19=ETA(20)/(((19*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC20=ETA(21)/(((20*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC21=ETA(22)/(((21*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC22=ETA(23)/(((22*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC23=ETA(24)/(((23*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC24=ETA(25)/(((24*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC25=ETA(26)/(((25*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1)

OCC26=ETA(27)/(((26*OCCL1-OCCL1/2.-TIME)**2/SW1**2)**GAM1+1);26th occasion

IOV=OCC1+OCC2+OCC3+OCC4+OCC5+OCC6+OCC7+OCC8+OCC9+OCC10+OCC11+OCC12

IOV=IOV+OCC13+OCC14+OCC15+OCC16+OCC17+OCC18+OCC19+OCC20+OCC21

IOV=IOV+OCC22+OCC23+OCC24+OCC25+OCC26

LAMB=(THETA(1))*EXP(ETA(1)+IOV)

LDVFAC=0

IF (DV.GT.0) LDVFAC=DV*LOG(DV)- DV+LOG(DV*(1+4*DV*(1+2*DV)))/6+LOG(3.1415)/2

LYY=-LAMB+DV*LOG(LAMB)-LDVFAC

Y = -2*LYY

$THETA (0,50) ; Lambda pdv == 1

$THETA (0,10) ; OCCL

$THETA (0,0.3) ; GAMMA

$OMEGA

0.09 ;

;; All dIOV eta share the same distribution

$OMEGA BLOCK(1) 0.1 ; IOV1

$OMEGA BLOCK(1) SAME ;

$ESTIMATION MAXEVAL=9999 METHOD=COND LAPLACE -2LL MSFO=msf3103

$COV PRINT=E

; Estimate using SDEs model

$SIZES NO=1000 LIM6=1000 ;Re-set pre-default value

$PROB control stream ;Estimate with SDE

$DATA screening_V4_sde5.csv IGNORE=@

$INPUT ID DROP TIME DV PCNT DROP DROPDROPDROPDROPDROPDROP

PDV DROP DROPDROPDROPDROPDROPMDUR EVID CMT SDE AMT=DROP

$MODEL

COMP(1) ;compartment for system

COMP(SDE) ;compartment for SDEs component

$SUBROUTINE ADVAN6 TOL=6

$THETA

(0, 0.5) ;LAMBDA 1

(0.01) ; SGW1

(1) FIX ; RV

$OMEGA

0.1

$PK

LAMB=THETA(1)*EXP(ETA(1))

SGW1=THETA(2)

RV=THETA(3)

; Initialize state and state covariance equations

IF(NEWIND.NE.2) THEN

AHT1=0

PHT1=0

A1=0

A2=0

ENDIF

;Store one-step predictions for EKF update

IF(EVID.NE.3) THEN

A1=A(1)

A2=A(2)

ELSE

A1=A1

A2=A2

ENDIF

;Store observations for EKF update

IF(EVID.EQ.0) OBS= DV

;EKF state update equations

IF(EVID.GT.2.AND.SDE.EQ.2) THEN

CMAT=1

RVAR=CMAT*A2*CMAT+THETA(3)**2

K1=A2*CMAT/(RVAR)

AHT1=A1+K1*(OBS-(A1+LAMB))

PHT1=A2-K1*RVAR*K1

ENDIF

;Update states

IF(A_0FLG.EQ.1) THEN

A_0(1)=AHT1

A_0(2)=PHT1

ENDIF

$DES

;Statepredicton eq.

DADT(1)=0

;State cov.

DADT(2)=SGW1*SGW1

$ERROR

LAMB1=LAMB+A1

LDVFAC=0

IF (DV.GT.0) LDVFAC=DV*LOG(DV)-DV+LOG(DV*(1+4*DV*(1+2*DV)))/6+LOG(3.1415)/2

LYY=-LAMB1+DV*LOG(LAMB1)-LDVFAC

Y = -2*LYY

$ESTIMATION MAXEVAL=9999 METHOD=COND LAPLACE -2LL MSFO=msf3104

$COV PRINT=E