;; 1. Based on: 127
$PROBLEM    DIG TIME on HAZ step BP est & AGE & SEX
;;Digoxin PK-HR-CONVERSION DATA 
$INPUT      ID TRET SEX AGE WT CRCL TIME DV AMT RATE FLG
; TRET is treatment, 1=Digoxin, 0=placebo

; FLG is FLAG, 0=HR data, 2=Dose 1=Digoxin conc, 3=Conversion, 4=Censored
$DATA      dig_all1_simulated.csv IGNORE=@

$SUBROUTINE ADVAN6 TRANS1 TOL=5
$MODEL      COMP=CENT COMP=PERI COMP=EFF COMP=HAZ
$PK 
;---------PK (parameters hardcoded)---------------
       CRC  = CRCL
         IF(CRCL.EQ.-99) CRC = 70     ;missing in treatment group
         IF(CRCL.EQ.0)   CRC = 70     ;missing in placebo group
       FCRC = 1+.01985*(CRC-70)
       CL   = 9.87*FCRC*EXP(ETA(3))
       V1   = 27.8
       Q    = 71.8
       V2   = 444
       S1   = V1
       K12  = Q/V1
       K21  = Q/V2
       K10  = CL/V1
;--------- Effect compartment--------------
       K13  = K10*.001  ;negliable loss to effect compt
       K31  = THETA(5) ; t1/2 for ditribution to the effect compt (=ln2/keo)
       S3   = S1*K13/K31 ;preserves Cess=Cpss

;---------Heart rate (HR)------------------
       BASE = THETA(1)*EXP(ETA(2))
       TEFF = THETA(6)*EXP(ETA(1)); placebo, time effect
       SLOP = THETA(4) ; digoxin effect
;---------Sinus conversion-----------------
$DES 
 DADT(1) = K21*A(2) - (K12+K10) * A(1) ; 3-comp kinetic model
 DADT(2) =-K21*A(2) +  K12      * A(1)
 DADT(3) =-K31*A(3) +  K13      * A(1)

    HR   = BASE*(1+TEFF*T)*(1+SLOP*A(3)/S3)  ; Model for heart rate

BP=THETA(12)
 HZ = THETA(7)
 IF(SEX.EQ.0) HZ=THETA(8)
 IF(T.GT.BP.AND.SEX.EQ.1) HZ = THETA(9)
 IF(T.GT.BP.AND.SEX.EQ.0) HZ = THETA(10)

 DADT(4) = HZ *EXP(-AGE*THETA(11))                    ; Integrate hazard

$ERROR 
;-------------PK and HR-----------------; continous data
 CONC = A(1)/S1
 CE   = A(3)/S3

 IHR   = BASE*(1+TEFF*TIME)*(1+SLOP*CE) ; VARIABLE DEFINED IN $DES MUST BE REDEFINED WITH NEW NAME IN $ERROR
 IPRED = IHR       ; For FLG=0, i.e. HR
 W = THETA(3)*IPRED
  
 IF(FLG.EQ.1) IPRED = CONC  ; For FLG=0, i.e. digoxin conc
 IF(FLG.EQ.1) W = THETA(2)

 IRES = DV-IPRED
 IWRES= IRES/W

 IF(FLG.LT.2) THEN          ; Digoxin conc and HR cont variables
    F_FLAG= 0
    Y  = IPRED+EPS(1)*W
 ENDIF

;------------TIME TO EVENT--------------   ; categorical data  

  CHZ   = A(4)              ; Cumulative hazard
  SUR   = EXP(-CHZ)         ; Survival

BRP=THETA(12)
 HAZ = THETA(7)
 IF(SEX.EQ.0) HAZ=THETA(8)
 IF(TIME.GT.BRP.AND.SEX.EQ.1) HAZ = THETA(9)
 IF(TIME.GT.BRP.AND.SEX.EQ.0) HAZ = THETA(10)

IF(FLG.EQ.3)    IPRED = SUR*(HAZ*EXP(-AGE*THETA(11)))  ; Survival * hazard (S*h), conversion
IF(FLG.EQ.4)    IPRED = SUR       ; Censored, no conversion during study interval
IF(FLG.GT.2) THEN
                F_FLAG= 1         ; Coversion and censoring events are categorical 
                Y     = IPRED
ENDIF

$THETA  (0,117) FIX ; 1 BASE HR
$THETA  (0,1.42) FIX ; 2 W CONC
$THETA  (0,0.1) FIX ; 3 W HR
$THETA  -0.1030 FIX ; 4 SLOP
$THETA  (0,0.174) FIX ; 5 KEO
$THETA  -0.00175 FIX ; 6 TIME-EFFECT
$THETA  (0,0.145483) ; 7 BHAZ FOR CONVERSION DIG+Male
$THETA  (0,0.285785) ; 8 BHAZ FOR CONVERSION DIG+female
$THETA  (0,0.068825) ; 9 BHAZ FOR CONVERSION DIG+M IF TIME >6
$THETA  (0,0.114181) ; 10 BHAZ FOR CONVERSION DIG+f IF TIME >6
$THETA  0.0166985 ; 11 AGE
$THETA  (0,4.16725) ; 12 BP

$OMEGA  1.510000  FIX  ; 1 TIME ON BASE IIV
$OMEGA  0.031300  FIX  ; 2 BASE IIV
$OMEGA  0.970000  FIX  ;   3 CL IIV

$SIGMA  1.000000  FIX

$ESTIMATION MAX=0 LAPLACIAN MET=1 PRINT=1 MSFO=msfb

;$COVARIANCE

$TABLE      ID TIME IPRED IWRES FLG CHZ SUR ONEHEADER NOPRINT
            FILE=sdtab127
$TABLE      ID SEX AGE WT CRCL NOPRINT ONEHEADER FILE=cotab127
$TABLE      ID BASE TEFF SLOP K31 HAZ CHZ SUR NOPRINT ONEHEADER
            FILE=patab127
$TABLE      ID TRET AMT NOPRINT ONEHEADER FILE=mytab127