$PROB Hidden Markov Process

$PROB Hidden Markov Process

; 0: sleep 1: awake

;##

$INPUT ID DOSE TIME TAD WAKE DAY DV PDV PTAD MDV CMT AMT RATE

;PDV is previous DV

;PTAD is previous TAD

$DATA NM.PD.csv IGNORE=#

$SUB ADVAN=6 TOL=3

$MODEL

COMP = (DOSE,DEFDOSE)

COMP = (CENTRAL)

COMP = (PERI)

COMP = (BIOPHASE)

$PK

;

CL = ;clearance

V2 = ;central volume

Q = ;intercompartmental clearance

V3 = ;peripheral volume

KA = ;absorption rate constant

F1 = 1

K20 = CL/V2

K23 = Q/V2

K32 = Q/V3

S2 = V2

KE0 = THETA(10)

$DES

DADT(1) = -KA*A(1)

DADT(2) = KA*A(1) - K20*A(2) - A(2)*K23 + K32*A(3)

DADT(3) = K23*A(2) - K32*A(3)

DADT(4) = KE0*(A(2)/V2-A(4))

$ERROR

CP = A(2)/V2

CE = A(4)

IF (DOSE.EQ.0) CP = 0

IF (DOSE.EQ.0) CE = 0

;------

;placebo max effect

PCBU = THETA(1) ; max PCB effect on acquisition of sleep

PCBV = THETA(2) ; max PCB effect on clearing sleep

;------

;Onset of PCB

KON = THETA(3)

;------

; Baseline acquisition of sleep

U0 = THETA(4)

; Baseline clearing sleep

V0 = THETA(5)

;------

;DRUG effect

EMAX_U = THETA(6)*EXP(ETA(2)) ;

EMAX_V = THETA(7)*EXP(ETA(2)) ;

EC50_U = THETA(8)*EXP(ETA(3)) ;

EC50_V = THETA(9)*EXP(ETA(3))

DRG_U = EMAX_U -EMAX_U/(1 + (CE/EC50_U)*EXP(CE-EC50_U))

DRG_V = EMAX_V -EMAX_V/(1 + (CP/EC50_V)*EXP(CP-EC50_V))

;------

; Random effect

ET = ETA(1)

;------

; Intensity of acquisition of sleep

u = EXP(U0 + PCBU*EXP(-KON*TAD) + DRG_U + ET)

; Intensity of clearing sleep

v = EXP(V0 + PCBV*EXP(-KON*TAD) + DRG_V + ET)

;------

; Probability of jumping from sleep to awake (P01)(sleep clearance)

P01=(v/(u+v))*(1-EXP(-(u+v)*(TAD-PTAD)))

; probability of acquisition odf sleep (P10)

P10=(u/(u+v))*(1-EXP(-(u+v)*(TAD-PTAD)))

IF (TAD.EQ.PTAD) P01 = 1

IF (TAD.EQ.PTAD) P10 = 0

;------

; Hidden state: PDV and DV are observed outcomes from states

; TPDV and TDV are the true unobserved states

P0 = 0;

P1 = 1

DELT = EXP(-ABS(WAKE-2.5))

IF (WAKE.GT.2.5.AND.WAKE.LT.4) P0 = EXP(THETA(11)*DELT +ETA(4))/(1+ EXP(THETA(11)*DELT+ETA(4))) ; P0 = P(DV =1| TDV = 0) (misclassification) 1-P0= P(DV=0|TDV=0); expect 1-P0>P0

IF (WAKE.LE.2.5.AND.WAKE.GT.1) P1 = EXP(THETA(12)*DELT+ETA(4))/(1+EXP(THETA(12)*DELT+ETA(4))) ; P1 = P(DV =1| TDV = 1) 1-P1= P(DV=0|TDV=1) (misclassification); expect P1>1-P1

;------

TP01=(1-P0)*P1*P01 + (1-P1)*P1*(1-P10) + (1-P0)*P0*(1-P01) + (1-P1)*P0*P10

TP10=P1*(1-P0)*P10 + P0*(1-P0)*(1-P01) + P1*(1-P1)*(1-P10) + P0*(1-P1)*P01

;Define likelihood - Prob(Y=m)

IF (PDV.EQ.0.AND.DV.EQ.1) Y=(1-P0)*P1*P01 + (1-P1)*P1*(1-P10) + (1-P0)*P0*(1-P01) + (1-P1)*P0*P10

IF (PDV.EQ.0.AND.DV.EQ.0) Y=(1-P0)*(1-P0)*(1-P01) + (1-P1)*(1-P0)*P10 + (1-P0)*(1-P1)*P01 + (1-P1)*(1-P1)*(1-P10)

IF (PDV.EQ.1.AND.DV.EQ.0) Y=P1*(1-P0)*P10 + P0*(1-P0)*(1-P01) + P1*(1-P1)*(1-P10) + P0*(1-P1)*P01

IF (PDV.EQ.1.AND.DV.EQ.1) Y=P1*P1*(1-P10) + P0*P1*P01 + P1*P0*P10 + P0*P0*(1-P01)

;------

$THETA

-8.47 ; 1 PCBU

5.16 ;2 PCBV

0.0452 ;3 KON

-2.95 ;4 U0

-4.17 ;5 V0

0.265 ;6 Emax U

0.442 ;7 Emax V

(0, 17.8) ;8 EC50 U

(0, 29.9) ;9 EC50 V

(0,0.0235) ; 10 KE0

-0.1 ;11 P0

0.1 ;12 P1

$OMEGA

0 FIX

0 FIX

0 FIX

0.01

$ESTIMATION MAX=9999 PRINT=10 NOABORT METHOD=1 LAPLACE LIKE

$COV MATRIX=R COMP

$TABLE ID TIME DOSE WAKE TAD PTAD DV PDV P01 P10 TP01 TP10 P0 P1 CP CE KE0 ETA(4)

NOAPPEND NOPRINT ONEHEADER

FILE=HMM.r.par