[R] function in nls argument

[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 at r-project.org 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.
>