hallowkm / RxODE

RxODE is an R package that facilitates easy simulations in R
20 stars 13 forks source link

simulating with variability #19

Open amitvtaneja opened 6 years ago

amitvtaneja commented 6 years ago

Dear Wenping

My query is with reference to the example presented in CPT:PSP Nov 2015 for simulating with variability. In the following example, the 0.09 and 0.08 are the Omegas for the CL and V2, while the 0.08 is the covariance between the two. In this case what is 0.25?

sigma <‐ matrix(c(0.09,0.08,0.08,0.25),2,2) Perhaps you could show how the covariance matrix in NONMEM is structured in this case?

Thanks in advance

Amit Taneja

wwang-at-github commented 6 years ago

Dear Amit:

In your example, the 0.09 and 0.25 are the Omegas for the CL and V2, while the 0.08 & 0.08 are the covariance between the two. Note that this is the FULL symmetric matrix. In NONMEM, it only outputs the lower half. In NONMEM, this same matrix would be denoted: $OMEGA BLOCK(2) 0.09 0.08 0.25

Wenping

On Tue, Mar 20, 2018 at 2:07 PM, amitvtaneja notifications@github.com wrote:

Dear Wenping

My query is with reference to the example presented in CPT:PSP Nov 2015 for simulating with variability. In the following example, the 0.09 and 0.08 are the Omegas for the CL and V2, while the 0.08 is the covariance between the two. In this case what is 0.25?

sigma <‐ matrix(c(0.09,0.08,0.08,0.25),2,2) Perhaps you could show how the covariance matrix in NONMEM is structured in this case?

Thanks in advance

Amit Taneja

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/hallowkm/RxODE/issues/19, or mute the thread https://github.com/notifications/unsubscribe-auth/AMdZdvFI3DGoPL-V1Ai5xiyW57pUlaC5ks5tgUV4gaJpZM4SyZLA .

amitvtaneja commented 6 years ago

Dear Wenping

I tried implementing 2 models with variability. In the first case I used a dataframe of parameters which was output by NONMEM(patab file). In this case, there is also a difference between using the loop and apply for looping thru individual parameters

In the second case I simulated the parameters from the variance covariance matrix. In both cases, however I get an error. Hereby the respective models with the respective error messages. I do hope you will be able to advise me in this matter. Thanking you in advance

Amit

model 1

ode.pk2222 <- "

A 2-compartment model with absorption compartment

d/dt(A1) =-K12A1; d/dt(A2) = K12A1 + K32A3 - (K20+K23)A2; d/dt(A3) = K23A2 - K32A3; "

evtab <- eventTable(amount.units="mg", time.units = "hours") ;#evtab$add.dosing(dose=params$DOSE*params$IF1,dosing.interval=24,nbr.doses=84,start.time=0.43)

evtab$add.dosing(dose=unique(params[params$DOSE==50,]$DOSE)mean(params[params$DOSE==50,]$IF1),dosing.interval=24,nbr.doses=84,start.time=0.43) evtab$add.sampling(8324+seq(0,24,l=101))

params[params$DOSE==50,]->p50 ;# empty parameter matrix for storing results res50<-NULL nsub<-32

for (i in 1:nsub) {

param50<-p50[i,]

x<-rxode.model$solve(param50,evtab,rep(0,3) )

res50<-cbind(res,x[,'A2']) }

===================================================================================================

;## ERROR MESSAGE

Error in data.frame(..., check.names = FALSE) :

;# arguments imply differing number of rows: 1999, 101

In addition: Warning message:

;# In rxInits(dll, inits, state, 0) : Assumed order of inputs: A1, A2, A3

=================================================================================================

;## Use of apply results in a 'strange' output

x<-apply(p50,1,function(th)rxode.model$run(th,evtab,rep(0,3)))

cbind(res50,x)->res50

;### on console

;## Reallocating to 64/128/256

;### output matrix x

;# dim(x)

[1] 404 32

729 730 731 732 745 746 747

[1,] 1992.0000000 1992.0000000 1.992000e+03 1.992000e+03 1992.000000 1992.000000 1992.0000000

[2,] 1992.2400000 1992.2400000 1.992240e+03 1.992240e+03 1992.240000 1992.240000 1992.2400000

[3,] 1992.4800000 1992.4800000 1.992480e+03 1.992480e+03 1992.480000 1992.480000 1992.4800000

[4,] 1992.7200000 1992.7200000 1.992720e+03 1.992720e+03 1992.720000 1992.720000 1992.7200000

[5,] 1992.9600000 1992.9600000 1.992960e+03 1.992960e+03 1992.960000 1992.960000 1992.9600000

[6,] 1993.2000000 1993.2000000 1.993200e+03 1.993200e+03 1993.200000 1993.200000 1993.2000000

[7,] 1993.4400000 1993.4400000 1.993440e+03 1.993440e+03 1993.440000 1993.440000 1993.4400000

[8,] 1993.6800000 1993.6800000 1.993680e+03 1.993680e+03 1993.680000 1993.680000 1993.6800000

[9,] 1993.9200000 1993.9200000 1.993920e+03 1.993920e+03 1993.920000 1993.920000 1993.9200000

[10,] 1994.1600000 1994.1600000 1.994160e+03 1.994160e+03 1994.160000 1994.160000 1994.1600000

;## the first few lines of the dataframe of individual parameters

head(p50)

IF1 ICL IVss IV IKA IQ E0 EMAX IC50 IALAG1 DOSE FORM FOOD FREQ SWCL_ORIG

729 1.0584 22.877 210.74 104.180 0.29577 5.9118 92.828 13.388 37.129 0.43 50 1 0 1 86.0

730 1.0584 22.877 210.74 104.180 0.29577 5.9118 92.828 13.388 37.129 0.43 50 1 0 1 85.5

731 1.0755 22.877 210.74 104.180 0.38270 5.9118 92.828 13.388 37.129 0.43 50 1 0 1 86.0

732 1.0755 22.877 210.74 104.180 0.38270 5.9118 92.828 13.388 37.129 0.43 50 1 0 1 85.5

745 1.0976 13.942 179.96 79.148 0.15116 6.2829 101.980 18.294 37.129 0.43 50 1 0 1 84.0

746 1.0976 13.942 179.96 79.148 0.15116 6.2829 101.980 18.294 37.129 0.43 50 1 0 1 91.5

K20 K12 K23 K32 AUC.TRUE SID

729 0.2195911 0.29577 0.02805258 0.02805258 2.313240 179

730 0.2195911 0.29577 0.02805258 0.02805258 2.313240 180

731 0.2195911 0.38270 0.02805258 0.02805258 2.350614 181

732 0.2195911 0.38270 0.02805258 0.02805258 2.350614 182

745 0.1761510 0.15116 0.03491276 0.03491276 3.936308 183

746 0.1761510 0.15116 0.03491276 0.03491276 3.936308 184

#####################################model 2#####################################

m1 <- RxODE(model = ode)

;# create dosing and observation (sampling) events ;# QD (once daily) dosing, 5 days

qd <- eventTable(amount.units="umg", time.units = "hours") qd$add.dosing(dose=10000, nbr.doses=5, dosing.interval = 168)

hourly sampling

qd$add.sampling(seq(from=0, to=2600, by=0.1) )

bid <- eventTable(amount.units="ug", time.units="days") # only dosing

bid$add.dosing(dose = 10000,

;# nbr.doses = 2*5, ;# dosing.interval = 12)

bid$add.sampling(0:24)

bid$add.sampling(96+0:24)

;# n subproblems

nsub<-1000 # 1000 subproblems ;# Cl,V2 sigma1<-matrix(c(0.0655334,0,0,0.020843),2,2)

Q,V3

sigma2<-matrix(c(0,0,0,0.077613),2,2) ;# KIn, Kout sigma3<-matrix(c(0,0,0,0.167372),2,2)

;# sample from the covariance matrix mvrnorm(n=nsub,rep(0,2),sigma1)->mv1 mvrnorm(n=nsub,rep(0,2),sigma2)->mv2 mvrnorm(n=nsub,rep(0,2),sigma3)->mv3

params<-cbind(CL=0.19, V2=51.28, Q=0.23, V3=28.69, Kin=0.0007727.09, Kout=0.00077, EC50=5.160.147)

ICL<-params[1]exp(mv1[,1]) IV2<-params[2]exp(mv1[,2]) IQ<-params[3]exp(mv2[,1]) IV3<-params[4]exp(mv2[,2]) IKin<-params[5]exp(mv3[,1]) IKout<-params[6]exp(mv3[,2]) EC50<-rep(params[7],nsub)

params.all<-cbind(CL=ICL,V2=IV2,V3=IV3,Q=IQ,Kin=IKin,Kout=IKout,EC50=EC50) res<-NULL

;## loop thru each row of parameter values

for (i in 1:nsub){

params<-params.all[i,] x<-m1$solve(params,qd,inits=c(0,0,27.09))

res<-cbind(res,x) }

res<-apply(params.all,1,function(theta)m1$run(theta,qd,inits = c(0,0,27.09)))

ERROR MESSAGE======================================================================================

;## consistent with use of both the loop and the apply function

Error in m1$run(theta, qd, inits = c(0, 0, 27.09)) :

;# Tried both LSODA and DOP853, but could not solve the system.

In addition: Warning message:

;# In rxInits(dll, inits, state, 0) :

=====================================================================================================