|
»
»
« Return to Thread: function in nls argument
Re: function in nls argument
by Fernando Moyano
::
Rate this Message:
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
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@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.
« Return to Thread: function in nls argument
| Free Forum Powered by Nabble | Forum Help |