Skip to content

nls convergence trouble

7 messages · Gabor Grothendieck, Benoit Boulinguiez, PIKAL Petr +1 more

Gabor Grothendieck
# 
Try squaring both sides of the formula.

On Wed, Sep 3, 2008 at 10:01 AM, Benoit Boulinguiez
<benoit.boulinguiez at ensc-rennes.fr> wrote:
PIKAL Petr
# 
Hi

Excel fit is not exceptionally good. Try

 fff<-function(a,b) (V + b * m * a + C0 * V * b - ((C0 * V * b)^2 + 2 * C0 
*
+     b * V^2 - 2 * C0 * V * m * a * b^2 + V^2 + 2 * V * m * a *
+     b + (b * m * a)^2)^(1/2))/(2 * b * m)

and with attached data frame

plot(Qe,fff(364,0.0126))
abline(0,1)

you clearly see linear relationship in smaller values but quite chaotic 
behaviour in bigger ones (or big deviation of experimental points from 
your model).

So it is up to you if you want any fit (like from Excel) or only a good 
one (like from R). 

Seems to me that simple linear could be quite a good choice although there 
is some nelinearity.

fit<-lm(Qe~Ce+C0+V+m)
summary(fit)

Call:
lm(formula = Qe ~ Ce + C0 + V + m)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.654  -8.653   2.426   9.971  11.912 

Coefficients:
              Estimate Std. Error t value Pr(>|t|) 
(Intercept) -8.148e+02  1.330e+03  -0.613 0.549254 
Ce          -6.894e-02  4.982e-03 -13.839 6.02e-10 ***
C0           3.284e-02  1.676e-03  19.589 4.26e-12 ***
V            2.153e+06  4.607e+05   4.674 0.000300 ***
m           -4.272e+04  1.218e+04  -3.509 0.003167 ** 
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

Residual standard error: 10.87 on 15 degrees of freedom
Multiple R-squared: 0.9903,     Adjusted R-squared: 0.9877 
F-statistic: 381.3 on 4 and 15 DF,  p-value: 6.91e-15 

plot(predict(fit), Qe)
abline(0,1)

Regards
Petr


r-help-bounces at r-project.org napsal dne 03.09.2008 16:01:36:
with
nls(formula=LgmAltFormula,data=bois.DATA,start=list(a=300,b=0.01),trace=TRUE
list(a =
the
http://www.R-project.org/posting-guide.html
Benoit Boulinguiez
# 
Hi,

I agree with you that Excel is not the best tool for fittings, that's why I
try to handle R.

But I need to use this specific model ("LgmAltFormula") and not a polynomial
expression with the different parameters even if your method produced
correct fitting.

The parameters "a" and "b" are the Langmuir parameters that describe the
adsorption of a compound onto activated carbon. I need to assess these
parameters.


Regards/Cordialement


Benoit Boulinguiez 


-----Message d'origine-----
De : Petr PIKAL [mailto:petr.pikal at precheza.cz] 
Envoy? : mercredi 3 septembre 2008 17:58
? : Benoit Boulinguiez
Cc : r-help at r-project.org
Objet : Odp: [R] nls convergence trouble

Hi

Excel fit is not exceptionally good. Try

 fff<-function(a,b) (V + b * m * a + C0 * V * b - ((C0 * V * b)^2 + 2 * C0 
*
+     b * V^2 - 2 * C0 * V * m * a * b^2 + V^2 + 2 * V * m * a *
+     b + (b * m * a)^2)^(1/2))/(2 * b * m)

and with attached data frame

plot(Qe,fff(364,0.0126))
abline(0,1)

you clearly see linear relationship in smaller values but quite chaotic 
behaviour in bigger ones (or big deviation of experimental points from 
your model).

So it is up to you if you want any fit (like from Excel) or only a good 
one (like from R). 

Seems to me that simple linear could be quite a good choice although there 
is some nelinearity.

fit<-lm(Qe~Ce+C0+V+m)
summary(fit)

Call:
lm(formula = Qe ~ Ce + C0 + V + m)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.654  -8.653   2.426   9.971  11.912 

Coefficients:
              Estimate Std. Error t value Pr(>|t|) 
(Intercept) -8.148e+02  1.330e+03  -0.613 0.549254 
Ce          -6.894e-02  4.982e-03 -13.839 6.02e-10 ***
C0           3.284e-02  1.676e-03  19.589 4.26e-12 ***
V            2.153e+06  4.607e+05   4.674 0.000300 ***
m           -4.272e+04  1.218e+04  -3.509 0.003167 ** 
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

Residual standard error: 10.87 on 15 degrees of freedom
Multiple R-squared: 0.9903,     Adjusted R-squared: 0.9877 
F-statistic: 381.3 on 4 and 15 DF,  p-value: 6.91e-15 

plot(predict(fit), Qe)
abline(0,1)

Regards
Petr


r-help-bounces at r-project.org napsal dne 03.09.2008 16:01:36:
with
nls(formula=LgmAltFormula,data=bois.DATA,start=list(a=300,b=0.01),trace=TRUE
list(a =
the
http://www.R-project.org/posting-guide.html
PIKAL Petr
# 
Well, the model you try to fit is almost independent of "b"

try

plot(Qe,fff1(364,0.0126))
points(Qe,fff1(364,10000), col=2)

So you can quite freely choose b without any substantial improvement to 
fit. I am not an expert in nonlinear fitting but seems to me that your 
data does not follow your model, so you need to change your model or 
change your data, depends on your preference.

Regards
Petr

r-help-bounces at r-project.org napsal dne 04.09.2008 12:11:27:
why I
polynomial
C0
there
data.
nls(formula=LgmAltFormula,data=bois.DATA,start=list(a=300,b=0.01),trace=TRUE
iterations
Excel
http://www.R-project.org/posting-guide.html
Christian Ritz
# 
Hi Benoit,

another way of making Petr's point is by looking at the profile log likelihood function
for b; that is, only estimating the a parameter for a grid of b values:


## Defining mean function for fixed b
lgma <- function(b){

    function(C0, m, V, a){ (V + b * m * a + C0 * V * b - ((C0 * V * b)^2 + 2 * C0 *
    b * V^2 - 2 * C0 * V * m * a * b^2 + V^2 + 2 * V * m * a *
    b + (b * m * a)^2)^(1/2))/(2 * b * m)}

}

## Defining profile log likelihood function
logLikb <- function(b)
{
    logLik(nls(Qe~(lgma(b))(C0, m, V, a),data = bois.DATA,start = list(a=300)))
}
logLikb2 <- Vectorize(logLikb, "b")  # vectorising the function

## Plotting the profile function
plot(x<-10^(seq(-4, 0, length.out=50)), logLikb2(x), type="l")


Essentially any b value from 0.2 upwards results in the same model fit.

Christian
Gabor Grothendieck
# 
+ e <- expression((V + b * m * a + C0 * V * b - ((C0 * V * b)^2 + 2 *
+ C0 * b * V^2 - 2 * C0 * V * m * a * b^2 + V^2 + 2 * V * m *
+ a * b + (b * m * a)^2))/(2 * b * m))
+ val <- eval(e)
+ attr(val, "gradient") <- cbind(a = eval(D(e, "a")), b = eval(D(e, "b")))
+ val
+ }
Nonlinear regression model
  model:  Qe ~ lgmg(a, b, C0, m, V)
   data:  bois.DATA
        a         b
337.74912   0.03864
 residual sum-of-squares: 15473

Number of iterations to convergence: 9
Achieved convergence tolerance: 3.16e-06
+ (V + b * m * a + C0 * V * b - ((C0 * V * b)^2 + 2 *
+ C0 * b * V^2 - 2 * C0 * V * m * a * b^2 + V^2 + 2 * V * m *
+ a * b + (b * m * a)^2))/(2 * b * m)
+ }
Nonlinear regression model
  model:  Qe^2 ~ lgm(a, b, C0, m, V)^2
   data:  bois.DATA
       a        b
225.6474   0.3568
 residual sum-of-squares: 9.98e+10

Number of iterations to convergence: 16
Achieved convergence tolerance: 6.096e-06
Nonlinear regression model
  model:  Qe ~ lgm(a, 1/b, C0, m, V)
   data:  bois.DATA
     a      b
337.75  25.88
 residual sum-of-squares: 15473

Number of iterations to convergence: 12
Achieved convergence tolerance: 1.722e-06
a          b
1 337.7492 0.03863738
On Wed, Sep 3, 2008 at 10:36 AM, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote: