[R] fitting log function: errors using nls and nlxb

1,582 views
Skip to first unread message

Elizabeth Webb

unread,
Jul 8, 2013, 10:27:24 PM7/8/13
to r-h...@r-project.org
Hi-
I am trying to fit a log function to my data, with the ultimate goal of
finding the second derivative of the function. However, I am stalled on
the first step of fitting a curve.

When I use the following code:
FG2.model<-(nls((CO2~log(a*Time)+b), start=setNames(coef(lm(CO2 ~
log(Time), data=FG2)), c("a", "b")),data=FG2))
I get the following error:
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In log(a * Time) : NaNs produced
4: In log(a * Time) : NaNs produced

When I fit the curve in Plot and use the coefficients as starting values:
start=c(a=68,b=400)
FG2.model<-(nls((CO2~log(a*Time)+b), start=start,data=FG2))
I get the following error:
Error in nls((CO2 ~ log(a * Time) + b), start = start, data = FG2) :
singular gradient
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

So then when I substituded nlxb for nls in the above two models, I got this
error:
Error in nlxb((CO2 ~ log(a * Time) + b), start = start, data = FG2) :
NaN in Jacobian


A few questions:
1.) How can I get R to fit my curve without returning errors?

2.) I am not sure that this data is log base 10. Is there a way I can ask
R to try for logs of different functions? For example,

FG2.model<-(nlxb((CO2~log(a*Time,c)+b), start=start,data=FG2)), where c is
an additional variable. When I try this, R tells me Non-numeric argument
to mathematical function

Thank you in advance,
Elizabeth

[[alternative HTML version deleted]]

______________________________________________
R-h...@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.

Adams, Jean

unread,
Jul 9, 2013, 8:16:20 AM7/9/13
to Elizabeth Webb, R help
Elizabeth,

It's difficult to troubleshoot without the data. Could you provide the
output from
dput(FG2)
or if your data set is quite large, perhaps
dput(FG2[1:50, ])

If you want to fit a third parameter to represent the base of the log, you
could use
nls(CO2 ~ log(a*Time) / log(c) + b, start=c(a=68, b=400, c=10),
data=FG2)
where c represents the base of the log in this relation:
CO2 = log_c(a*Time) + b

Jean




On Mon, Jul 8, 2013 at 9:27 PM, Elizabeth Webb
<webb.eli...@gmail.com>wrote:

Adams, Jean

unread,
Jul 9, 2013, 9:35:30 AM7/9/13
to Elizabeth Webb, R help
Elizabeth,

You should cc rhelp on all correspondence so other readers can follow the
thread of conversation.

Now that I have some data to play with, I see where I went wrong in my
previous e-mail.

First of all, you can't fit a model
CO2 ~ log(a*Time) + b
because log(a*Time) can be rewritten as log(a) + log(Time) and there's no
way that the two parameters, log(a) and b, can be uniquely estimated.

You could change your model to
CO2 ~ a*log(Time) + b
I did this and fit the model using both the base 10 and base e logs ...
same result. So the base of the logs doesn't matter. However, this model
doesn't do a great job of fitting the curve.

I tried another model
CO2 ~ a*(1-exp(-b*Time))
which seemed to do better. Still not great, though. So, I tried it with
one more parameter
CO2 ~ c + a*(1-exp(-b*Time))
and that improved it further.

Jean


fit1 <- nls(CO2 ~ a*log(Time) + b, start=c(a=68, b=400), data=FG2)
plot(FG2$Time, FG2$CO2)
lines(FG2$Time, predict(fit1), col="red")

fit2 <- nls(CO2 ~ a*log10(Time) + b, start=c(a=68, b=400), data=FG2)
lines(FG2$Time, predict(fit1), col="blue", lty=2, lwd=2)

fit3 <- nls(CO2 ~ a*(1-exp(-b*Time)), start=c(a=500, b=0.03), data=FG2)
lines(FG2$Time, predict(fit3), col="green", lwd=2)

fit4 <- nls(CO2 ~ c + a*(1-exp(-b*Time)), start=c(a=500, b=0.03, c=0),
data=FG2)
lines(FG2$Time, predict(fit4), col="brown", lwd=2)





On Tue, Jul 9, 2013 at 8:10 AM, Elizabeth Webb
<webb.eli...@gmail.com>wrote:

> Hi Jean-
> Thanks for responding. Below I have provided dput(FG2[1:50, ]). I have
> also attached a graph of my data.
> Thanks for the hint on working a third parameter into my model. I will
> certainly try that once I get the model working.
> Elizabeth
>
> dput(FG2[1:50, ])
> structure(list(CO2 = c(383.29, 392, 394.38, 392.85, 413.14, 394.56,
> 405.83, 409.61, 408.15, 412.63, 414.62, 423.19, 422.39, 426.81,
> 433.34, 433.95, 438.02, 438.21, 442.84, 441.81, 444.09, 444.59,
> 446.35, 447.11, 450.03, 452.03, 452.69, 453.7, 455.17, 456.65,
> 458.72, 458.88, 459.25, 459.88, 464.06, 461.34, 464.66, 465.19,
> 466.96, 466.86, 468.41, 469.49, 471.08, 471.61, 472.95, 473.94,
> 474.63, 475.79, 477.07, 476.53), Time = c(53, 54, 55, 56, 57,
> 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73,
> 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
> 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102)), .Names = c("CO2",
> "Time"), row.names = c(NA, 50L), class = "data.frame")

Bert Gunter

unread,
Jul 9, 2013, 9:51:45 AM7/9/13
to Adams, Jean, R help, Elizabeth Webb
1. This is a statistical, not an R issue. So I am keeping this offlist.

2. Without a prior family of models, you should **not** be fitting a
parametric nonlinear model.

3. An interpolating smooth is what is wanted. FIt one (e.g. splines,
kernel smooth, etc.)

4. You are out of your depth statistically, Elizabeth. FInd someone
local to work with.

Cheers,
Bert
--

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

Prof J C Nash (U30A)

unread,
Jul 10, 2013, 10:55:42 AM7/10/13
to r-h...@r-project.org, webb.eli...@gmail.com
This reply only addresses the NaN in Jacobian matter. I believe it is a
result of getting a perfect fit (0 sum of squares). I have amended the
r-forge version of nlmrt package in routines nlfb and nlxb and did not
get the error running Elizabeth's example. This only answers the
software issue, of course, not the statistical one.

Use the version of nlmrt from the SCM repository on
https://r-forge.r-project.org/R/?group_id=395

or email me for a tarball of this.

JN
Reply all
Reply to author
Forward
0 new messages