[R] Beginners Question: Make nlm work

[Johannes Graumann]

Aaaahhhh - there was my conceptual problem ... Thanks!
... but 'Oh, God!'! R sucks at fitting (when directed to do so by the
simple minded)! Or am I really that off?
I attach a plot that contains my data as well as two modeled curves:
- the nice one fitted by gnuplot (nonlinear least-squares (NLLS)
Marquardt-Levenberg)
- the nasty one fitted by the code below.

Both algorithms were started with the same values for p[1] and p[2].

Comments?

Is there any way to access the 2 elements of 'out$estimate' from the
program?

Joh

setwd("~/Biology/R_versuch")
mydata<-read.table("YJG45-7_Growth.dat")
#plot(mydata$V1,mydata$V2,xlab="Time
(h)",ylab=expression(OD[600][~nm]),las=1) x<-mydata$V1
y<-mydata$V2
VH <- function(p) sum(
	(y - p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x))))^2) 
plot(x, y,
	axes=FALSE,
	type="p",
	xlab="Time (min)",
	ylab=expression(OD[600][~nm]),
	las=1,
	pch=20,
	tck=0.015,
	mgp=c(1.5,0.25,0),
)
axis(1,
	at=NULL,
	tick=TRUE,
	tck=0.015,
	mgp=c(1.5,0.25,0)
)
axis(2,
	at=NULL,
	tick=TRUE,
	tck=0.015,
	mgp=c(1.5,0.25,0)
)
axis(3,
	at=c(300,600,900,1200,1500,1800,2100),
	labels=c(5,10,15,20,25,30,35),
	tick=TRUE,
	tck=0.015,
	mgp=c(1.5,0.25,0))
axis(4,
	labels=FALSE,
	at=NULL,
	tick=TRUE,
	tck=0.015,
	mgp=c(1.5,0.25,0)
)
box() 
out <- nlm(VH, p = c(3, 4e-3), hessian = TRUE)
#Plot Results from Gnuplot fit
yfit1 <- 4.21632/(1+((4.21632-0.008)/0.008)*exp(-(0.00414403*x)))
lines(spline(x, yfit1))
#Plot Results from 'out'
yfit2 <-
3.000002050/(1+((3.000002050-0.008)/0.008)*exp(-(0.004587924*x)))
lines(spline(x,yfit2))

On Wed, 25 Aug 2004 21:06:24 -0400
"Jim Brennan" <jfbrennan at rogers.com> wrote:

> 
> ----- Original Message -----
> From: "Johannes Graumann" <johannes_graumann at web.de>
> To: <r-help at stat.math.ethz.ch>
> Sent: Wednesday, August 25, 2004 7:07 PM
> Subject: [R] Beginners Question: Make nlm work
> 
> 
> > Hello,
> >
> > I'm new to this and am trying to teach myself some R by plotting
> > biological data. The growth curve in question is supposed to be
> > fitted to the Verhulst equation, which may be transcribed as
> > follows: f(x)=a/(1+((a-0.008)/0.008)*exp(-(b*x)))
> > - for a known population density (0.008) at t(0).
> >
> > I am trying to rework the example from "An Introduction to R" (p.
> > 72) for my case and am failing miserably. Could somebody glance over
> > the code below and nudge me into the right direction - I must have
> > some major conceptual problem which can't be solved by me staring at
> > it ... Since I'm repeating something I have done with gnuplot I know
> > that 3 and 4e-3 as starting values for the fit are appropriate ...
> >
> > Thanks for any hint,
> >
> > Joh
> >
> > setwd("~/Biology/R_versuch")
> > mydata<-read.table("YJG45-7_Growth.dat")
> > x<-mydata$V1
> > y<-mydata$V2
> > VH <- function(p) y ~ p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x)))
> 
> Shouldn't this be
> VH <- function(p) sum((y -
> p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x))))^2) Where you are
> minimizing the squared deviations from what your given equation is
> when you pass it to nlm?
> 
> Maybe if you send some data that would make it more clear.
> 
> Jim
> 
> 
> > plot(x, y, xlab="Time (h)",ylab=expression(OD[600][~nm]),las=1)
> > out <- nlm(VH, p = c(3, 4e-3), hessian = TRUE)
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
> 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot1.eps
Type: application/postscript
Size: 6160 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20040825/acb8a8b9/plot1.eps