Re: [R] function in nls argument

From: Fernando Moyano <nanomoyano_at_yahoo.com>
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 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.
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.

list of date sections of archive