$PROBLEM    CMS+Colistin PK 
 ;$INPUT      ID TIME AMT RAT0 TAG CONC BQL DV ADDL II EVID CMT OCC TAD
 ;           GEND=DROP AGE=DROP WT IWT=DROP SCR CRCL SALB=DROP NAD=DROP
 ;           LACT=DROP OSAT=DROP FLBA=DROP APAC=DROP DIAB=DROP
 ;           SEPC=DROP ARF=DROP OEDM=DROP HCT=DROP HB=DROP STUD DOSE
 ;           COLA=DROP CPLB=DROP COAB=DROP DURU INF TYPE RATE
 ;$DATA      data.csv   
$INPUT ID TIME AMT BQL DV ADDL II EVID CMT OCC TAD SCRCL INF TYPE RATE

$DATA       Simulated_CMS_colistin_PK_data.csv
            IGNORE=(ID.EQ.#) 

$SUBROUTINE ADVAN5
$MODEL      COMP=CMS1 COMP=CMS2 COMP=COL1 COMP=INTM COMP=INTM2  
$PK  

OCC1=0
OCC4=0
OCC6=0
OCC42 = 0
IF(OCC.EQ.1) OCC1 = 1
IF(OCC.EQ.2.OR.OCC.EQ.4) OCC4 = 1
IF(OCC.EQ.6.OR.OCC.EQ.7.OR.OCC.EQ.8) OCC6=1
IF(OCC.EQ.5.OR.OCC.EQ.42) OCC42=1

IOVCL = OCC1*ETA(8)  + OCC4*ETA(9) + OCC6*ETA(10) + OCC42*ETA(11)
IOVFM = OCC1*ETA(12) + OCC4*ETA(13) + OCC6*ETA(14) + OCC42*ETA(15)
IOVV2 = OCC1*ETA(16) + OCC4*ETA(17) + OCC6*ETA(18) + OCC42*ETA(19)

CLCR = CRCL
IF(CRCL.GT.150) CLCR = 150  ; assume higher CRCL is not reflecting CLR
IF(CRCL.EQ.-99) CLCR = 80   ; some IDs have missing CRCL

CLR = THETA(15)*CLCR*60/1000*EXP(ETA(1)+IOVCL) ; renal CMS
CLH = THETA(1)*EXP(ETA(1)+IOVCL)  ; hydrolysis/non-renal CMS
V1  = THETA(2)*EXP(ETA(2))        ; V1 CMS
Q2  = THETA(14)                   ; Q2 CMS
Q   = THETA(4)*EXP(ETA(3))        ; Q CMS
V2  = THETA(5)*EXP(ETA(5)+IOVV2)  ; V2 CMS
CLCO = THETA(6)*EXP(ETA(1)*THETA(11)+IOVFM) ; CL colistin  
VCO = THETA(7)*EXP(IOVFM) ; V colistin
CL2 = CLH             ; hydrolysis peripheral 

SCALE=THETA(11)

K12 = Q/V1
K21 = Q/V2
K43 = CLH/V1
K10 = CLR/V1
K40 = CLR/V1
K31 = 0
K30 = CLCO/VCO
K14 = K43
K45 = Q2/V1
K54 = Q2/V2
K25 = CL2/V1

F1 = THETA(12)
F4 = 1 -F1
IF(TYPE.LE.2) F1 = THETA(13)*THETA(16) 
IF(TYPE.LE.2) F4 = THETA(13)*(1-THETA(16))
D1 = INF  ; infusion duration
D4 = INF

$ERROR  
IPR = 0.00001
IF(CMT.EQ.3) IPR = A(3)/VCO
IF(CMT.EQ.1) IPR= A(1)/V1 + A(4)/V1

W = SQRT(THETA(3)**2+THETA(9)**2/(IPR**2+0.00001))*EXP(ETA(4))
IF(CMT.EQ.1) W = SQRT(THETA(8)**2+THETA(10)**2/(IPR**2+0.00001))*EXP(ETA(7))

LLOQ = LOG(0.030)
IPRED = LOG(IPR)        ; Ln-transformed data/IPRED
DUM  = (LLOQ-IPRED)/W   ; Increase the lower the IPRED

;-- PHI(X) is the integral from .INF to X of N(0,1) .
DUM2 = PHI(DUM)	            ; Probability < LLOQ

IWRES=0
Y = 1

  F_FLAG=0
  IRES = DV-IPRED
  IWRES = IRES/W
  Y     = IPRED+EPS(1)*W

IF(BQL.EQ.1) THEN
  F_FLAG=1
  IWRES = 0
  Y=DUM2
ENDIF

STRA = 10+OCC
IF(CMT.EQ.3) STRA=30+OCC

$THETA  
 (0,5.69,100)  ; 1 CL
 (0,1.74,200)  ; 2 V1
 (0,0.0894,.5) ; 3 Prop Err Col
 (1,588,4000)  ; 4 Q
 (0,12.1,400)  ; 5 V2
 (0,4.69,300)  ; 6 CL colistin
 (0,76.6,3000) ; 7 V1 colistin
 (0,0.154,.7)  ; 8 Prop error CMS
 (0,0.0618,1)  ; 9 Add err Col
 (0,0.168,1)   ; 10 Add err CMS
 (0,4.02)      ; 11 Scale corr between CLCMS and CLCol
 1 FIX         ; 12 Total F (and F1) Study 3
 (0,0.69,1)    ; 13 Total F Study 1&2
 (0,7.18)      ; 14 Q2
 (0,0.528,10)  ; 15 CRCL-CL cov relationship 
 (0,0.89)      ; 16 F1 study 1&2 (prop of Tot F)

$OMEGA  0.032  ; 1 CL IIV and CL colistin
$OMEGA  0  FIX ; 2 V1
$OMEGA  0  FIX ; 3 Q
$OMEGA  0.212  ; 4 Res err col
$OMEGA  0  FIX ; 5 V2
$OMEGA  0  FIX ; 6 V col
$OMEGA  0  FIX ; 7
$OMEGA  BLOCK(1)
 0.151         ; 8-11 IOVCL
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1)
 0.181         ; 12-15 IOVFM
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1)
 0.094         ; 16-19 IOVV2
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$OMEGA  BLOCK(1) SAME
$SIGMA  1  FIX
$ESTIMATION MAXEVALS=9999 PRINT=2 METHOD=1 LAPLACE INTER NOABORT
$COVARIANCE UNCOND MATRIX=S
$TABLE      ID TIME TAD IPRED IWRES CWRES DV PRED NOAPPEND ONEHEADER NOPRINT
            FILE=sdtab1099_ddmore
$TABLE      ID TIME CLH CLR V1 Q V2 CLCO VCO SCALE ETA(1)
            ETA(3) ONEHEADER NOPRINT FILE=patab1099_ddmore
$TABLE      ID OCC TAG STRA CMT ONEHEADER NOPRINT FILE=catab1099_ddmore
$TABLE      ID WT CLCR ONEHEADER NOPRINT FILE=cotab1099_ddmore