Difference between revisions of "Simulation"
m |
m (→Examples) |
||
(14 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
== Simulation of a hierarchical model == | == Simulation of a hierarchical model == | ||
− | In [[What is a model? A joint probability distribution!]] | + | In the [[What is a model? A joint probability distribution!]] section we proposed to consider a model as a joint probability distribution of random variables. Furthermore, the hierarchical structure of a model allows us to decompose this joint distribution into a product |
of conditional and marginal distributions. | of conditional and marginal distributions. | ||
Line 11: | Line 11: | ||
</math> }} | </math> }} | ||
− | Simulating this model means successively drawing the variables of the model using their associated | + | Simulating this model means successively drawing the variables of the model using their associated submodels: |
Line 36: | Line 36: | ||
Evidently, if certain variables are not considered random and are given, the associated simulation step is not performed. | Evidently, if certain variables are not considered random and are given, the associated simulation step is not performed. | ||
− | We have already seen that for other tasks (model exploration, estimation) | + | We have already seen that for other tasks (model exploration, estimation) $\mlxtran$ is well-adapted to implementing this type of hierarchical model because it can deal with this structure by defining the various submodels in easy-to-read blocks. |
It is also possible to use $\mlxtran$ like a function from within R and Matlab. This capability means that we can write our own R or Matlab script in which the design is defined, the simulation function is called and the results are used for plotting figures, performing tests, etc. This solution combines the flexibility of the R and Matlab environments with the ability of $\mlxtran$ to easily encode complex models. | It is also possible to use $\mlxtran$ like a function from within R and Matlab. This capability means that we can write our own R or Matlab script in which the design is defined, the simulation function is called and the results are used for plotting figures, performing tests, etc. This solution combines the flexibility of the R and Matlab environments with the ability of $\mlxtran$ to easily encode complex models. | ||
Line 43: | Line 43: | ||
== Examples == | == Examples == | ||
− | The joint model ${\cal M}_2$ detailed in the [[Model evaluation#Model selection| Model selection]] | + | The joint model ${\cal M}_2$ detailed in the [[Model evaluation#Model selection| Model selection]] section can be implemented with $\mlxtran$. We suppose that there can be either simple or multiple drug administrations either orally ({{Verbatim|adm{{-}}1}}) or by |
iv ({{Verbatim|adm{{-}}2}}). The model is "ready" to accept these different dose regimens, which are defined either within the data (for estimation) or the simulation script. | iv ({{Verbatim|adm{{-}}2}}). The model is "ready" to accept these different dose regimens, which are defined either within the data (for estimation) or the simulation script. | ||
Line 93: | Line 93: | ||
|name= | |name= | ||
|code=<pre style="background-color: #EFEFEF; border:none"> | |code=<pre style="background-color: #EFEFEF; border:none"> | ||
− | + | %% dose regimen | |
− | admin1 = struct('adm',1,'amount',50,'time',0:24:144); | + | admin1 = struct('adm',1,'amount',50,'time',0:24:144); |
− | + | %% group | |
group1 = struct('size',10); | group1 = struct('size',10); | ||
− | + | %% parameter | |
parameterName = {'ka_pop','V_pop','Cl_pop','h0_pop','omega_ka','omega_V',... | parameterName = {'ka_pop','V_pop','Cl_pop','h0_pop','omega_ka','omega_V',... | ||
'omega_Cl','omega_h0','beta_V','beta_Cl','gamma','a','b'}; | 'omega_Cl','omega_h0','beta_V','beta_Cl','gamma','a','b'}; | ||
Line 105: | Line 105: | ||
param1 = struct('name',parameterName,'value',parameterValue); | param1 = struct('name',parameterName,'value',parameterValue); | ||
− | + | %% output | |
outputName = {'Cc','h','Concentration','Hemorrhaging'}; | outputName = {'Cc','h','Concentration','Hemorrhaging'}; | ||
outputTime = {0:1:200 , 0:1:200 , [1 2 4 6 12 18 28 52 76 100 124 148 172]}; | outputTime = {0:1:200 , 0:1:200 , [1 2 4 6 12 18 28 52 76 100 124 148 172]}; | ||
output1 = struct('name',outputName,'time',outputTime); | output1 = struct('name',outputName,'time',outputTime); | ||
− | + | %% tasks | |
data=simulix('model','joint1_model.txt','administration',admin1,'group',group1,... | data=simulix('model','joint1_model.txt','administration',admin1,'group',group1,... | ||
'parameter',param1,'output',output1); | 'parameter',param1,'output',output1); | ||
Line 124: | Line 124: | ||
|name= | |name= | ||
|code=<pre style="background-color: #EFEFEF; border:none"> | |code=<pre style="background-color: #EFEFEF; border:none"> | ||
− | + | ## dose regimen | |
admin1 <-list(adm = 1, amount = 50, time = seq(0, 144, by=24)) | admin1 <-list(adm = 1, amount = 50, time = seq(0, 144, by=24)) | ||
− | + | ##group | |
group1 <- list(size=10) | group1 <- list(size=10) | ||
− | + | ## parameter | |
parameterName <- c("ka_pop","V_pop","Cl_pop","h0_pop","omega_ka","omega_V", | parameterName <- c("ka_pop","V_pop","Cl_pop","h0_pop","omega_ka","omega_V", | ||
"omega_Cl","omega_h0","beta_V","beta_Cl","gamma","a","b") | "omega_Cl","omega_h0","beta_V","beta_Cl","gamma","a","b") | ||
Line 136: | Line 136: | ||
param1 <- list(name=parameterName, value=parameterValue) | param1 <- list(name=parameterName, value=parameterValue) | ||
− | + | ## output | |
outputName <- c("Cc","h","Concentration","Hemorrhaging"); | outputName <- c("Cc","h","Concentration","Hemorrhaging"); | ||
outputTime <- c(seq(0,200), seq(0,200), c(1 2 4 6 12 18 28 52 76 100 124 148 172)) | outputTime <- c(seq(0,200), seq(0,200), c(1 2 4 6 12 18 28 52 76 100 124 148 172)) | ||
output1 <- list(name=outputName, time=outputTime) | output1 <- list(name=outputName, time=outputTime) | ||
− | + | ## tasks | |
data <- simulix(model="joint1_model.txt",administration=admin1, group=group1, | data <- simulix(model="joint1_model.txt",administration=admin1, group=group1, | ||
parameter=param1,output=output1); | parameter=param1,output=output1); | ||
Line 151: | Line 151: | ||
{{MLXTran | {{MLXTran | ||
− | |name= | + | |name=joint2_model.txt |
|text=<pre style="background-color: #EFEFEF; border:none"> | |text=<pre style="background-color: #EFEFEF; border:none"> | ||
[POPULATION] | [POPULATION] | ||
Line 169: | Line 169: | ||
− | {{ | + | {{Remarks |
|title=Remark | |title=Remark | ||
− | The {{Verbatim|simulix}} function and necessary R and Matlab connectors are still under development and will be available soon. | + | |text=The {{Verbatim|simulix}} function and necessary R and Matlab connectors are still under development and will be available soon. |
}} | }} | ||
+ | |||
+ | |||
+ | {{Back | ||
+ | |link=Model evaluation }} |
Latest revision as of 15:10, 5 June 2013
Simulation of a hierarchical model
In the What is a model? A joint probability distribution! section we proposed to consider a model as a joint probability distribution of random variables. Furthermore, the hierarchical structure of a model allows us to decompose this joint distribution into a product of conditional and marginal distributions.
Thus, if a model includes observations $\by$, individual parameters $\bpsi$, a design made up of observation times $\bt$ and source terms (e.g., doses) $\bu$, individual covariates $\bc$ and population parameters $\theta$, the joint distribution of these random variables can be written:
Simulating this model means successively drawing the variables of the model using their associated submodels:
- The population parameters $\theta$ from the distribution $\qth$.
- The individual covariates $\bc$ from the distribution $\qc$.
- The individual parameters $\bpsi$ from the distribution $\qcpsithc$ using the values of $\theta$ and $\bc$ obtained in steps 1 and 2.
- The dose regimen $\bu$ from the distribution $\qu$.
- The measurement times $\bt$ from the distribution $\qt$.
- Lastly, observations $\by$ from the distribution $\qcypsiut$ using the values of $\bpsi$, $\bu$ and $\bt$ obtained at steps 3, 4 and 5.
Evidently, if certain variables are not considered random and are given, the associated simulation step is not performed.
We have already seen that for other tasks (model exploration, estimation) $\mlxtran$ is well-adapted to implementing this type of hierarchical model because it can deal with this structure by defining the various submodels in easy-to-read blocks. It is also possible to use $\mlxtran$ like a function from within R and Matlab. This capability means that we can write our own R or Matlab script in which the design is defined, the simulation function is called and the results are used for plotting figures, performing tests, etc. This solution combines the flexibility of the R and Matlab environments with the ability of $\mlxtran$ to easily encode complex models.
Examples
The joint model ${\cal M}_2$ detailed in the Model selection section can be implemented with $\mlxtran$. We suppose that there can be either simple or multiple drug administrations either orally (adm=1) or by iv (adm=2). The model is "ready" to accept these different dose regimens, which are defined either within the data (for estimation) or the simulation script.
Here is an example of a Matlab script that uses the joint1_model.txt model to simulate data for 10 patients. In this example, a dose of 50mg is administered every 24 hours over a 7 day period. The simulated data is stored in a cell array. We can then use exploredatapx for visualizing it.
If we preferred, we could instead use R to run exactly the same simulations using exactly the same $\mlxtran$ model:
It is also possible to include in the model a statistical model for the population parameters or the covariates (here, the weight) by adding sections [POPULATION] and/or [COVARIATE] to the $\mlxtran$ model. For example, we could add the following code to the previous $\mlxtran$ code: