$PROB Steve's hierarchical interaction model
; Using TCI
;1 react to shake and shout 2 no react to shake and shout; 3 no react to int.

$DATA binary.csv

$INPUT  ID OP DROP=IREM HY DROP=IPRO INT=DV

$PRED 
; Note: Drug OP = remifentanil (opioid), Drug HY = propofol (hypnotic)

;define parameters
OC50=THETA(1)*EXP(ETA(1)) ; Opioid C50
HC50=THETA(2)*EXP(ETA(2)) ; Hypnotic C50
OG=THETA(3)*EXP(ETA(3))   ; Opioid Gamma
HG=THETA(4)*EXP(ETA(4))   ; Hypnotic Gamma
PN1=THETA(5)*EXP(ETA(5))  ; Pain of Shake/Shout
PN2=THETA(6)*EXP(ETA(6))  ; Pain of Laryngoscopy

OPG=OP**OG
HYG=HY**HG

;Compute Opioid Out for each stimulus
OO1 = PN1*(1 - OPG/((OC50*PN1)**OG + OPG))
OO2 = PN2*(1 - OPG/((OC50*PN2)**OG + OPG))

;Compute probabilities of each stimulus
P1=HYG/((HC50*OO1)**HG + HYG)
P2=HYG/((HC50*OO2)**HG + HYG)

Y1=(1-P1)*(1-P2)
Y2=P1*(1-P2)
Y3=P1*P2

;Combine probababilities
IF (INT.EQ.1) THEN
Y=Y1
ENDIF

IF (INT.EQ.2) THEN
Y=Y2
ENDIF

IF (INT.EQ.3) THEN
Y=Y3
ENDIF

$THETA 
(1,1.07,30)    ; ec50 remi
(1,8.9,20)     ; ec50 prop 
(.1, .86, 5)   ; Opioid gamma
(.1, 5, 7)     ; Hypnotic gamma
(.1,.53,10)    ; Pain of shake/shout
(.1,.923, 10)  ; Pain of intubation
$OMEGA  .07 0 FIX .05 .8 0 FIX 0 FIX

$ESTIMATION SIG=4 MAX=2000 PRINT=5 NOABORT METHOD=1 LAPLACE LIKELIHOOD
$COVARIANCE