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