Skip to content

Error with nls

6 messages · 1-27206531-0-90000491, Shawn Hornsby, Douglas Bates

#
Dear R-group members,
 
 I use:
 
 platform i386-pc-mingw32
 arch     x86            
 os       Win32          
 system   x86, Win32     
 status                  
 major    1              
 minor    4.1            
 year     2002           
 month    01             
 day      30             
 language R
 
 I try to fit a 2 compartment model. The compartments are open, connected
 to each other and are filled via constant input and a time depended
 function as well. Data describes increasing of Apo B after dialysis. Aim
 of the analysis is to test the hypothesis whether the data could described
 by two simple disconnected one compartment modes ore the "saturated
 model" holds? The first order differential equation for the saturated
 model:
 
 db5 = - (k50+k56)*b5 + k56*b6 + c*g(t) + h
 db6 = + k65*b5 - (k60+k65)*b6 + d
 
 db5, db6 are the first derivatives, b5, b6 are the functions to be
 fitted. The remaining parameters are unknown and should follow from the
 fit.
 
 assuming that g(t) has the functional form: b4i + (b40-b4i)*exp(-k4*t)
 
 (after calculations of 2 papers of A4) follows the solution:
 
 L5L6 <- function(b40, b4i, k4, t, p50, p56, p60, p65, pc, ph, pd, pb50,
 pb60) {
 
 	k50 <- exp(p50)
 	k56 <- exp(p56)
 	k60 <- exp(p60)
 	k65 <- exp(p65)
 	c   <- exp(pc)
 	h   <- exp(ph)
 	d   <- exp(pd)
 	b50 <- exp(pb50)
 	b60 <- exp(pb60)
 	a <- (k50+k56)
 	b <- k65
 	e <- k56
 	f <- (k60+k65)
 	z1 <- (-(a+f)/2 - sqrt((a+f)^2/4 - a*f + b*e))
 	z2 <- (-(a+f)/2 + sqrt((a+f)^2/4 - a*f + b*e))
 	K  <- ((z1+a)/(z2-z1))
 	B1 <- (b/(z2-z1)*b60 - K*b50)
 	A1 <- (b50-B1)
 	X1 <- (b*d/(z2-z1)-K*(c*b4i+h))
 	X2 <- (K*c*(b4i-b40))
 	X3 <- (c*b4i + h - X1)
 	X4 <- (c*(b40-b4i)- X2)
 	C1E <- (X3/(-z1)*(1-exp(z1*t)) +
 X4/(-(k4+z1))*(exp(-k4*t)-exp(z1*t)))
 	C2E <- (X1/(-z2)*(1-exp(z2*t)) +
 X2/(-(k4+z2))*(exp(-k4*t)-exp(z2*t)))
 	b5 <- (A1*exp(z1*t) + B1*exp(z2*t) + C1E + C2E)
 	b6 <- ((z1+a)/b * A1*exp(z1*t) + (z2+a)/b * B1*exp(z2*t) +
 (z1+a)/b * C1E + (z2+a)/b * C2E)
 	y <- f5*b5 + f6*b6
 	return(y)
 }
 
 I am in the lucky circumstances having starting values, because a nlr-fit
 succeeds, the graphical presentation of the fits looks quite nice. The nlr
 function is part of Lindsey's library(gnlm), but now I would like to apply
 Pinheiro and Bates library(nlme) and I have got an error:
If somebody feel that he can help me, I could send him my R- code and
 data file as well.
 
 Kind regards,
 
 Dominik
 

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
Dear Kind,

I would like to let you know up front that I am not a mathematician nor do I
want to insult you or your intelligence. I am a humble person of simple
education and means and I am offering a suggestion, you may have already
resolved this issue. Here are my suggestions:

In the expression, db5 = - (k50+k56)*b5 + k56*b6 + c*g(t) + h, I notice that
function g(t) is not explicitly defined as an expression. I do see you make
reference to it in, assuming that g(t) has the functional form: b4i +
(b40-b4i)*exp(-k4*t), maybe this is the expression and I am overlooking it.
While, I continued to examine the function g(t), I saw no explicit
definition of variables b40, b4i, k4, and t. They are shown to be part of
the equation but they are not explicitly defined with values.

Again, I am hoping my observation spark an idea that would lead you to your
resolve.

If there is someone in the R-Project that has already helped Kind, I would
like to thank you.

Regards,
Shawn

-----Original Message-----
From: owner-r-help at stat.math.ethz.ch
[mailto:owner-r-help at stat.math.ethz.ch]On Behalf Of
1-27206531-0-90000491
Sent: Wednesday, March 27, 2002 4:48 AM
To: r-help at stat.math.ethz.ch
Subject: [R] Error with nls


 Dear R-group members,

 I use:

 platform i386-pc-mingw32
 arch     x86
 os       Win32
 system   x86, Win32
 status
 major    1
 minor    4.1
 year     2002
 month    01
 day      30
 language R

 I try to fit a 2 compartment model. The compartments are open, connected
 to each other and are filled via constant input and a time depended
 function as well. Data describes increasing of Apo B after dialysis. Aim
 of the analysis is to test the hypothesis whether the data could described
 by two simple disconnected one compartment modes ore the "saturated
 model" holds? The first order differential equation for the saturated
 model:

 db5 = - (k50+k56)*b5 + k56*b6 + c*g(t) + h
 db6 = + k65*b5 - (k60+k65)*b6 + d

 db5, db6 are the first derivatives, b5, b6 are the functions to be
 fitted. The remaining parameters are unknown and should follow from the
 fit.

 assuming that g(t) has the functional form: b4i + (b40-b4i)*exp(-k4*t)

 (after calculations of 2 papers of A4) follows the solution:

 L5L6 <- function(b40, b4i, k4, t, p50, p56, p60, p65, pc, ph, pd, pb50,
 pb60) {

 	k50 <- exp(p50)
 	k56 <- exp(p56)
 	k60 <- exp(p60)
 	k65 <- exp(p65)
 	c   <- exp(pc)
 	h   <- exp(ph)
 	d   <- exp(pd)
 	b50 <- exp(pb50)
 	b60 <- exp(pb60)
 	a <- (k50+k56)
 	b <- k65
 	e <- k56
 	f <- (k60+k65)
 	z1 <- (-(a+f)/2 - sqrt((a+f)^2/4 - a*f + b*e))
 	z2 <- (-(a+f)/2 + sqrt((a+f)^2/4 - a*f + b*e))
 	K  <- ((z1+a)/(z2-z1))
 	B1 <- (b/(z2-z1)*b60 - K*b50)
 	A1 <- (b50-B1)
 	X1 <- (b*d/(z2-z1)-K*(c*b4i+h))
 	X2 <- (K*c*(b4i-b40))
 	X3 <- (c*b4i + h - X1)
 	X4 <- (c*(b40-b4i)- X2)
 	C1E <- (X3/(-z1)*(1-exp(z1*t)) +
 X4/(-(k4+z1))*(exp(-k4*t)-exp(z1*t)))
 	C2E <- (X1/(-z2)*(1-exp(z2*t)) +
 X2/(-(k4+z2))*(exp(-k4*t)-exp(z2*t)))
 	b5 <- (A1*exp(z1*t) + B1*exp(z2*t) + C1E + C2E)
 	b6 <- ((z1+a)/b * A1*exp(z1*t) + (z2+a)/b * B1*exp(z2*t) +
 (z1+a)/b * C1E + (z2+a)/b * C2E)
 	y <- f5*b5 + f6*b6
 	return(y)
 }

 I am in the lucky circumstances having starting values, because a nlr-fit
 succeeds, the graphical presentation of the fits looks quite nice. The nlr
 function is part of Lindsey's library(gnlm), but now I would like to apply
 Pinheiro and Bates library(nlme) and I have got an error:
If somebody feel that he can help me, I could send him my R- code and
 data file as well.

 Kind regards,

 Dominik


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
1-27206531-0-90000491 <domi at sun11.ukl.uni-freiburg.de> writes:
Thank you for providing that information.
It is likely that the iterative algorithm is progressing to values of
the parameters that don't make sense physically.  I suggest that you
add trace = TRUE to your call to nls.  This will provide a record of the
parameter values, the residual sum of squares, and the convergence
criterion throughout the iterations.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
1 day later
#
On Wed, 27 Mar 2002, Shawn Hornsby wrote:

            
In the origin, it was a compartment model with 6 compartments. Data was
available for every compartment. But of scientific interest are specially
the last 2 compartments, #5 and #6. The function g(t) is the approach to
fill compartment  5 via a forcing function. The parameters b4i and b40 and
k4 results from a fit of compartment 4. In my model of compartment 5 and 6
the parameters b40, b4i and k4 are like time constant covariables.
Greetings,

Dominik

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
Thank for the advice.

I repeated the call of nls with option: trace=TRUE:
pb50, pb60), 
+ data=help, start=c(p50=0.008678954, p56=-0.595153967, p60=-4.602990518,
p65=-0.625732096, 
+ pc=-0.128657978, ph=0.708033556, pd=1.140357461, pb50=1.311141424,
pb60=1.270852258),
+ trace=TRUE)
4189.237 :   0.008678954 -0.595153967 -4.602990518 -0.625732096
-0.128657978  0.708033556  1.140357461  1.311141424  1.270852258 
Error in numericDeriv(form[[3]], names(ind), env) : 
        Missing value or an Infinity produced when evaluating the model

Excuse me, but I did not know, how I should interpretate it, or what
consequences I have to do.

Kind regards,

Dominik
On 27 Mar 2002, Douglas Bates wrote:

            
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
Dear Dominik,

It has been in my experience that with global variables sometimes it has to
be defined and initialized within the local module before the variables can
use to manipulate data from another source.

Something else came to mind with function g(t), there are times when the
space causes the operation to be read incorrectly, and it see it as, do "g"
then do "t" and not, do function g(t).

When combining modules sometimes you may have to check each module manually
to make sure that it's algorithm routine is functioning correctly. I have
been told some time ago by a wise person, do the mathematical operation
manually then, use a calculator, then, use a spreadsheet and finally enter
it into you code. Here is where you will have additional references to check
you result(s). Give it a try and let me know how you made out...

Regards,
Shawn

-----Original Message-----
From: 1-27206531-0-90000491 [mailto:domi at sun11.ukl.uni-freiburg.de]
Sent: Thursday, March 28, 2002 1:27 PM
To: Shawn Hornsby
Cc: r-help at stat.math.ethz.ch
Subject: RE: [R] Error with nls
On Wed, 27 Mar 2002, Shawn Hornsby wrote:

            
I
that
make
it.
In the origin, it was a compartment model with 6 compartments. Data was
available for every compartment. But of scientific interest are specially
the last 2 compartments, #5 and #6. The function g(t) is the approach to
fill compartment  5 via a forcing function. The parameters b4i and b40 and
k4 results from a fit of compartment 4. In my model of compartment 5 and 6
the parameters b40, b4i and k4 are like time constant covariables.
your
described
nlr
apply
-.
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
Greetings,

Dominik


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._