multiple categorical response simulation

IQR Tools supports and increases the efficiency, quality, and compliance of model-based analyses in pharmacometrics, systems parmacology, and systems biology by incorporating and extending the capabilities of existing tools. The user-friendliness of the package considerably lowers the threshold for the conduct of modeling & simulation based analyses, which also makes it useful both for complex analyses and for educational purposes.

Here you can ask questions and get answers ...

Moderators: daniel.kaschek, intiquan

Post Reply
poppharm@gmail.com
Posts: 2
Joined: Thu Jun 20, 2019 10:20 am

multiple categorical response simulation

Post by poppharm@gmail.com » Thu Jun 20, 2019 10:39 am

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*(1-exp(-kdecay*t))+btongue4+btongue43+btongue432+btongue4321
logit(P(tongue >=2)) = di0+ai*(1-exp(-kdecay*t))+btongue4+btongue43+btongue432
logit(P(tongue >=3)) = di0+ai*(1-exp(-kdecay*t))+btongue4+btongue43
logit(P(tongue >=4)) = di0+ai*(1-exp(-kdecay*t))+btongue4


}
upper = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(upper>=1)) = di0+ai*(1-exp(-kdecay*t))+bupper4+bupper43+bupper432+bupper4321
logit(P(upper>=2)) = di0+ai*(1-exp(-kdecay*t))+bupper4+bupper43+bupper432
logit(P(upper>=3)) = di0+ai*(1-exp(-kdecay*t))+bupper4+bupper43
logit(P(upper>=4)) = di0+ai*(1-exp(-kdecay*t))+bupper4

}


lower = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(lower>=1)) = di0+ai*(1-exp(-kdecay*t))+blower4+blower43+blower432+blower4321
logit(P(lower>=2)) = di0+ai*(1-exp(-kdecay*t))+blower4+blower43+blower432
logit(P(lower>=3)) = di0+ai*(1-exp(-kdecay*t))+blower4+blower43
logit(P(lower>=4)) = di0+ai*(1-exp(-kdecay*t))+blower4
}



jaw = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(jaw>=1)) = di0+ai*(1-exp(-kdecay*t))+bjaw4+bjaw43+bjaw432+bjaw4321
logit(P(jaw>=2)) = di0+ai*(1-exp(-kdecay*t))+bjaw4+bjaw43+bjaw432
logit(P(jaw>=3)) = di0+ai*(1-exp(-kdecay*t))+bjaw4+bjaw43
logit(P(jaw>=4)) = di0+ai*(1-exp(-kdecay*t))+bjaw4
}

lips = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(lips>=1)) = di0+ai*(1-exp(-kdecay*t))+blips4+blips43+blips432+blips4321
logit(P(lips>=2)) = di0+ai*(1-exp(-kdecay*t))+blips4+blips43+blips432
logit(P(lips>=3)) = di0+ai*(1-exp(-kdecay*t))+blips4+blips43
logit(P(lips>=4)) = di0+ai*(1-exp(-kdecay*t))+blips4
}


muscles = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(muscles>=1)) = di0+ai*(1-exp(-kdecay*t))+bmuscles4+bmuscles43+bmuscles432+bmuscles4321
logit(P(muscles>=2)) = di0+ai*(1-exp(-kdecay*t))+bmuscles4+bmuscles43+bmuscles432
logit(P(muscles>=3)) = di0+ai*(1-exp(-kdecay*t))+bmuscles4+bmuscles43
logit(P(muscles>=4)) = di0+ai*(1-exp(-kdecay*t))+bmuscles4
}


neck = { type = categorical, categories = {0, 1, 2, 3,4},
logit(P(neck>=1)) = di0+ai*(1-exp(-kdecay*t))+bneck4+bneck43+bneck432+bneck4321
logit(P(neck>=2)) = di0+ai*(1-exp(-kdecay*t))+bneck4+bneck43+bneck432
logit(P(neck>=3)) = di0+ai*(1-exp(-kdecay*t))+bneck4+bneck43
logit(P(neck>=4)) = di0+ai*(1-exp(-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
intiquan
Site Admin
Posts: 7
Joined: Mon Sep 03, 2018 11:50 am
Location: Basel, Switzerland
Contact:

Re: multiple categorical response simulation

Post by intiquan » Thu Jul 04, 2019 10:48 am

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

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*(1-X)
 
# 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

 
intiquan
Site Admin
Posts: 7
Joined: Mon Sep 03, 2018 11:50 am
Location: Basel, Switzerland
Contact:

Re: multiple categorical response simulation

Post by intiquan » Fri Jul 05, 2019 9:52 am

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

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
The example.txt file contents

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*(1-exp(-kdecay*time))+btongue4+btongue43+btongue432+btongue4321
logit_P_tongue_ge_2 = di0+ai*(1-exp(-kdecay*time))+btongue4+btongue43+btongue432
logit_P_tongue_ge_3 = di0+ai*(1-exp(-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*(1-X)+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 = 1-P_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

Post Reply