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

Tue Jul 9 15:51:45 CEST 2013

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

On Tue, Jul 9, 2013 at 6:35 AM, Adams, Jean <jvadams at usgs.gov> wrote:
> 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.elizabeth.e at 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")
>>
>>
>> On Tue, Jul 9, 2013 at 4:16 AM, Adams, Jean <jvadams at usgs.gov> wrote:
>>
>>> 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.elizabeth.e at gmail.com> wrote:
>>>
>>>> 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-help at 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.
>>>>
>>>
>>>
>>
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.

-- 

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