[R] solving linear equations

Ben Bolker bbolker at gmail.com
Wed Sep 14 03:48:26 CEST 2011


varunshivashankar <varunshivashankar <at> gmail.com> writes:

> 
> I have a dataset 
> 
> X	Y1
> 1200	1.375
> 4000	0.464
> 1333.33	0.148
> 444.44	0.047
> 148.148	0.014
> 49.383	0.005
> 16.461	0.004
> 
> I have to find a curve fit for the above dataset based on a 4-parameter
> logistic equation viz.
> 
> Y1 = d + ((a-d)/(1+(X/cc)^b)), where X and Y1 are the values above.
> 
> I need to know how to solve the above equation for values a, b, c, d.
> 


  1. They're not linear equations.
  2. ?nls
  3. You're going to need very good starting conditions, and a bit
of good luck, to fit a 4-parameter nonlinear equation to 7 data points ...
Looking at a plot of these data, the point at (1200,1.375) is a huge
outlier, which is going to make it almost impossible.
 4. Your equation doesn't look logistic to me, it looks hyperbolic.
Are you sure you have it right?

  This can be done but it's quite delicate:

z <- read.table(textConnection(" X	Y1
1200	1.375
4000	0.464
1333.33	0.148
444.44	0.047
148.148	0.014
49.383	0.005
16.461	0.004"),
header=TRUE)                

with(z,plot(X,Y1))
with(z[-1,],plot(X,Y1))
tmpf <- function(X,a,b,cc,d) d + ((a-d)/(1+(X/cc)^b))
curve(tmpf(x,0.5,-1,1000,0),add=TRUE)

## fix a and d to reasonable values, fit b and cc
n1 <- nls(Y1~tmpf(X,a=0.7,b,cc,d=0),
    data=z[-1,],
    start=list(cc=1000,b=-1))

with(as.list(coef(n1)),curve(tmpf(x,a=0.7,b,cc=1000,d=0),add=TRUE,col=2))

## now use b and cc as starting values, fit a and d
n2 <- nls(Y1~tmpf(X,a,b=-1.721,cc=2738,d),
    data=z[-1,],
    start=list(a=0.7,d=0))

with(as.list(c(coef(n1),coef(n2))),curve(tmpf(x,a,b,cc,d),add=TRUE,col=4))

## now fit all simultaneously
n3 <- nls(Y1~tmpf(X,a,b,cc,d),
    data=z[-1,],
    start=as.list(c(coef(n1),coef(n2))))
coef(n3)
with(as.list(coef(n3)),curve(tmpf(x,a,b,cc,d),add=TRUE,col=5))

## looks OK

  Note that I didn't try to accommodate your outlier. That will
make the problem nearly impossible.



More information about the R-help mailing list