From: Fernando Moyano <nanomoyano_at_yahoo.com>

Date: Thu, 08 May 2008 10:07:52 -0700 (PDT)

SM <- SM [!is.na(R)] # soil moisture data

>> 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.
*

*>>
*

Date: Thu, 08 May 2008 10:07:52 -0700 (PDT)

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.

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 makenls.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 # convertresiduals 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

> > ______________________________________________ > 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.Received on Thu 08 May 2008 - 17:40: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 Thu 08 May 2008 - 20:31:17 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.
*