Re: [R] deSolve question
Dear Andras
In algebra (and computing) terms are operated upon in order of priority. Terms
in brackets are evaluated first. In your original code, - pars$k * state[1] is
not divided by pars$v since division happens before addition. When you use the
brackets that expression becomes part of the term that is divided.
Therefore the outputs of these functions should not be the same.
Best wishes
Chris
Chris Campbell, PhD
Tel. +44 (0) 1249 705 450 | Mobile. +44 (0) 7929 628 349
mailto:ccampb...@mango-solutions.com | http://www.mango-solutions.com
Mango Solutions, 2 Methuen Park, Chippenham, Wiltshire , SN14 OGB UK
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Andras Farkas
Sent: 19 June 2013 00:46
To: r-help@r-project.org; r-sig-dynamic-mod...@r-project.org
Subject: [R] deSolve question
Dear All
wonder if you could provide some insights on the following: currently I have
this code which produces the expected results:
require(deSolve)
pars - list(k = 0.08,v=15)
intimes - c(0,0.5,12)
input - c(800,0,0)
forc - approxfun(intimes, input, method=constant, rule=2)
derivs - function(t, state, pars) {
inp - forc(t)
dy1 - - pars$k * state[1] + inp/pars$v
return(list(c(dy1)))
}
model - function(pars, times=seq(0, 13, by = 0.01)) {
state - c(y = 0)
return(ode(y = state, times = times, func = derivs, parms = pars,
method=lsoda))
}
plot(model(pars))
intuitively it would make sense to me that if I change the derivs as:
derivs - function(t, state, pars) {
inp - forc(t)
#note the added ()
dy1 - (- pars$k * state[1] + inp)/pars$v
return(list(c(dy1)))
}
than I would get the same results, but that is not the case. I need to
relocate pars$v to another position within derivs so that I could implement
it in a forcing function to be able to change its value during derivation.
Appreciate your help,
Andras
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
LEGAL NOTICE\ \ This message is intended for the use of ...{{dropped:18}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] deSolve question
Dear All
wonder if you could provide some insights on the following: currently I have
this code which produces the expected results:
require(deSolve)
pars - list(k = 0.08,v=15)
intimes - c(0,0.5,12)
input - c(800,0,0)
forc - approxfun(intimes, input, method=constant, rule=2)
derivs - function(t, state, pars) {
inp - forc(t)
dy1 - - pars$k * state[1] + inp/pars$v
return(list(c(dy1)))
}
model - function(pars, times=seq(0, 13, by = 0.01)) {
state - c(y = 0)
return(ode(y = state, times = times, func = derivs, parms = pars,
method=lsoda))
}
plot(model(pars))
intuitively it would make sense to me that if I change the derivs as:
derivs - function(t, state, pars) {
inp - forc(t)
#note the added ()
dy1 - (- pars$k * state[1] + inp)/pars$v
return(list(c(dy1)))
}
than I would get the same results, but that is not the case. I need to
relocate pars$v to another position within derivs so that I could implement
it in a forcing function to be able to change its value during derivation.
Appreciate your help,
Andras
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] deSolve question
insun nam wrote: Dear All, I like to simulate a physiologically based pharmacokinetics model using R but am having a problem with the daspk routine. The same problem has been implemented in Berkeley madonna and Winbugs so that I know that it is working. However, with daspk it is not, and the numbers are everywhere! Please see the following and let me know if I am missing something... Thanks a lot in advance, In-Sun [ ... long example deleted, see original post ...] Hi In-Sun, as far as I see the *first* of your problems is your high demand on precision. Double precision allows 16 digits, so 1e-10 may lead internally to values that are close to the maximum number of digits in double precision arithmetics, because you set relative *and* absolute tolerance of all states to such a small value. Reduce atol and rtol to reasonable values (e.g. the default of 1e-6) and it should work (it does on my machine). Another possibility would be to increase maxsteps, but this will not help if the requested precision is higher than would be possible with double precision. If you try the following: ODE - as.data.frame(daspk(y = y, times = times, func = Fun_ODE, parms = pars, atol = 1e-6, rtol = 1e-6)) ... then the simulation goes through in a technical sense, however state variables reach very high values which, again, are beyond what is possible with double precision arithmetics -- is this what you mean with the numbers are everywhere? Are you sure that your model formulation is correct? Please check your equations, especially parentheses and signs. Thomas Petzoldt __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] deSolve question
In-Sun,
It is very simple. You define your state variables in the following order:
y - c(Agi = 0,Alu = 0, Abr = 0, Ah = 0, Ali = 0, Ak = 0, Am = 0, Ask =
0, Aad = 0, Apa = 0, Asp = 0, Aar = 0, Ave = 0)
and your rates of change in another order:
dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi, dAk,
dAm, dAad, dAsk)
So, The solver uses the rate of change of Aar to update Agi etc... and you
get nonsense in the end.
I hope correcting this will solve your problem.
I have seen this multiple times; it is definitely the most common type of
mistake in using R for solving differential equations.
Btw, why use daspk ? It is really meant for solving DAEs, not ODEs.
lsoda or lsode or vode might be a better choice.
Cheers,
Karline
insun nam wrote:
Dear All,
I like to simulate a physiologically based pharmacokinetics model using R
but am having a problem with the daspk routine.
The same problem has been implemented in Berkeley madonna and Winbugs so
that I know that it is working. However, with daspk it is not, and the
numbers are everywhere!
Please see the following and let me know if I am missing something...
Thanks a lot in advance,
In-Sun
#---
library(deSolve)
y - c(Agi = 0,Alu = 0, Abr = 0, Ah = 0, Ali = 0, Ak = 0, Am = 0, Ask
=
0, Aad = 0, Apa = 0, Asp = 0, Aar = 0, Ave = 0)
times = seq(0, 100, length=100)
pars - c(
dose = 80 * 0.26,
doseduration = 10,
Vmax = 1.44,
Km = 0.96,
s = 1.33,
fp = 0.236,
Kpfgi=0.324,
Kpflu = 1.092,
Kpfbr= 0.155 ,
Kpfh=0.767,
Kpfli = 0.551,
Kpfk=0.537,
Kpfm=0.339,
Kpfsk=0.784,
Kpfad=0.465,
Kpfpa=0.595,
Kpfsp=0.410,
Qar = 51.9,
Qve = 51.9,
Qgi = 12.3,
Qlu = 51.9,
Qbr = 3.2,
Qh = 6.4,
Qli = 16.5,
Qk = 12.8,
Qm = 7.6,
Qsk = 5.0,
Qad = 0.4,
Qpa = 1.0,
Qsp = 1.0,
Var = 7.0,
Vve = 14.1,
Vgi = 12.4,
Vlu = 1.3,
Vbr = 1.3,
Vh = 1.2,
Vli = 12.4,
Vk = 2.2,
Vm = 140.0,
Vsk = 49.0,
Vad = 11.2,
Vpa = 1.0,
Vsp = 1.0
)
Fun_ODE - function(t,y, pars){
with (as.list(c(y, pars)), {
It - dose/doseduration
Car - Aar/Var
Cve - Ave/Vve
Clu - Alu/Vlu
Cli - Ali/Vli
Cbr - Abr/Vbr
Ch - Ah/Vh
Cpa - Apa/Vpa
Csp - Asp/Vsp
Cgi - Agi/Vgi
Ck - Ak/Vk
Cm - Am/Vm
Cad - Aad/Vad
Csk - Ask/Vsk
Kpbbr - s*fp*Kpfbr
Kpbli - s*fp*Kpfli
Kpbh - s*fp*Kpfh
Kpbpa - s*fp*Kpfpa
Kpbsp - s*fp*Kpfsp
Kpbgi - s*fp*Kpfgi
Kpbk - s*fp*Kpfk
Kpbm - s*fp*Kpfm
Kpbad - s*fp*Kpfad
Kpbsk - s*fp*Kpfsk
Kpblu - s*fp*Kpflu
dAar - (Clu/Kpblu - Car)*Qar
dAve - if (t 10) It + Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli +
Ck*Qk/Kpbk + Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
else Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli + Ck*Qk/Kpbk +
Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
dAlu - (Cve-Clu/Kpblu)*Qlu
dAli - ((Qli - Qgi- Qpa-Qsp)*Car + Cgi*Qgi/Kpbgi + Csp*Qsp/Kpbsp +
Cpa*Qpa/Kpbpa - Cli*Qli/Kpbli) - Vmax*Cli/Kpfli/(Km + Cli/Kpfli)
dAbr - (Car - Cbr/Kpbbr)*Qbr
dAh - (Car - Ch/Kpbh)*Qh
dApa - (Car - Cpa/Kpbpa)*Qpa
dAsp - (Car - Csp/Kpbsp)*Qsp
dAgi - (Car - Cgi/Kpbgi)*Qgi
dAk - (Car - Ck/Kpbk)*Qk
dAm - (Car - Cm/Kpbm)*Qm
dAad - (Car - Cad/Kpbad)*Qad
dAsk - (Car - Csk/Kpbsk)*Qsk
return(list(dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi,
dAk,
dAm, dAad, dAsk),
Car = Car, Cve=Cve, Clu=Clu, Cli=Cli, Cbr=Cbr, Ch=Ch, Cpa=Cpa,
Csp=Csp, Cgi=Cgi, Ck=Ck, Cm=Cm, Cad=Cad, Csk=Csk))
})
}
ODE - as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
parms = pars, atol = 1e-10, rtol = 1e-10))
--
Dr In-Sun Nam Knutsson
Research Associate
The Centre for Applied Pharmacokinetic Research (CAPKR)
School of Pharmacy and Pharmaceutical Sciences
University of Manchester
Stopford Building
Oxford Road
Manchester
U.K.
Phone: +44 161 275 2355
Email: in-sunnam.knuts...@manchester.ac.uk
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
View this message in context:
http://www.nabble.com/deSolve-question-tp23985008p23994033.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] deSolve question
Dear All,
I like to simulate a physiologically based pharmacokinetics model using R
but am having a problem with the daspk routine.
The same problem has been implemented in Berkeley madonna and Winbugs so
that I know that it is working. However, with daspk it is not, and the
numbers are everywhere!
Please see the following and let me know if I am missing something...
Thanks a lot in advance,
In-Sun
#---
library(deSolve)
y - c(Agi = 0,Alu = 0, Abr = 0, Ah = 0, Ali = 0, Ak = 0, Am = 0, Ask =
0, Aad = 0, Apa = 0, Asp = 0, Aar = 0, Ave = 0)
times = seq(0, 100, length=100)
pars - c(
dose = 80 * 0.26,
doseduration = 10,
Vmax = 1.44,
Km = 0.96,
s = 1.33,
fp = 0.236,
Kpfgi=0.324,
Kpflu = 1.092,
Kpfbr= 0.155 ,
Kpfh=0.767,
Kpfli = 0.551,
Kpfk=0.537,
Kpfm=0.339,
Kpfsk=0.784,
Kpfad=0.465,
Kpfpa=0.595,
Kpfsp=0.410,
Qar = 51.9,
Qve = 51.9,
Qgi = 12.3,
Qlu = 51.9,
Qbr = 3.2,
Qh = 6.4,
Qli = 16.5,
Qk = 12.8,
Qm = 7.6,
Qsk = 5.0,
Qad = 0.4,
Qpa = 1.0,
Qsp = 1.0,
Var = 7.0,
Vve = 14.1,
Vgi = 12.4,
Vlu = 1.3,
Vbr = 1.3,
Vh = 1.2,
Vli = 12.4,
Vk = 2.2,
Vm = 140.0,
Vsk = 49.0,
Vad = 11.2,
Vpa = 1.0,
Vsp = 1.0
)
Fun_ODE - function(t,y, pars){
with (as.list(c(y, pars)), {
It - dose/doseduration
Car - Aar/Var
Cve - Ave/Vve
Clu - Alu/Vlu
Cli - Ali/Vli
Cbr - Abr/Vbr
Ch - Ah/Vh
Cpa - Apa/Vpa
Csp - Asp/Vsp
Cgi - Agi/Vgi
Ck - Ak/Vk
Cm - Am/Vm
Cad - Aad/Vad
Csk - Ask/Vsk
Kpbbr - s*fp*Kpfbr
Kpbli - s*fp*Kpfli
Kpbh - s*fp*Kpfh
Kpbpa - s*fp*Kpfpa
Kpbsp - s*fp*Kpfsp
Kpbgi - s*fp*Kpfgi
Kpbk - s*fp*Kpfk
Kpbm - s*fp*Kpfm
Kpbad - s*fp*Kpfad
Kpbsk - s*fp*Kpfsk
Kpblu - s*fp*Kpflu
dAar - (Clu/Kpblu - Car)*Qar
dAve - if (t 10) It + Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli +
Ck*Qk/Kpbk + Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
else Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli + Ck*Qk/Kpbk +
Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
dAlu - (Cve-Clu/Kpblu)*Qlu
dAli - ((Qli - Qgi- Qpa-Qsp)*Car + Cgi*Qgi/Kpbgi + Csp*Qsp/Kpbsp +
Cpa*Qpa/Kpbpa - Cli*Qli/Kpbli) - Vmax*Cli/Kpfli/(Km + Cli/Kpfli)
dAbr - (Car - Cbr/Kpbbr)*Qbr
dAh - (Car - Ch/Kpbh)*Qh
dApa - (Car - Cpa/Kpbpa)*Qpa
dAsp - (Car - Csp/Kpbsp)*Qsp
dAgi - (Car - Cgi/Kpbgi)*Qgi
dAk - (Car - Ck/Kpbk)*Qk
dAm - (Car - Cm/Kpbm)*Qm
dAad - (Car - Cad/Kpbad)*Qad
dAsk - (Car - Csk/Kpbsk)*Qsk
return(list(dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi, dAk,
dAm, dAad, dAsk),
Car = Car, Cve=Cve, Clu=Clu, Cli=Cli, Cbr=Cbr, Ch=Ch, Cpa=Cpa,
Csp=Csp, Cgi=Cgi, Ck=Ck, Cm=Cm, Cad=Cad, Csk=Csk))
})
}
ODE - as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
parms = pars, atol = 1e-10, rtol = 1e-10))
--
Dr In-Sun Nam Knutsson
Research Associate
The Centre for Applied Pharmacokinetic Research (CAPKR)
School of Pharmacy and Pharmaceutical Sciences
University of Manchester
Stopford Building
Oxford Road
Manchester
U.K.
Phone: +44 161 275 2355
Email: in-sunnam.knuts...@manchester.ac.uk
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] deSolve question
Dear In-Sun:
I have not seen a reply, so I will offer some suggestions, hoping
I can help without understanding all the details.
1. Have you run that code with options(warn=2)? It
produced over 50 warnings for me, and options(warn=2) will convert the
warnings to errors, making it easier for you to find the exact problem
lines of code. With this, have you used the debug function to trace
through your code line by line to provide more detail about what the
code is doing? I use that routinely.
2. Have you considered the PKfit package? I don't know if
it includes your model, but it includes code for many pharmacokinetic
models. If your interest in deSolve is as a means to solve this
problem, you might consider PKfit.
3. Have you tried writing directly to the authors? Names
and email addresses are available in help(pac=deSolve).
Hope this helps.
Spencer Graves
insun nam wrote:
Dear All,
I like to simulate a physiologically based pharmacokinetics model using R
but am having a problem with the daspk routine.
The same problem has been implemented in Berkeley madonna and Winbugs so
that I know that it is working. However, with daspk it is not, and the
numbers are everywhere!
Please see the following and let me know if I am missing something...
Thanks a lot in advance,
In-Sun
#---
library(deSolve)
y - c(Agi = 0,Alu = 0, Abr = 0, Ah = 0, Ali = 0, Ak = 0, Am = 0, Ask =
0, Aad = 0, Apa = 0, Asp = 0, Aar = 0, Ave = 0)
times = seq(0, 100, length=100)
pars - c(
dose = 80 * 0.26,
doseduration = 10,
Vmax = 1.44,
Km = 0.96,
s = 1.33,
fp = 0.236,
Kpfgi=0.324,
Kpflu = 1.092,
Kpfbr= 0.155 ,
Kpfh=0.767,
Kpfli = 0.551,
Kpfk=0.537,
Kpfm=0.339,
Kpfsk=0.784,
Kpfad=0.465,
Kpfpa=0.595,
Kpfsp=0.410,
Qar = 51.9,
Qve = 51.9,
Qgi = 12.3,
Qlu = 51.9,
Qbr = 3.2,
Qh = 6.4,
Qli = 16.5,
Qk = 12.8,
Qm = 7.6,
Qsk = 5.0,
Qad = 0.4,
Qpa = 1.0,
Qsp = 1.0,
Var = 7.0,
Vve = 14.1,
Vgi = 12.4,
Vlu = 1.3,
Vbr = 1.3,
Vh = 1.2,
Vli = 12.4,
Vk = 2.2,
Vm = 140.0,
Vsk = 49.0,
Vad = 11.2,
Vpa = 1.0,
Vsp = 1.0
)
Fun_ODE - function(t,y, pars){
with (as.list(c(y, pars)), {
It - dose/doseduration
Car - Aar/Var
Cve - Ave/Vve
Clu - Alu/Vlu
Cli - Ali/Vli
Cbr - Abr/Vbr
Ch - Ah/Vh
Cpa - Apa/Vpa
Csp - Asp/Vsp
Cgi - Agi/Vgi
Ck - Ak/Vk
Cm - Am/Vm
Cad - Aad/Vad
Csk - Ask/Vsk
Kpbbr - s*fp*Kpfbr
Kpbli - s*fp*Kpfli
Kpbh - s*fp*Kpfh
Kpbpa - s*fp*Kpfpa
Kpbsp - s*fp*Kpfsp
Kpbgi - s*fp*Kpfgi
Kpbk - s*fp*Kpfk
Kpbm - s*fp*Kpfm
Kpbad - s*fp*Kpfad
Kpbsk - s*fp*Kpfsk
Kpblu - s*fp*Kpflu
dAar - (Clu/Kpblu - Car)*Qar
dAve - if (t 10) It + Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli +
Ck*Qk/Kpbk + Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
else Cbr*Qbr/Kpbbr + Ch *Qh/Kpbh + Cli*Qli/Kpbli + Ck*Qk/Kpbk +
Cm*Qm/Kpbm + Csk * Qsk /Kpbsk + Cad*Qad/Kpbad - Cve*Qve
dAlu - (Cve-Clu/Kpblu)*Qlu
dAli - ((Qli - Qgi- Qpa-Qsp)*Car + Cgi*Qgi/Kpbgi + Csp*Qsp/Kpbsp +
Cpa*Qpa/Kpbpa - Cli*Qli/Kpbli) - Vmax*Cli/Kpfli/(Km + Cli/Kpfli)
dAbr - (Car - Cbr/Kpbbr)*Qbr
dAh - (Car - Ch/Kpbh)*Qh
dApa - (Car - Cpa/Kpbpa)*Qpa
dAsp - (Car - Csp/Kpbsp)*Qsp
dAgi - (Car - Cgi/Kpbgi)*Qgi
dAk - (Car - Ck/Kpbk)*Qk
dAm - (Car - Cm/Kpbm)*Qm
dAad - (Car - Cad/Kpbad)*Qad
dAsk - (Car - Csk/Kpbsk)*Qsk
return(list(dy = c(dAar, dAve, dAlu, dAli, dAbr, dAh, dApa, dAsp, dAgi, dAk,
dAm, dAad, dAsk),
Car = Car, Cve=Cve, Clu=Clu, Cli=Cli, Cbr=Cbr, Ch=Ch, Cpa=Cpa,
Csp=Csp, Cgi=Cgi, Ck=Ck, Cm=Cm, Cad=Cad, Csk=Csk))
})
}
ODE - as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
parms = pars, atol = 1e-10, rtol = 1e-10))
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
