function in nls argument

[Fernando Moyano]

Greetings R users, maybe there is someone who can help
me with this problem:

I define a function "optim.fun" and want as output the
sum of squared errors between predicted and measured
values, as follows:

optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
        {
        predR <-
(a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
        abserr <- abs(ch2no-predR)
        qnum <- quantile(abserr, probs=0.95, na.rm=T)
     
        is.na(abserr) <- (abserr > qnum)
        errsq <- sum(abserr^2, na.rm=T)
        errsq
        }

Then I want to optimize parameters a,b,d and E as to
minimize the function output with the following:

optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
ch2no, a, b, d, E), data=tower,
start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
na.exclude )

were nulldat=0
At this point I get the following error message:

Error in qr.default(.swts * attr(rhs, "gradient")) :
  NA/NaN/Inf in foreign function call (arg 1)
Warning messages:
1: In if (na.rm) x <- x[!is.na(x)] else if
(any(is.na(x))) stop("missing values and NaN's not
allowed if 'na.rm' is FALSE") ... :
  the condition has length > 1 and only the first
element will be used
(this warning is repeated 12 times)

Question: what does the error mean? What am I doing
wrong?
Thanks a bunch.
Nano
Jen, Germany
Max Planck for Biogeochemistry




      ____________________________________________________________________________________

[[elided Yahoo spam]]

______________________________________________
R-help@... 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.

[Katharine Mullen]

The error message means that the gradient (first derivative of residual
vector with respect to the parameter vector) is not possible to work with;
calling the function qr on the gradient multiplied by the square root of
the weight vector .swts (in your case all 1's) fails.

If you want concrete advice it would be helpful to provide the commented,
minimal, self-contained, reproducible code that the posting guide
requests.  what are the values of ST04, SM08b, ch2no, and tower?

General comments:  If your goal is to minimize sum( (observed -
predicted)^2), the function you give nls to minimize (optim.fun in your
case) should return the vector (observed - predicted), not the scalar sum(
(observed - predicted)^2).  You may want to see the nls.lm function in
package minpack.lm for another way to deal with the problem.

On Wed, 7 May 2008, Fernando Moyano wrote:

> Greetings R users, maybe there is someone who can help
> me with this problem:
>
> I define a function "optim.fun" and want as output the
> sum of squared errors between predicted and measured
> values, as follows:
>
> optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
>         {
>         predR <-
> (a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
>         abserr <- abs(ch2no-predR)
>         qnum <- quantile(abserr, probs=0.95, na.rm=T)
>
>         is.na(abserr) <- (abserr > qnum)
>         errsq <- sum(abserr^2, na.rm=T)
>         errsq
>         }
>
> Then I want to optimize parameters a,b,d and E as to
> minimize the function output with the following:
>
> optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
> ch2no, a, b, d, E), data=tower,
> start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
> na.exclude )
>
> were nulldat=0
> At this point I get the following error message:
>
> Error in qr.default(.swts * attr(rhs, "gradient")) :
>   NA/NaN/Inf in foreign function call (arg 1)
> Warning messages:
> 1: In if (na.rm) x <- x[!is.na(x)] else if
> (any(is.na(x))) stop("missing values and NaN's not
> allowed if 'na.rm' is FALSE") ... :
>   the condition has length > 1 and only the first
> element will be used
> (this warning is repeated 12 times)
>
> Question: what does the error mean? What am I doing
> wrong?
> Thanks a bunch.
> Nano
> Jen, Germany
> Max Planck for Biogeochemistry
>
>
>
>
>       ____________________________________________________________________________________
>
> [[elided Yahoo spam]]
>
> ______________________________________________
> R-help@... 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.
>

______________________________________________
R-help@... 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.

[Fernando Moyano]

I've basically solved the problem using the nls.lm function from the minpack.lm (thanks Katharine) with some modifications for ignoring residuals above a given percentile. This is to avoid the strong influence of points which push my modeled vs. measured values away from the 1:1 line.
I based it on the example given for nls.lm. Here it is:

R                               # soil respiration data
ST <- ST [!is.na(R)]     # soil temeprature data. Had to remove na to make nls.lm work
SM <- SM [!is.na(R)]    # soil moisture data
R <- R [!is.na(R)]        
q <- 0.95                    # quantile
   
p <- c("a"=-0.003, "b"=0.13, "c"=0.50, E=400)    # model parameters with estimated values

Rf <- function(ST, SM, a, b, c, E)
    {
    expr <- expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
    eval(expr)
    }
   
    optim.f <- function(p, ST, SM, R, Rfcall)
    {
    res <- R - do.call("Rfcall", c(list(ST = ST, SM = SM), as.list(p)))   # the nls.lm example divides this by sqrt(R), I don't know why. I removed that.
    abserr <- abs(res)
    qnum <- quantile(abserr, probs=q, na.rm=T)    # calculate the "q" quantile of the absolute errors
    res[abserr > qnum]=0                                  # convert residuals above qnum to  0
    res
    }
   
    Rmodel<- nls.lm(par = p, fn = optim.f, Rfcall = Rf, ST = ST, SM = SM, R = R)
   
    summary(Rmodel)

The only error I still get is when using lower q values. A q of around 0.95 or less (depending on the dataset) gives me completely wrong parameter estimates resulting in negative predicted values. Maybe someone has a suggestion here.

Fernando
   

Katharine Mullen wrote:

The error message means that the gradient (first derivative of residual
vector with respect to the parameter vector) is not possible to work with;
calling the function qr on the gradient multiplied by the square root of
the weight vector .swts (in your case all 1's) fails.

If you want concrete advice it would be helpful to provide the commented,
minimal, self-contained, reproducible code that the posting guide
requests.  what are the values of ST04, SM08b, ch2no, and tower?

General comments:  If your goal is to minimize sum( (observed -
predicted)^2), the function you give nls to minimize (optim.fun in your
case) should return the vector (observed - predicted), not the scalar sum(
(observed - predicted)^2).  You may want to see the nls.lm function in
package minpack.lm for another way to deal with the problem.

On Wed, 7 May 2008, Fernando Moyano wrote:

> Greetings R users, maybe there is someone who can help
> me with this problem:
>
> I define a function "optim.fun" and want as output the
> sum of squared errors between predicted and measured
> values, as follows:
>
> optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
>         {
>         predR <-
> (a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
>         abserr <- abs(ch2no-predR)
>         qnum <- quantile(abserr, probs=0.95, na.rm=T)
>
>         is.na(abserr) <- (abserr > qnum)
>         errsq <- sum(abserr^2, na.rm=T)
>         errsq
>         }
>
> Then I want to optimize parameters a,b,d and E as to
> minimize the function output with the following:
>
> optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
> ch2no, a, b, d, E), data=tower,
> start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
> na.exclude )
>
> were nulldat=0
> At this point I get the following error message:
>
> Error in qr.default(.swts * attr(rhs, "gradient")) :
>   NA/NaN/Inf in foreign function call (arg 1)
> Warning messages:
> 1: In if (na.rm) x <- x[!is.na(x)] else if
> (any(is.na(x))) stop("missing values and NaN's not
> allowed if 'na.rm' is FALSE") ... :
>   the condition has length > 1 and only the first
> element will be used
> (this warning is repeated 12 times)
>
> Question: what does the error mean? What am I doing
> wrong?
> Thanks a bunch.
> Nano
> Jen, Germany
> Max Planck for Biogeochemistry
>
>
>
>
>       ____________________________________________________________________________________
>
> [[elided Yahoo spam]]
>
> ______________________________________________
> [e-mail] 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.
>

______________________________________________
[e-mail] 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.

[Zaihra T]

______________________________________________
R-help@... 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.

[Fernando Moyano]

I solved the problem arising from using certain quantile values simply by printing the iterations with the argument nprint. This gave me correct estimates. I have no idea why.

----- Original Message ----
From: elnano <nanomoyano@...>
To: r-help@...
Sent: Thursday, May 8, 2008 5:43:31 PM
Subject: Re: [R] function in nls argument


I've basically solved the problem using the nls.lm function from the
minpack.lm (thanks Katharine) with some modifications for ignoring residuals
above a given percentile. This is to avoid the strong influence of points
which push my modeled vs. measured values away from the 1:1 line.
I based it on the example given for nls.lm. Here it is:

R                               # soil respiration data
ST <- ST [!is.na(R)]     # soil temeprature data. Had to remove na to make
nls.lm work
SM <- SM [!is.na(R)]    # soil moisture data
R <- R [!is.na(R)]        
q <- 0.95                    # quantile
 
p <- c("a"=-0.003, "b"=0.13, "c"=0.50, E=400)    # model parameters with
estimated values

Rf <- function(ST, SM, a, b, c, E)
    {
    expr <-
expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
    eval(expr)
    }
   
    optim.f <- function(p, ST, SM, R, Rfcall)
    {
    res <- R - do.call("Rfcall", c(list(ST = ST, SM = SM), as.list(p)))   #
the nls.lm example divides this by sqrt(R), I don't know why. I removed
that.
    abserr <- abs(res)
    qnum <- quantile(abserr, probs=q, na.rm=T)    # calculate the "q"
quantile of the absolute errors
    res[abserr > qnum]=0                                  # convert
residuals above qnum to  0
    res
    }
   
    Rmodel<- nls.lm(par = p, fn = optim.f, Rfcall = Rf, ST = ST, SM = SM, R
= R)
   
    summary(Rmodel)

The only error I still get is when using lower q values. A q of around 0.95
or less (depending on the dataset) gives me completely wrong parameter
estimates resulting in negative predicted values. Maybe someone has a
suggestion here.

Fernando
   


Katharine Mullen wrote:

>
> The error message means that the gradient (first derivative of residual
> vector with respect to the parameter vector) is not possible to work with;
> calling the function qr on the gradient multiplied by the square root of
> the weight vector .swts (in your case all 1's) fails.
>
> If you want concrete advice it would be helpful to provide the commented,
> minimal, self-contained, reproducible code that the posting guide
> requests.  what are the values of ST04, SM08b, ch2no, and tower?
>
> General comments:  If your goal is to minimize sum( (observed -
> predicted)^2), the function you give nls to minimize (optim.fun in your
> case) should return the vector (observed - predicted), not the scalar sum(
> (observed - predicted)^2).  You may want to see the nls.lm function in
> package minpack.lm for another way to deal with the problem.
>
> On Wed, 7 May 2008, Fernando Moyano wrote:
>
>> Greetings R users, maybe there is someone who can help
>> me with this problem:
>>
>> I define a function "optim.fun" and want as output the
>> sum of squared errors between predicted and measured
>> values, as follows:
>>
>> optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
>>         {
>>         predR <-
>> (a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
>>         abserr <- abs(ch2no-predR)
>>         qnum <- quantile(abserr, probs=0.95, na.rm=T)
>>
>>        is.na(abserr) <- (abserr > qnum)
>>         errsq <- sum(abserr^2, na.rm=T)
>>         errsq
>>         }
>>
>> Then I want to optimize parameters a,b,d and E as to
>> minimize the function output with the following:
>>
>> optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
>> ch2no, a, b, d, E), data=tower,
>> start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
>> na.exclude )
>>
>> were nulldat=0
>> At this point I get the following error message:
>>
>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>>   NA/NaN/Inf in foreign function call (arg 1)
>> Warning messages:
>> 1: In if (na.rm) x <- x[!is.na(x)] else if
>> (any(is.na(x))) stop("missing values and NaN's not
>> allowed if 'na.rm' is FALSE") ... :
>>   the condition has length > 1 and only the first
>> element will be used
>> (this warning is repeated 12 times)
>>
>> Question: what does the error mean? What am I doing
>> wrong?
>> Thanks a bunch.
>> Nano
>> Jen, Germany
>> Max Planck for Biogeochemistry
>>
>>
>>
>>
>>      
>> ____________________________________________________________________________________
>>
>> [[elided Yahoo spam]]
>>
>> ______________________________________________
>> R-help@... 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.
>>
>
> ______________________________________________
> R-help@... 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.
>
>

--
View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17127514.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help@... 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.



      ____________________________________________________________________________________

[[elided Yahoo spam]]

______________________________________________
R-help@... 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.

[Katharine Mullen]

To make your example reproducible you have to provide the data somehow;
I am pretty sure nprint doesn't effect the solution, but if it does this
would be a bug and I would appreciate a reproducible report.

The example in nls.lm is a little complicated in order to show how to use
an analytical expression for the gradient (possible to provide in the
argument jac); since you don't need this, note that your residual function
can be simplified into something along the lines of

p <- list(a=-0.003, b=0.13, c=0.50, E=400)

optim.f <- function(p, ST, SM, R, q) {
 res <- R -(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
 abserr <- abs(res)
 qnum <- quantile(abserr, probs=q, na.rm=T)
 res[abserr > qnum] <- 0
 res
}

Rmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)

On Thu, 8 May 2008, Fernando Moyano wrote:

> I solved the problem arising from using certain quantile values simply
> by printing the iterations with the argument nprint. This gave me
> correct estimates. I have no idea why.
>
> ----- Original Message ----
> From: elnano <nanomoyano@...>
> To: r-help@...
> Sent: Thursday, May 8, 2008 5:43:31 PM
> Subject: Re: [R] function in nls argument
>
>
> I've basically solved the problem using the nls.lm function from the
> minpack.lm (thanks Katharine) with some modifications for ignoring residuals
> above a given percentile. This is to avoid the strong influence of points
> which push my modeled vs. measured values away from the 1:1 line.
> I based it on the example given for nls.lm. Here it is:
>
> R                               # soil respiration data
> ST <- ST [!is.na(R)]     # soil temeprature data. Had to remove na to make
> nls.lm work
> SM <- SM [!is.na(R)]    # soil moisture data
> R <- R [!is.na(R)]
> q <- 0.95                    # quantile
>
> p <- c("a"=-0.003, "b"=0.13, "c"=0.50, E=400)    # model parameters with
> estimated values
>
> Rf <- function(ST, SM, a, b, c, E)
>     {
>     expr <-
> expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
>     eval(expr)
>     }
>
>     optim.f <- function(p, ST, SM, R, Rfcall)
>     {
>     res <- R - do.call("Rfcall", c(list(ST = ST, SM = SM), as.list(p)))   #
> the nls.lm example divides this by sqrt(R), I don't know why. I removed
> that.
>     abserr <- abs(res)
>     qnum <- quantile(abserr, probs=q, na.rm=T)    # calculate the "q"
> quantile of the absolute errors
>     res[abserr > qnum]=0                                  # convert
> residuals above qnum to  0
>     res
>     }
>
>     Rmodel<- nls.lm(par = p, fn = optim.f, Rfcall = Rf, ST = ST, SM = SM, R
> = R)
>
>     summary(Rmodel)
>
> The only error I still get is when using lower q values. A q of around 0.95
> or less (depending on the dataset) gives me completely wrong parameter
> estimates resulting in negative predicted values. Maybe someone has a
> suggestion here.
>
> Fernando
>
>
>
> Katharine Mullen wrote:
> >
> > The error message means that the gradient (first derivative of residual
> > vector with respect to the parameter vector) is not possible to work with;
> > calling the function qr on the gradient multiplied by the square root of
> > the weight vector .swts (in your case all 1's) fails.
> >
> > If you want concrete advice it would be helpful to provide the commented,
> > minimal, self-contained, reproducible code that the posting guide
> > requests.  what are the values of ST04, SM08b, ch2no, and tower?
> >
> > General comments:  If your goal is to minimize sum( (observed -
> > predicted)^2), the function you give nls to minimize (optim.fun in your
> > case) should return the vector (observed - predicted), not the scalar sum(
> > (observed - predicted)^2).  You may want to see the nls.lm function in
> > package minpack.lm for another way to deal with the problem.
> >
> > On Wed, 7 May 2008, Fernando Moyano wrote:
> >
> >> Greetings R users, maybe there is someone who can help
> >> me with this problem:
> >>
> >> I define a function "optim.fun" and want as output the
> >> sum of squared errors between predicted and measured
> >> values, as follows:
> >>
> >> optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
> >>         {
> >>         predR <-
> >> (a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
> >>         abserr <- abs(ch2no-predR)
> >>         qnum <- quantile(abserr, probs=0.95, na.rm=T)
> >>
> >>        is.na(abserr) <- (abserr > qnum)
> >>         errsq <- sum(abserr^2, na.rm=T)
> >>         errsq
> >>         }
> >>
> >> Then I want to optimize parameters a,b,d and E as to
> >> minimize the function output with the following:
> >>
> >> optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
> >> ch2no, a, b, d, E), data=tower,
> >> start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
> >> na.exclude )
> >>
> >> were nulldat=0
> >> At this point I get the following error message:
> >>
> >> Error in qr.default(.swts * attr(rhs, "gradient")) :
> >>   NA/NaN/Inf in foreign function call (arg 1)
> >> Warning messages:
> >> 1: In if (na.rm) x <- x[!is.na(x)] else if
> >> (any(is.na(x))) stop("missing values and NaN's not
> >> allowed if 'na.rm' is FALSE") ... :
> >>   the condition has length > 1 and only the first
> >> element will be used
> >> (this warning is repeated 12 times)
> >>
> >> Question: what does the error mean? What am I doing
> >> wrong?
> >> Thanks a bunch.
> >> Nano
> >> Jen, Germany
> >> Max Planck for Biogeochemistry
> >>
> >>
> >>
> >>
> >>
> >> ____________________________________________________________________________________
> >>
> >> [[elided Yahoo spam]]
> >>
> >> ______________________________________________
> >> R-help@... 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.
> >>
> >
> > ______________________________________________
> > R-help@... 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.
> >
> >
>
> --
> View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17127514.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@... 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.
>
>
>
>       ____________________________________________________________________________________
>
> [[elided Yahoo spam]]
>
> ______________________________________________
> R-help@... 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.
>

______________________________________________
R-help@... 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.

[Fernando Moyano]

Thank you Katharine. I am certain nprint is affecting my solution. Let me know how I can send the data (~300Kb). The script I used it:

ST1 <- ST04
SM1 <- SM08            
SR1 <- SRch2

ST <- ST1 [!is.na(SR1)]
SM <- SM1 [!is.na(SR1)]
SR <- SR1 [!is.na(SR1)]
q <- 0.90
   
p <- c("a"=-0.003, "b"=0.13, "c"=0.50, "E"=400)

SRf <- function(ST, SM, a, b, c, E)
    {
    expr <- expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
    eval(expr)
    }
   
    optim.f <- function(p, ST, SM, SR, SRfcall)
    {
    res <- SR - do.call("SRfcall", c(list(ST = ST, SM = SM), as.list(p)))
    abserr <- abs(res)
    qnum <- quantile(abserr, probs=q, na.rm=TRUE)
    res[abserr > qnum]=0    
    res
    }
   
    SRmodel<- nls.lm(par = p, fn = optim.f, SRfcall = SRf, ST = ST, SM = SM, SR = SR, control = nls.lm.control(nprint=1))
   
    summary(SRmodel)
    SRpred <- do.call(SRf, c(list(ST = ST1, SM = SM1), SRmodel$par))
   
    plot(SR1~SRpred)
    a=c(0,2,4,6)
    b=c(0,2,4,6)
    lines(a,b,col=4)

Changing the nprint argument to 0 (or removing nprint) results in different and completely wrong solutions. The same is true for this equivalent simplified script from your example.

ST1 <- ST04
SM1 <- SM08            
SR1 <- SRch2

ST <- ST1 [!is.na(SR1)]
SM <- SM1 [!is.na(SR1)]
SR <- SR1 [!is.na(SR1)]
q <- 0.9

p <- list(a=-0.003, b=0.13, c=0.50, E=400)

optim.f <- function(p, ST, SM, SR, q) {
 res <- SR - (p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))
 abserr <- abs(res)
 qnum <- quantile(abserr, probs=q, na.rm=TRUE)
 res[abserr > qnum] <- 0
 res
}

SRmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, SR = SR, q = q, control = nls.lm.control(nprint=1))

summary(SRmodel)
SRfun <- function(ST, SM, a, b, c, E) {(a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))}
SRpred <- do.call(SRfun, c(list(ST = ST1, SM = SM1), SRmodel$par))
   
plot(SR1~SRpred)
a=seq(0,6,1)
b=seq(0,6,1)
lines(a,b,col=4)

Here, however, I get the following additional error after summary(SRmodel):
Error in param/se : non-numeric argument to binary operator
although the solution is same as for the first script.

Note that a q of 95 is OK, a q of 90 is not, but a q of 50 is OK again...



To make your example reproducible you have to provide the data somehow;
I am pretty sure nprint doesn't effect the solution, but if it does this
would be a bug and I would appreciate a reproducible report.

The example in nls.lm is a little complicated in order to show how to use
an analytical expression for the gradient (possible to provide in the
argument jac); since you don't need this, note that your residual function
can be simplified into something along the lines of

p <- list(a=-0.003, b=0.13, c=0.50, E=400)

optim.f <- function(p, ST, SM, R, q) {
 res <- R -(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
 abserr <- abs(res)
 qnum <- quantile(abserr, probs=q, na.rm=T)
 res[abserr > qnum] <- 0
 res
}

Rmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)

On Thu, 8 May 2008, Fernando Moyano wrote:

> I solved the problem arising from using certain quantile values simply
> by printing the iterations with the argument nprint. This gave me
> correct estimates. I have no idea why.

[Katharine Mullen]

Thanks for the details - it sounds like a bug.  You can either send me the
data in an email off-list or make it available on-line somewhere, so that
I and other people can download it.

On Fri, 9 May 2008, elnano wrote:

>
> Thank you Katharine. I am certain nprint is affecting my solution. Let me
> know how I can send the data (~300Kb). The script I used it:
>
> ST1 <- ST04
> SM1 <- SM08
> SR1 <- SRch2
>
> ST <- ST1 [!is.na(SR1)]
> SM <- SM1 [!is.na(SR1)]
> SR <- SR1 [!is.na(SR1)]
> q <- 0.90
>
> p <- c("a"=-0.003, "b"=0.13, "c"=0.50, "E"=400)
>
> SRf <- function(ST, SM, a, b, c, E)
>     {
>     expr <-
> expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
>     eval(expr)
>     }
>
>     optim.f <- function(p, ST, SM, SR, SRfcall)
>     {
>     res <- SR - do.call("SRfcall", c(list(ST = ST, SM = SM), as.list(p)))
>     abserr <- abs(res)
>     qnum <- quantile(abserr, probs=q, na.rm=TRUE)
>     res[abserr > qnum]=0
>     res
>     }
>
>     SRmodel<- nls.lm(par = p, fn = optim.f, SRfcall = SRf, ST = ST, SM = SM,
> SR = SR, control = nls.lm.control(nprint=1))
>
>     summary(SRmodel)
>     SRpred <- do.call(SRf, c(list(ST = ST1, SM = SM1), SRmodel$par))
>
>     plot(SR1~SRpred)
>     a=c(0,2,4,6)
>     b=c(0,2,4,6)
>     lines(a,b,col=4)
>
> Changing the nprint argument to 0 (or removing nprint) results in different
> and completely wrong solutions. The same is true for this equivalent
> simplified script from your example.
>
> ST1 <- ST04
> SM1 <- SM08
> SR1 <- SRch2
>
> ST <- ST1 [!is.na(SR1)]
> SM <- SM1 [!is.na(SR1)]
> SR <- SR1 [!is.na(SR1)]
> q <- 0.9
>
> p <- list(a=-0.003, b=0.13, c=0.50, E=400)
>
> optim.f <- function(p, ST, SM, SR, q) {
>  res <- SR -
> (p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))
>  abserr <- abs(res)
>  qnum <- quantile(abserr, probs=q, na.rm=TRUE)
>  res[abserr > qnum] <- 0
>  res
> }
>
> SRmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, SR = SR, q = q,
> control = nls.lm.control(nprint=1))
>
> summary(SRmodel)
> SRfun <- function(ST, SM, a, b, c, E)
> {(a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))}
> SRpred <- do.call(SRfun, c(list(ST = ST1, SM = SM1), SRmodel$par))
>
> plot(SR1~SRpred)
> a=seq(0,6,1)
> b=seq(0,6,1)
> lines(a,b,col=4)
>
> Here, however, I get the following additional error after summary(SRmodel):
> Error in param/se : non-numeric argument to binary operator
> although the solution is same as for the first script.
>
> Note that a q of 95 is OK, a q of 90 is not, but a q of 50 is OK again...
>
>
>
> To make your example reproducible you have to provide the data somehow;
> I am pretty sure nprint doesn't effect the solution, but if it does this
> would be a bug and I would appreciate a reproducible report.
>
> The example in nls.lm is a little complicated in order to show how to use
> an analytical expression for the gradient (possible to provide in the
> argument jac); since you don't need this, note that your residual function
> can be simplified into something along the lines of
>
> p <- list(a=-0.003, b=0.13, c=0.50, E=400)
>
> optim.f <- function(p, ST, SM, R, q) {
>  res <- R
> -(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
>  abserr <- abs(res)
>  qnum <- quantile(abserr, probs=q, na.rm=T)
>  res[abserr > qnum] <- 0
>  res
> }
>
> Rmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)
>
> On Thu, 8 May 2008, Fernando Moyano wrote:
>
> > I solved the problem arising from using certain quantile values simply
> > by printing the iterations with the argument nprint. This gave me
> > correct estimates. I have no idea why.
>
> --
> View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17145801.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@... 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.
>

______________________________________________
R-help@... 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.

[Fernando Moyano]

Find the data (data_nls.lm_moyano.txt) here:
ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano

Katharine Mullen wrote:

Thanks for the details - it sounds like a bug.  You can either send me the
data in an email off-list or make it available on-line somewhere, so that
I and other people can download it.


______________________________________________
[e-mail] 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.

[Katharine Mullen]

You indeed found a bug.  I can reproduce it (which I should have tried to
do on other examples in the first place!).  Thanks for finding it.

It will be fixed in version 1.1-0 which I will submit to CRAN soon.

On Fri, 9 May 2008, elnano wrote:

>
> Find the data (data_nls.lm_moyano.txt) here:
> ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano
>
>
>
> Katharine Mullen wrote:
> >
> > Thanks for the details - it sounds like a bug.  You can either send me the
> > data in an email off-list or make it available on-line somewhere, so that
> > I and other people can download it.
> >
> >
> > ______________________________________________
> > R-help@... 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.
> >
> >
>
> --
> View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17146812.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@... 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.
>

______________________________________________
R-help@... 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.

[Katharine Mullen]

You can take minpack.lm_1.1-0 (source code and MS Windows build,
respectively) from here:

http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.tar.gz
http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.zip

The bug that occurs when nprint = 0 is fixed.  Also fixed is another
problem suggested your example: when the argument par is a list, calling
summary on the output of nls.lm was not working.

I'll submit the new version to CRAN soon.

This disscusion has been fruitful - thanks for it.

On Fri, 9 May 2008, Katharine Mullen wrote:

> You indeed found a bug.  I can reproduce it (which I should have tried to
> do on other examples in the first place!).  Thanks for finding it.
>
> It will be fixed in version 1.1-0 which I will submit to CRAN soon.
>
> On Fri, 9 May 2008, elnano wrote:
>
> >
> > Find the data (data_nls.lm_moyano.txt) here:
> > ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano
> >
> >
> >
> > Katharine Mullen wrote:
> > >
> > > Thanks for the details - it sounds like a bug.  You can either send me the
> > > data in an email off-list or make it available on-line somewhere, so that
> > > I and other people can download it.
> > >
> > >
> > > ______________________________________________
> > > R-help@... 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.
> > >
> > >
> >
> > --
> > View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17146812.html
> > Sent from the R help mailing list archive at Nabble.com.
> >
> > ______________________________________________
> > R-help@... 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.
> >
>
> ______________________________________________
> R-help@... 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.
>

______________________________________________
R-help@... 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.

[Martin Maechler]

Hi Kate and Fernando,

I'm late into this thread,
but from reading it I get the impression that Fernando really
wants to do *robust* (as opposed to least-squares) non-linear
model fitting.  His proposal to set residuals to zero when they
are outside a given bound is a very special case of an
M-estimator, namely (if I'm not mistaken) the so-called "Huber
skipped-mean", an M-estimator with psi-function
   psi <- function(x, k) ifelse(abs(x) <= k, x, 0)
It is known that this can be far from optimal, and either using
Huber-psi or "a redescender" such as Tukey's biweight can be
considerably better.
Also note that the standard inference (std.errors, P-values, ...)
that you'd get from summary(nlsfit) or anova(nls1, nl2) is
*invalid* here, since you are effectively using *random* weighting.

The nlrob() function in package 'robustbase'
implements M-estimation of nonlinear models directly.
Unfortunately, how to do correct inference in this situation
is a hard problem, probably even an open research question in
parts. I would expect that "the" bootstrap should work if you only
have a few outliers.

I don't have time at the moment to look at the example data and
the model, and show you how to use it for nlrob();
if you find a way to you it for nls() , then the same should
work for nlrob().

I'm CCing this to the specialists for "Robust Stats with R"
mailing list, R-SIG-robust.

Best regards,
Martin Maechler
ETH Zurich

>>>>> "KateM" == Katharine Mullen <kate@...>
>>>>>     on Fri, 9 May 2008 15:50:08 +0200 (CEST) writes:

    KateM> You can take minpack.lm_1.1-0 (source code and MS Windows build,
    KateM> respectively) from here:

    KateM> http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.tar.gz
    KateM> http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.zip

    KateM> The bug that occurs when nprint = 0 is fixed.  Also fixed is another
    KateM> problem suggested your example: when the argument par is a list, calling
    KateM> summary on the output of nls.lm was not working.

    KateM> I'll submit the new version to CRAN soon.

    KateM> This disscusion has been fruitful - thanks for it.

    KateM> On Fri, 9 May 2008, Katharine Mullen wrote:

    >> You indeed found a bug.  I can reproduce it (which I should have tried to
    >> do on other examples in the first place!).  Thanks for finding it.
    >>
    >> It will be fixed in version 1.1-0 which I will submit to CRAN soon.
    >>
    >> On Fri, 9 May 2008, elnano wrote:
    >>
    >> >
    >> > Find the data (data_nls.lm_moyano.txt) here:
    >> > ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano
    >> >
    >> >
    >> >
    >> > Katharine Mullen wrote:
    >> > >
    >> > > Thanks for the details - it sounds like a bug.  You can either send me the
    >> > > data in an email off-list or make it available on-line somewhere, so that
    >> > > I and other people can download it.
    >> > >
    >> > >
    >> > > ______________________________________________
    >> > > R-help@... 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.
    >> > >
    >> > >
    >> >
    >> > --
    >> > View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17146812.html
    >> > Sent from the R help mailing list archive at Nabble.com.
    >> >
    >> > ______________________________________________
    >> > R-help@... 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.
    >> >
    >>
    >> ______________________________________________
    >> R-help@... 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.
    >>

    KateM> ______________________________________________
    KateM> R-help@... mailing list
    KateM> https://stat.ethz.ch/mailman/listinfo/r-help
    KateM> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
    KateM> and provide commented, minimal, self-contained, reproducible code.

______________________________________________
R-help@... 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.