[R] Beginners Question: Make nlm work
Johannes Graumann
graumann at its.caltech.edu
Thu Aug 26 20:18:12 CEST 2004
Due to the combined nudging of the newsgroup I get it done VERY nicely
now - thanks guys!
The one problem remaining (similar to what Jim explained earlier) is:
how to get to the fitted parameters?
For the code below (were the nls was assigned to 'fit'), 'coef(fit)'
gives me:
> coef(fit)
a b
4.216478980 0.004143935
and 'coef(fit)[1]' results in:
> coef(fit)[1]
a
4.216479
How now to extract just '4.216479'?
Thanks for the help again, everybody!
Joh
Here is the code:
setwd("~/Biology/R_versuch")
mydata<-read.table("YJG45-7_Growth.dat")
x<-mydata$V1
y<-mydata$V2
plot(x, y,
axes=FALSE,
type="p",
xlab="Time (min)",
ylab=expression(OD[600][~nm]),
las=1,
pch=20,
# log="y",
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,las=1,mgp=c(1.5,0.5,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()
fit <- nls(y ~ a/(1+((a-0.008)/0.008)*exp(-(b*x))),
start = list(a = 3, b = 4e-3),
# trace = TRUE
)
yfit1 <- 4.216479/(1+((4.216479-0.008)/0.008)*exp(-(0.004143935*x)))
lines(spline(x, yfit1))
And the data:
# YJG45, YJG46, YJG47 Growth curves 08/12/-08/x/04
# Inoculation with an over night culture to OD600=0.008
# Inoculums:
# YJG45 (A): 2.67 --> use 74.9 myl in 25 ml
# YJG46 (B): 2.65 --> use 75.5 myl in 25 ml
# YJG47 (C): 2.26 --> use 88.5 myl in 25 ml
#
# t OD600
# (min) YJG45 YJG46 YJG47
# A B C
#0 0.008 0.008 0.008
94 0.011 0.013 0.012
188 0.013 0.016 0.018
281 0.018 0.02 0.023
376 0.027 0.026 0.028
454 0.037 0.041 0.039
599 0.078 0.076 0.075
695 0.113 0.121 0.113
814 0.204 0.210 0.197
914 0.3 0.33 0.3
1009 0.47 0.51 0.48
1136 0.75 0.8 0.77
1265 1.18 1.28 1.16
1394 1.78 1.74 1.73
1525 2.13 2.16 2.20
1652 2.64 2.64 2.67
1771 3.08 3.01 2.99
1897 3.38 3.43 3.41
2036 3.78 3.76 3.72
2164 4.1 4.14 4.17
On Wed, 25 Aug 2004 22:39:12 -0400
"Jim Brennan" <jfbrennan at rogers.com> wrote:
> I can't open your file it comes up missing or damaged.
> Maybe you could send me your data as I am curious now to see the
> result. :-) probably something still wrong if it is not fitting the
> points properly...
>
> Jim
> ----- Original Message -----
> From: "Johannes Graumann" <graumann at its.caltech.edu>
> To: "Jim Brennan" <jfbrennan at rogers.com>
> Cc: <r-help at stat.math.ethz.ch>
> Sent: Wednesday, August 25, 2004 9:57 PM
> Subject: Re: [R] Beginners Question: Make nlm work
>
>
> > 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
> > >
> >
>
>
More information about the R-help
mailing list