# Re: [R] nonlinear fitting on many voxels

From: Charles C. Berry <cberry_at_tajo.ucsd.edu>
Date: Mon, 09 Jun 2008 19:43:29 -0700

On Mon, 9 Jun 2008, Ayman Oweida wrote:

>
> Below is a sample code using random generated numbers that represents
> what I'm trying to do. I ran it a few times until I got the same
> error.&nbsp; I hope this can help in defining where I'm going
> wrong.&nbsp;
>

The first item under 'Technical Details' in the posting guide says

'No HTML...'

This includes stuff like "&lt;-" (ampersand - ell - tee - ; - minus).

You get the error because you have NaN showing up in the parameter vector.

My suggestion to set options( error = recover ) and try to trace this through yourself still holds.

Looking at your code I wonder if 'p' needs to obey some constraints. Does p > 0 need to hold?? p[1] < p[2]??

Perhaps a reparameterization or use of boundary constraints would straighten things out.

HTH, Chuck

> library(nlme)
>
> Ref&lt;-runif(10) # data for reference region
> R1Tn&lt;-runif(10) # data for region of interest
> R10T&lt;-rep(0.1,10) # initial value for each point
> t&lt;- c(0:9) #time
> KR&lt;-0.1 #constant
> VR&lt;-0.1 #constant
>
> f.int&lt;-array(0,length(t)) #array to contain integral values of function f
>
> f&lt;-function(p){
> for (i in 1:(length(Ref)-1)){
>
> f.int[i+1] &lt;- sintegral(t[1:(i+1)], (Ref[1:(i+1)] * (exp((-p[1]/p[2])*(t[i+1]-t)))[1:(i+1)]))
> }
>
> F &lt;- (p[1]/KR)*(Ref) + (p[1]/KR)*( (KR/VR)-(p[1]/p[2]) )* f.int
> F
> }
>
> cf&lt;-array(0,length(t)) #array to contain the first estimated parameter for each voxel
> for (j in 1:length(t)){
> res&lt;- function(p) sum(((R1Tn[j]-R10T[j])-f(p))^2)&nbsp;&nbsp;
>
> out &lt;-&nbsp; ( nlminb(c(0.1,0.25),res) )
>
> cf[j]&lt;-out\$par[1] # fill the 1st estimated parameter for each voxel into cf
> }
>
> Thank you,
> Ayman
> --- On Mon, 6/9/08, Charles C. Berry &lt;cberry_at_tajo.ucsd.edu&gt; wrote:
> From: Charles C. Berry &lt;cberry_at_tajo.ucsd.edu&gt;
> Subject: Re: [R] nonlinear fitting on many voxels
> To: "Ayman Oweida" &lt;ajoweida_at_yahoo.com&gt;
> Cc: r-help_at_r-project.org
> Date: Monday, June 9, 2008, 12:59 PM
>
> On Mon, 9 Jun 2008, Ayman Oweida wrote:
>
> &gt; After many months, I am now banging my head against the wall because I
> can't find a solution to this seemingly trivial problem.&amp;nbsp; Any help
> would be appreciated:
> &gt;
> &gt; I am trying to apply a nonlinear fitting routine to a 3D MR image on a
> voxel-by-voxel basis.&amp;nbsp; I've tested the routine using simulated
> data and things went well.&amp;nbsp; As for the real data, the fitting routine
> works variably.&amp;nbsp; By variably, I mean the following: &amp;nbsp; with a
> specific set of starting parameters the routine would work for the first 10
> voxels, for another set of starting values the routine would work for the first
> 1000 voxels and so on.&amp;nbsp; NEVER was I able to get starting values that
> would allow the routine to run entirely!&amp;nbsp; The error I would get after
> fitting the limited number of voxels is:&amp;nbsp;
> &gt;
> &gt; "Error in approx(x, fx, n = 2 * n.pts + 1) :
> &gt; &amp;nbsp;&amp;nbsp;&amp;nbsp;&amp;nbsp;&amp;nbsp;&amp;nbsp;&amp;nbsp;
> need at least two non-NA values to interpolate"
> &gt;
> &gt; I think the error is from the library Bolstad since I'm using the
> &gt; sintegral function as part of the fitting equation.&amp;nbsp; I've
> tried
> &gt; numerous nonlinear functions including: nlm, optim and nlminb, but all
> &gt; stop after performing the fit on a limited number of voxels.&amp;nbsp; I
> &gt; took the values for which the fitting routine stops and I applied
> &gt; different starting values to them and it works!&amp;nbsp; I'm still
> using
> &gt; the same voxel/parameter values and I am sure they are not NA values as
> &gt; per the error!
>
> Reread the error message. It does NOT say that there are NA values. It
> only says that you 'need at least two non-NA values'.
>
> &gt;
> &gt; Ok, I think the problem is clear.&amp;nbsp; any solutions?
>
> The problem is not clear.
>
>
> http://www.R-project.org/posting-guide.html and provide commented,
> minimal, self-contained, reproducible code.
>
> Often, developing the 'minimal' example helps you to perceive the
> underlying difficulty. Without it, you will only get very general advice
> and off-hand guesses about where your problem lies, which may not move you
> closer to a solution.
>
> ---
>
> So here is some of that general advice: find out why approx() sends you
> that error message.
>
> One way to do this is to set
>
> options(error=recover)
>
> before running your function. And then inspect objects in the frame in
> which the error occurred and in the frames leading up to the one in which
> the error was triggered. See
>
> ?recover
>
> HTH,
>
> Chuck
>
>
> &amp;nbsp; I am new to this complex world of statistical analysis.
> &gt; Your help is very very much appreciated.
> &gt;
> &gt; Thanks,
> &gt; Ayman
> &gt;
> &gt; Montreal Neurological Institute
> &gt;
> &gt;
> &gt;
> &gt;
> &gt; [[alternative HTML version deleted]]
> &gt;
> &gt; ______________________________________________
> &gt; R-help_at_r-project.org mailing list
> &gt; https://stat.ethz.ch/mailman/listinfo/r-help
> http://www.R-project.org/posting-guide.html
> &gt; and provide commented, minimal, self-contained, reproducible code.
> &gt;
>
> Charles C. Berry (858) 534-2098
> Dept of Family/Preventive Medicine
> E mailto:cberry_at_tajo.ucsd.edu UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.
>

```Charles C. Berry                            (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cberry_at_tajo.ucsd.edu	            UC San Diego
```
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901

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. Received on Tue 10 Jun 2008 - 03:25:23 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 10 Jun 2008 - 03:30:43 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.