I hope I can get some help.
I know very well NONMEM and Phoenix NLME.
I am looking for an example that shows how to simulate with IQRtools multiple categorical response model.
The MONOLIX model is as follows:
;monolix code
[LONGITUDINAL]
input = { btongue4,btongue43,btongue432,btongue4321, bupper4,bupper43,bupper432,bupper4321, blower4,blower43,blower432,blower4321,
bjaw4,bjaw43,bjaw432,bjaw4321,
blips4,blips43,blips432,blips4321,
bmuscles4,bmuscles43,bmuscles432,bmuscles4321,
bneck4,bneck43,bneck432,bneck4321,
ai,kdecay,di0
}
DEFINITION:
tongue = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(tongue>=1)) = di0+ai*(1exp(kdecay*t))+btongue4+btongue43+btongue432+btongue4321
logit(P(tongue >=2)) = di0+ai*(1exp(kdecay*t))+btongue4+btongue43+btongue432
logit(P(tongue >=3)) = di0+ai*(1exp(kdecay*t))+btongue4+btongue43
logit(P(tongue >=4)) = di0+ai*(1exp(kdecay*t))+btongue4
}
upper = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(upper>=1)) = di0+ai*(1exp(kdecay*t))+bupper4+bupper43+bupper432+bupper4321
logit(P(upper>=2)) = di0+ai*(1exp(kdecay*t))+bupper4+bupper43+bupper432
logit(P(upper>=3)) = di0+ai*(1exp(kdecay*t))+bupper4+bupper43
logit(P(upper>=4)) = di0+ai*(1exp(kdecay*t))+bupper4
}
lower = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(lower>=1)) = di0+ai*(1exp(kdecay*t))+blower4+blower43+blower432+blower4321
logit(P(lower>=2)) = di0+ai*(1exp(kdecay*t))+blower4+blower43+blower432
logit(P(lower>=3)) = di0+ai*(1exp(kdecay*t))+blower4+blower43
logit(P(lower>=4)) = di0+ai*(1exp(kdecay*t))+blower4
}
jaw = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(jaw>=1)) = di0+ai*(1exp(kdecay*t))+bjaw4+bjaw43+bjaw432+bjaw4321
logit(P(jaw>=2)) = di0+ai*(1exp(kdecay*t))+bjaw4+bjaw43+bjaw432
logit(P(jaw>=3)) = di0+ai*(1exp(kdecay*t))+bjaw4+bjaw43
logit(P(jaw>=4)) = di0+ai*(1exp(kdecay*t))+bjaw4
}
lips = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(lips>=1)) = di0+ai*(1exp(kdecay*t))+blips4+blips43+blips432+blips4321
logit(P(lips>=2)) = di0+ai*(1exp(kdecay*t))+blips4+blips43+blips432
logit(P(lips>=3)) = di0+ai*(1exp(kdecay*t))+blips4+blips43
logit(P(lips>=4)) = di0+ai*(1exp(kdecay*t))+blips4
}
muscles = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(muscles>=1)) = di0+ai*(1exp(kdecay*t))+bmuscles4+bmuscles43+bmuscles432+bmuscles4321
logit(P(muscles>=2)) = di0+ai*(1exp(kdecay*t))+bmuscles4+bmuscles43+bmuscles432
logit(P(muscles>=3)) = di0+ai*(1exp(kdecay*t))+bmuscles4+bmuscles43
logit(P(muscles>=4)) = di0+ai*(1exp(kdecay*t))+bmuscles4
}
neck = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(neck>=1)) = di0+ai*(1exp(kdecay*t))+bneck4+bneck43+bneck432+bneck4321
logit(P(neck>=2)) = di0+ai*(1exp(kdecay*t))+bneck4+bneck43+bneck432
logit(P(neck>=3)) = di0+ai*(1exp(kdecay*t))+bneck4+bneck43
logit(P(neck>=4)) = di0+ai*(1exp(kdecay*t))+bneck4
}
OUTPUT:
output={jaw,lips,lower,muscles,neck,tongue,upper}
7 responses, each with 5 levels.
If I can get a simple example that explains me step by step how to proceed with this kind of simulation, I would really appreciate.
first of all, is it implemented? I did not see any examples of multi level, multi response categorical responses.
If it can be done, then more questions will come up about how to define the fixed effects, the random effects, the probability distributions etc..
Thanks a lot in advance
Serge Guzy
multiple categorical response simulation
Moderators: daniel.kaschek, intiquan

 Posts: 2
 Joined: Thu Jun 20, 2019 10:20 am

 Site Admin
 Posts: 7
 Joined: Mon Sep 03, 2018 11:50 am
 Location: Basel, Switzerland
 Contact:
Re: multiple categorical response simulation
Hi!
great question!
Yes you can do it. I will not provide examples for all of the code you pasted. But I will provide a simple runnable example for one  although simplified to one level ... as you are familiar with ordinal categorical models you will immediately see how to expand to multiple levels.
R code using IQR Tools comes first. Below is also the code for the "example.txt" model that you would need to copy into a file "example.txt".
Cheers,
Henning
great question!
Yes you can do it. I will not provide examples for all of the code you pasted. But I will provide a simple runnable example for one  although simplified to one level ... as you are familiar with ordinal categorical models you will immediately see how to expand to multiple levels.
R code using IQR Tools comes first. Below is also the code for the "example.txt" model that you would need to copy into a file "example.txt".
Cheers,
Henning
Code: Select all
library(IQRtools)
# Load model
model < IQRmodel("example.txt")
# 
# Single simulation
# 
# Simulate model
simres < sim_IQRmodel(model,simtime = 1000)
# Show results graphically
plot(simres)
# Get probability at defined endpoint time
p24 < simres$p_u_gt_1[simres$TIME==24]
# Now you have got the probability and thus you can sample
runif(1)<p24
# 
# Population simulation
# 
# Normally if you have fitted the model in NONMEM or MONOLIX you
# would get a list of individual parameters and IQR Tools would
# allow to directly import them as a data.frame.
# Here in this constructed example I will just create such a data.frame manually
# assuming 50% CV variability on all parameters
Nsubjects <20
indiv < data.frame(
kdecay = 0.01*exp(rnorm(Nsubjects,0,0.5)),
di0 = 0.1*exp(rnorm(Nsubjects,0,0.5)),
ai = 5*exp(rnorm(Nsubjects,0,0.5))
)
# Lets simulate the model for each of the parameters  again assuming we are
# interested in the outcome after 24 time units  here we first ccalculate the probabilities
probpopsim < sapply(1:Nsubjects, function (k) {
simres < sim_IQRmodel(model,simtime = 24,parameters = indiv[k,])
p24 < simres$p_u_gt_1[nrow(simres)]
})
# Next step  sampling from the probabilities  you can do that as often as you like
NREP < 10
do.call(rbind,lapply(1:NREP, function (krep) {
probpopsim > runif(length(probpopsim))
}))
# Otherwise the sky is the limit ... you can take into account uncertainty by just
# sampling with IQRtools from a MONOLIX or NONMEM fit.
Code: Select all
********** MODEL NAME
Provide Model Name
********** MODEL NOTES
Provide Model Information
********** MODEL STATES
d/dt(X) = kdecay*X
X(0) = 1
********** MODEL PARAMETERS
kdecay = 0.01
di0 = 0.1
ai = 5
********** MODEL VARIABLES
# Calculate the logit of your probability
# the RHS elements are parameters or states, etc.
# Only thing you need to change is t=>time
logit_p_u_gt_1 = di0+ai*(1X)
# Calculate your probability
p_u_gt_1 = exp(logit_p_u_gt_1)/(1+exp(logit_p_u_gt_1))
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

 Site Admin
 Posts: 7
 Joined: Mon Sep 03, 2018 11:50 am
 Location: Basel, Switzerland
 Contact:
Re: multiple categorical response simulation
Realized the example for a simple 0/1 thing might not be very convincing and depending on prior knowledge challenging to extrapolate to multiple levels. So here is an example a little bit extended to 5 categories (0,1,2,3,4) but only for "tongue" as for the others it really is copy and paste. Skipping also the population simulation step  it is shown in the previous example.
You could claim that for this simple example an underlying simulator is not really needed ... after all it just can be written down easily and quickly as there is no integration or anything else remotely complex present in the structural model. But if there is a more complex dynamic model driving the probabilities ... then you would like to implement it this way. Or think about  you have estimated a popPK model in NONMEM and have then built the ordered logistic model in MONOLIX and subsequently you would like to simulate the whole thing based on NONMEM and MONOLIX results. IQR Tools will do it. But we do not provide an interface to Phoenix. Again you could argue (and many do) "why using NONMEM for one model and MONOLIX for the other"? Well, it just happens that depending on the task at hand one tool as advantages and disadvantages for a specific task and vice versa.
Not sure if Phoenix provides an interface to both NONMEM and MONOLIX? Might have a look into it to see what it can do.
The R code
The example.txt file contents
You could claim that for this simple example an underlying simulator is not really needed ... after all it just can be written down easily and quickly as there is no integration or anything else remotely complex present in the structural model. But if there is a more complex dynamic model driving the probabilities ... then you would like to implement it this way. Or think about  you have estimated a popPK model in NONMEM and have then built the ordered logistic model in MONOLIX and subsequently you would like to simulate the whole thing based on NONMEM and MONOLIX results. IQR Tools will do it. But we do not provide an interface to Phoenix. Again you could argue (and many do) "why using NONMEM for one model and MONOLIX for the other"? Well, it just happens that depending on the task at hand one tool as advantages and disadvantages for a specific task and vice versa.
Not sure if Phoenix provides an interface to both NONMEM and MONOLIX? Might have a look into it to see what it can do.
The R code
Code: Select all
library(IQRtools)
model < IQRmodel("example.txt")
res < sim_IQRmodel(model,simtime = 1000)
plot(res)
# Simulate for defined time vector to get the probabilities
simtime < c(0,1,2,3,4,5,6)*168
res < sim_IQRmodel(model,simtime = simtime)
# Get cumulative probs per time point
propsT < cbind(res$P_tongue_0,res$P_tongue_1,res$P_tongue_2,res$P_tongue_3,res$P_tongue_4)
cumprobsT < apply(propsT,1,cumsum)
# Just sample one realization
result < sapply(1:length(simtime), function (k) {
test < runif(1) < cumprobsT[,k]
if (all(test)) return(0)
max(which(!test))
})
# Output on realization
result
# Remember! Probabilities are determined deterministically.
# "result" then is sampled based on the probs
Code: Select all
********** MODEL NAME
Provide Model Name
********** MODEL NOTES
Provide Model Information
********** MODEL STATES
# Define your differential equations  if needed
# Here we just use this to replace the "exp(kdecay*t)" in your example
# in some places
d/dt(X) = kdecay*X
X(0) = 1
********** MODEL PARAMETERS
# Define the model parameters
kdecay = 0.01
di0 = 2
ai = 0.3
# Limiting example to tongue only ... you can copy and paste the same things again
# for the other ones  you will get the idea
btongue4 = 1
btongue43 = 1
btongue432 = 1
btongue4321 = 1
********** MODEL VARIABLES
# Determine the logits of your probabilities
logit_P_tongue_ge_1 = di0+ai*(1exp(kdecay*time))+btongue4+btongue43+btongue432+btongue4321
logit_P_tongue_ge_2 = di0+ai*(1exp(kdecay*time))+btongue4+btongue43+btongue432
logit_P_tongue_ge_3 = di0+ai*(1exp(kdecay*time))+btongue4+btongue43
# In the last one we use the ODE instead of the exp(kdecay*time) ... but same result
logit_P_tongue_ge_4 = di0+ai*(1X)+btongue4
# Calculate the probabilities. Note P_tongue_ge_0 by definition of the problem is 1
# Note that: "ge": means "greater or equal"
# P_tongue_ge_0 = 1 # Just a reminder
P_tongue_ge_1 = exp(logit_P_tongue_ge_1)/(1+exp(logit_P_tongue_ge_1))
P_tongue_ge_2 = exp(logit_P_tongue_ge_2)/(1+exp(logit_P_tongue_ge_2))
P_tongue_ge_3 = exp(logit_P_tongue_ge_3)/(1+exp(logit_P_tongue_ge_3))
P_tongue_ge_4 = exp(logit_P_tongue_ge_4)/(1+exp(logit_P_tongue_ge_4))
# Calculate the probabilities for the individual levels
P_tongue_0 = 1P_tongue_ge_1
P_tongue_1 = P_tongue_ge_1  P_tongue_ge_2
P_tongue_2 = P_tongue_ge_2  P_tongue_ge_3
P_tongue_3 = P_tongue_ge_3  P_tongue_ge_4
P_tongue_4 = P_tongue_ge_4
# That's it! Model is defined!
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS