From: elnano <nanomoyano_at_yahoo.com>

Date: Fri, 09 May 2008 04:07:52 -0700 (PDT)

SM <- SM1 [!is.na(SR1)]

SR <- SR1 [!is.na(SR1)]

q <- 0.90

SM <- SM1 [!is.na(SR1)]

SR <- SR1 [!is.na(SR1)]

q <- 0.9

Date: Fri, 09 May 2008 04:07:52 -0700 (PDT)

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.Received on Fri 09 May 2008 - 11:10:00 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 09 May 2008 - 12:30:37 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*