From: Katharine Mullen <kate_at_few.vu.nl>

Date: Thu, 08 May 2008 21:24:22 +0200 (CEST)

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 Thu 08 May 2008 - 19:27:55 GMT

Date: Thu, 08 May 2008 21:24:22 +0200 (CEST)

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_at_yahoo.com>
**> To: r-help_at_r-project.org
**> 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_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.
**> >>
**> >
**> > ______________________________________________
**> > 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.
**> >
**> >
**>
**> --
**> 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_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.
**>
**>
**>
**> ____________________________________________________________________________________
**>
**> [[elided Yahoo spam]]
**>
**> ______________________________________________
**> 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.
**>
*

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 Thu 08 May 2008 - 19:27:55 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.
*