Continuous-time Markov model for Ordinal Categorical Score Data
Hi All,
I was trying to implement the similar model structure with some
modification as this publication below.
CPT Pharmacometrics Syst Pharmacol. 2014 Oct 29;3:e143.
Simultaneous Exposure-Response Modeling of ACR20, ACR50, and ACR70
Improvement Scores in Rheumatoid Arthritis Patients Treated With
Certolizumab Pegol.
Lacroix BD1, Karlsson MO2, Friberg LE2.
My data consist of six scores 0,1,2,3,4,5.
In the draft model, I assumes all the transition rate constants are the
same for upward or downward direction except for the Score 4 to Score 5,
where I have another upward rate constant.
I have reasonable estimation of the three THETAs, since the data is the
vehicle group so the score is increasing.
1.37E+00 ; Kup ; worsen
8.05E-01 ; Kdown ; improve
2.11E-01 ; Kup2 ; worsen2
However, I ran into some error and questions when I use the code from the
publication for simulation to do predictive check.
The code was attached as followed.
I could have ran the simulation without termination but have few warnings,
and I would like to make sure it is not due to the code error.
Here are the error messages.
(WARNING 2) NM-TRAN INFERS THAT THE DATA ARE POPULATION.
(WARNING 102) $ERROR: Y IS NOT SET AT ICALL=2. VALUES OF PRED DISPLAYED BY
$TABLE OR $SCATTER MAY BE INCORRECT.
(WARNING 85) EPS VARIABLES REQUIRED WITH ANALYSIS OF POPULATION DATA.
PRESENT ABBREVIATED CODE COULD CAUSE AN ERROR IN A SUBSEQUENT RUN IF
CHECKDATA OR ONLYSIM IS NOT USED.
(WARNING 37) $ERROR: TO SIMULATE THE DV ITEM, Y SHOULD BE SET TO THE
SIMULATED VALUE.
(WARNING 52) $ERROR: TO SIMULATE THE DV ITEM, DV SHOULD BE SET TO THE
SIMULATED VALUE.
Also, I would like to understand better for this part of code. How does
NONMEM scale the amount in each compartment into probability which less
than 1.
;scale to 1
P1 = AA1
P2 = P1+AA2
P3 = P2+AA3
P4 = P3+AA4
P5 = P4+AA5
P6 = 1
Thank you very much,
Any suggestion is much appreciated,
Wu
------------------Code of Control file ------------------------------------
$PROBLEM PROJECT Markovsim
$INPUT C ID STD TIME DV AMT CMT TYPE EVID
$SUBROUTINE ADVAN6 TOL=5
$DATA Markovsim.csv IGNORE=C
$MODEL
COMP=(PR0) ; score=0
COMP=(PR1) ; score=1
COMP=(PR2) ; score=2
COMP=(PR3) ; score=3
COMP=(PR4) ; score=4
COMP=(PR5) ; score=5
$PK
REP=IREP
IF(NEWIND.NE.2) THEN
PDV=0
F1=1
F2=0
F3=0
F4=0
F5=0
F6=0
ENDIF
PDV0=0
PDV1=0
PDV2=0
PDV3=0
PDV4=0
PDV5=0
;TURN on/off the correct CPT according to simulated DV2
IF(PDV.EQ.0) THEN
PDV0=1
F1=1
F2=0
F3=0
F4=0
F5=0
F6=0
ENDIF
IF(PDV.EQ.1) THEN
PDV1=1
F1=0
F2=1
F3=0
F4=0
F5=0
F6=0
ENDIF
IF(PDV.EQ.2) THEN
PDV2=1
F1=0
F2=0
F3=1
F4=0
F5=0
F6=0
ENDIF
IF(PDV.EQ.3) THEN
PDV3=1
F1=0
F2=0
F3=0
F4=1
F5=0
F6=0
ENDIF
IF(PDV.EQ.4) THEN
PDV4=1
F1=0
F2=0
F3=0
F4=0
F5=1
F6=0
ENDIF
IF(PDV.EQ.5) THEN
PDV5=1
F1=0
F2=0
F3=0
F4=0
F5=0
F6=1
ENDIF
;MODEL parameters--------------
;--Rate constants---
TVKup = THETA(1)
Kup = TVKup*EXP(ETA(1))
TVKdown = THETA(2)
Kdown = TVKdown*EXP(ETA(2))
TVKup2 = THETA(3)
Kup2 = TVKup2*EXP(ETA(3))
$DES
KTup = Kup
KTdown = Kdown
KTup2 = Kup2
DADT(1) = KTdown*A(2) - KTup*A(1)
DADT(2) = (KTup*A(1)+KTdown*A(3)) - (KTup*A(2)+KTdown*A(2))
DADT(3) = (KTup*A(2)+KTdown*A(4)) - (KTup*A(3)+KTdown*A(3))
DADT(4) = (KTup*A(3)+KTdown*A(5)) - (KTup*A(4)+KTdown*A(4))
DADT(5) = (KTup*A(4)+KTdown*A(6)) - (KTup2*A(5)+KTdown*A(5))
DADT(6) = KTup2*A(5) - KTdown*A(6)
$ERROR
; OUTPUT THE AMOUNT IN EACH COMPARTMENT
AA1 = A(1)
AA2 = A(2)
AA3 = A(3)
AA4 = A(4)
AA5 = A(5)
AA6 = A(6)
; check the sum
SUM = A(1)+A(2)+A(3)+A(4)+A(5)+A(6)
;scale to 1
P1 = AA1
P2 = P1+AA2
P3 = P2+AA3
P4 = P3+AA4
P5 = P4+AA5
P6 = 1
;simulation block --------------
IF(ICALL.EQ.4) THEN
CALL RANDOM(2,R)
RAN=R
;initialized internal variables----
IF(NEWIND.NE.2) THEN
PDV=0
;PREV=0
DV2=0
DVA=0
;PDVA=0
PS=0
ENDIF
;simulate socres------
IF(TYPE.EQ.1)THEN
IF(R.LE.P1)DV2=0
IF(R.GT.P1.AND.R.LE.P2)DV2=1
IF(R.GT.P2.AND.R.LE.P3)DV2=2
IF(R.GT.P3.AND.R.LE.P4)DV2=3
IF(R.GT.P4.AND.R.LE.P5)DV2=4
IF(R.GT.P5.AND.R.LE.P6)DV2=5
ENDIF
ENDIF ; END of simulation block
; track the previous Score
IF(TYPE.EQ.1)THEN
DVA=DV2
PS=PDV
PDV=DVA
ENDIF
$THETA
1.37E+00 ; Kup ; worsen
8.05E-01 ; Kdown ; improve
2.11E-01 ; Kup2 ; worsen2
$OMEGA BLOCK(3)
4.27E-02
3.61E-02 ,3.05E-02
3.22E-02 ,2.73E-02 ,2.44E-02
$SIGMA 1 FIX
$SIMULATION (0238721)(14201 UNIFORM) ONLYSIM SUBPROBS=100
$TABLE REP ID TIME TYPE DV2 PS NOAPPEND NOHEADER
NOPRINT FILE=Markovsim.txt