Re: [R] non-linear regression and root finding
On Mon, 6 Nov 2023 14:55:39 -0500 J C Nash wrote: > However, I'm wondering if this approach is worth writing up, at least > as a vignette or blog post. It does need a shorter example and some > explanation of the "why" and some testing perhaps. Do you mean using this problem as a basis to illustrate ordering constraints on parameters? Weird constraints do come up every now and then in regression problems. I could definitely offer my help with at least some of the text. > If there's interest, I'll be happy to join in. And my own posting > suggests how the ordering is enforced by bounding the "delta" > parameters from below. I have just tried nlsr::nlxb for a slightly larger dataset shared by Troels off-list, and it worked great with the delta parameters as you suggested, thank you! It's interesting that nlxb and nlsLM give slightly different answers, differing in 0.5 pK units for pK1 and (pK2-pK1) but not (pK3-pK2). Then again, they both agree that the standard error for pK1 and (pK2-pK1) is very large, so perhaps the problem is just very ill-conditioned. -- Best regards, Ivan __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
G'day Troels, On Tue, 7 Nov 2023 07:14:02 +0100 Troels Ring wrote: > Be as it may, I wonder if not your method might work if only we KNOW > that pK1 is either positive OR negative, in which case we have pK1 = > -exp(theta1)? If pK1 can be either negative or positive (or 0 :-) ), and it is just the ordering that you want to have enforced, then I would try the parameterisation: pK1 = pK1 pK2 = pK1 + exp(theta2) pK3 = pk2 + exp(theta3) and optimise over pK1, theta2 and theta3. As long as you want to know the estimates only. Asking for standard errors of the original estimates would open another can of worms. :-) Cheers, Berwin __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
Thanks a lot, Berwin. Unfortunately, pK1 may well be negative and as I understand the literature it may be poorly defined as such, and also seems to be at a boundary, since when lower is set to say rep(-4,3) pK1 is returned as -4 while pK2 and pK3 are undisturbed. Perhaps the point is that pK1 is not carrying any information at the pH around 5. Fair enough, I guess. Only, I believe I need stick to including all three pK values to be in agreement with the molecular information about HEPES - although even this is contentious. Be as it may, I wonder if not your method might work if only we KNOW that pK1 is either positive OR negative, in which case we have pK1 = -exp(theta1)? Best wishes Troels Den 07-11-2023 kl. 05:10 skrev Berwin A Turlach: G'day Troels, On Mon, 6 Nov 2023 20:43:10 +0100 Troels Ring wrote: Thanks a lot! This was amazing. I'm not sure I see how the conditiion pK1 < pK2 < pK3 is enforced? One way of enforcing such constraints (well, in finite computer arithemtic only "<=" can be enforced) is to rewrite the parameters as: pK1 = exp(theta1) ## only if pK1 > 0 pK2 = pK1 + exp(theta2) pK3 = pk2 + exp(theta3) And then use your optimiser to optimise over theta1, theta2 and theta3. There might be better approaches depending on the specific problem. Cheers, Berwin __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
G'day Troels, On Mon, 6 Nov 2023 20:43:10 +0100 Troels Ring wrote: > Thanks a lot! This was amazing. I'm not sure I see how the conditiion > pK1 < pK2 < pK3 is enforced? One way of enforcing such constraints (well, in finite computer arithemtic only "<=" can be enforced) is to rewrite the parameters as: pK1 = exp(theta1) ## only if pK1 > 0 pK2 = pK1 + exp(theta2) pK3 = pk2 + exp(theta3) And then use your optimiser to optimise over theta1, theta2 and theta3. There might be better approaches depending on the specific problem. Cheers, Berwin __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
I won't send to list, but just to the two of you, as I don't have anything to add at this time. However, I'm wondering if this approach is worth writing up, at least as a vignette or blog post. It does need a shorter example and some explanation of the "why" and some testing perhaps. If there's interest, I'll be happy to join in. And my own posting suggests how the ordering is enforced by bounding the "delta" parameters from below. Note that nls() can only handle bounds in the "port" algorithm, and the man page rather pours cold water on using that algorithm. Best, JN On 2023-11-06 14:43, Troels Ring wrote: Thanks a lot! This was amazing. I'm not sure I see how the conditiion pK1 < pK2 < pK3 is enforced? - it comes from the derivation via generalized Henderson-Hasselbalch but perhaps it is not really necessary. Anyway, the use of Vectorize did the trick! Best wishes Troels Den 06-11-2023 kl. 19:19 skrev Ivan Krylov: В Mon, 6 Nov 2023 17:53:49 +0100 Troels Ring пишет: Hence I wonder if I could somehow have non linear regression to find the 3 pK values. Below is HEPESFUNC which delivers charge in the fluid for known pKs, HEPTOT and SID. Is it possible to have root-finding in the formula with nls? Sure. Just reformulate the problem in terms of a function that takes a vector of predictors (your independent variable SID) and the desired parameters (pK1, pK2, pK3) as separate arguments and then returns predicted values of the dependent variable (to compare against pHobs): kw <- 1e-14 # I'm assuming pHm <- Vectorize(function(SID, pK1, pK2, pK3) -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4, HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root)) (Yes, Vectorize() doesn't make the function any faster, but I don't see an easy way to rewrite this function to make it truly vectorised.) Unfortunately, nls() seems to take a step somewhere where crossprod() of the Jacobian of the model function cannot be inverted and fails with the message "Singular gradient". I wish that R could have a more reliable built-in nonlinear least squares solver. (I could also be holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and minpack.lm: minpack.lm::nlsLM( pHobs ~ pHm(SID, pK1, pK2, pK3), data.frame(pHobs = pHobs, SID = SID), start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3), # the following is also needed to avoid MINPACK failing to fit lower = rep(-1, 3), upper = rep(9, 3) ) # Nonlinear regression model # model: pHobs ~ pHm(SID, pK1, pK2, pK3) # data: data.frame(pHobs = pHobs, SID = SID) # pK1 pK2 pK3 # -1.000 2.966 7.606 # residual sum-of-squares: 0.001195 # # Number of iterations to convergence: 15 # Achieved convergence tolerance: 1.49e-08 (Unfortunately, your code seemed to have a lot of non-breakable spaces, which confuse R's parser and make it harder to copy&paste it.) I think that your function can also be presented as a degree-5 polynomial in H, so it should also be possible to use polyroot() to obtain your solutions in a more exact manner. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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 -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
Thanks a lot! This was amazing. I'm not sure I see how the conditiion pK1 < pK2 < pK3 is enforced? - it comes from the derivation via generalized Henderson-Hasselbalch but perhaps it is not really necessary. Anyway, the use of Vectorize did the trick! Best wishes Troels Den 06-11-2023 kl. 19:19 skrev Ivan Krylov: В Mon, 6 Nov 2023 17:53:49 +0100 Troels Ring пишет: Hence I wonder if I could somehow have non linear regression to find the 3 pK values. Below is HEPESFUNC which delivers charge in the fluid for known pKs, HEPTOT and SID. Is it possible to have root-finding in the formula with nls? Sure. Just reformulate the problem in terms of a function that takes a vector of predictors (your independent variable SID) and the desired parameters (pK1, pK2, pK3) as separate arguments and then returns predicted values of the dependent variable (to compare against pHobs): kw <- 1e-14 # I'm assuming pHm <- Vectorize(function(SID, pK1, pK2, pK3) -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4, HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root)) (Yes, Vectorize() doesn't make the function any faster, but I don't see an easy way to rewrite this function to make it truly vectorised.) Unfortunately, nls() seems to take a step somewhere where crossprod() of the Jacobian of the model function cannot be inverted and fails with the message "Singular gradient". I wish that R could have a more reliable built-in nonlinear least squares solver. (I could also be holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and minpack.lm: minpack.lm::nlsLM( pHobs ~ pHm(SID, pK1, pK2, pK3), data.frame(pHobs = pHobs, SID = SID), start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3), # the following is also needed to avoid MINPACK failing to fit lower = rep(-1, 3), upper = rep(9, 3) ) # Nonlinear regression model # model: pHobs ~ pHm(SID, pK1, pK2, pK3) #data: data.frame(pHobs = pHobs, SID = SID) #pK1 pK2pK3 # -1.000 2.966 7.606 # residual sum-of-squares: 0.001195 # # Number of iterations to convergence: 15 # Achieved convergence tolerance: 1.49e-08 (Unfortunately, your code seemed to have a lot of non-breakable spaces, which confuse R's parser and make it harder to copy&paste it.) I think that your function can also be presented as a degree-5 polynomial in H, so it should also be possible to use polyroot() to obtain your solutions in a more exact manner. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] non-linear regression and root finding
Your script is missing something (in particular kw).
I presume you are trying to estimate the pK values. You may have more success
with package nlsr than nls(). nlsr::nlxb() tries to get the Jacobian of the
model specified by a formula and do so by applying symbolic or automatic
differentiation. The multi expression function would probably not work.
(That is one advantage of nls(), but it has an unstabilized solver and unless
carefully called will use a simple forward derivative approximation.)
There is also nlsr::nlfb() that can use functions for the residuals and
jacobian. Your function is messy, but likely the jacobian can be developed
with a bit of work.
Some of the symbolic derivative features of R can help here. Alternatively,
there are several numerical approximations in nlsr and the vignette "Intro
to nlsr" explains how to use these. Note that simple forward and backward
approximations are generally not worth using, but central approximation is
quite good.
The solver in the nlsr package allows of bounds on the parameters. Since you
want them ordered, you may want to transform
pK1 < pK2 <- pK3
to pk1, deltapk2, deltapk3 so pK2 = pk1+deltapk2
and pk3 = pk2 + deltapk3 = pk1 + deltapk2 + deltapk3
and all values bounded below by 0 or possibly some small numbers to
keep the paramters apart.
I won't pretend any of this is trivial. It is VERY easy to make small errors
that still allow for output that is disastrously incorrect. If it is possible
to get a single "formula" as one expression even if spread over multiple lines,
then nlxb() might be able to handle it.
J Nash (maintainer of nlsr and optimx)
On 2023-11-06 11:53, Troels Ring wrote:
HEPESFUNC <-
function(H,SID,HEPTOT,pK1,pK2,pK3) {
XX <- (H^3/(10^-pK3*10^-pK2*10^-pK1)+H^2/(10^-pK3*10^-pK2)+H/(10^-pK3))
IV <- HEPTOT/(1+XX)
I <- IV*H^3/(10^-pK3*10^-pK2*10^-pK1)
II <- IV*H^2/(10^-pK3*10^-pK2)
III <- IV*H/10^-pK3
H - kw/H + SID + I*2 + II - IV
}
HEPTOT <- 0.050
SID <- c(-seq(10,1)*1e-4,0,seq(1,10)*1e-4)
pHobs <- c(4.63,4.68,4.72,4.77,4.83,4.9,4.96,5.04,5.12,5.21,5.3,
5.39,5.48,5.55,5.63,5.69,5.74,5.8,5.85,5.89,5.93)
pK1 <- -1; pK2 <- 3; pK3 <- 7.55 # literature values
pK3 <- 7.6; pK2 <- 2.96 # values eye-balled to be better
pH <- c()
for (i in 1:length(SID)) {
pH[i] <- -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4,
HEPTOT=HEPTOT,SID = SID[i],
pK1=pK1,pK2=pK2,pK3=pK3)$root)
}
plot(SID,pH)
points(SID,pHobs,col="red")
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
Re: [R] non-linear regression and root finding
В Mon, 6 Nov 2023 17:53:49 +0100 Troels Ring пишет: > Hence I wonder if I could somehow have non linear regression to find > the 3 pK values. Below is HEPESFUNC which delivers charge in the > fluid for known pKs, HEPTOT and SID. Is it possible to have > root-finding in the formula with nls? Sure. Just reformulate the problem in terms of a function that takes a vector of predictors (your independent variable SID) and the desired parameters (pK1, pK2, pK3) as separate arguments and then returns predicted values of the dependent variable (to compare against pHobs): kw <- 1e-14 # I'm assuming pHm <- Vectorize(function(SID, pK1, pK2, pK3) -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4, HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root)) (Yes, Vectorize() doesn't make the function any faster, but I don't see an easy way to rewrite this function to make it truly vectorised.) Unfortunately, nls() seems to take a step somewhere where crossprod() of the Jacobian of the model function cannot be inverted and fails with the message "Singular gradient". I wish that R could have a more reliable built-in nonlinear least squares solver. (I could also be holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and minpack.lm: minpack.lm::nlsLM( pHobs ~ pHm(SID, pK1, pK2, pK3), data.frame(pHobs = pHobs, SID = SID), start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3), # the following is also needed to avoid MINPACK failing to fit lower = rep(-1, 3), upper = rep(9, 3) ) # Nonlinear regression model # model: pHobs ~ pHm(SID, pK1, pK2, pK3) #data: data.frame(pHobs = pHobs, SID = SID) #pK1 pK2pK3 # -1.000 2.966 7.606 # residual sum-of-squares: 0.001195 # # Number of iterations to convergence: 15 # Achieved convergence tolerance: 1.49e-08 (Unfortunately, your code seemed to have a lot of non-breakable spaces, which confuse R's parser and make it harder to copy&paste it.) I think that your function can also be presented as a degree-5 polynomial in H, so it should also be possible to use polyroot() to obtain your solutions in a more exact manner. -- Best regards, Ivan __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Non-Linear Regression Help
I have a fair bit of experience with both nls and rating curves. This is not a nls() problem, this is a model problem. The power law rating curve favored by hydrologists would not apply to your data as it's based on the idea that a log-log plot of discharge vs. stage, or state+a in your case is a straight line, statistical assumptions notwithstanding. A log-log plot of your data, plot(discharge~stage,data=yourdata,pch=19,log='xy') clear is not a straight line. The very large discharge at stage=6.53 vs. stage=6.32 says one of two things: 1) there is an error in the data (perhaps the 2592.05 should be 592.05) or 2) the river channel geometry has changed dramatically, as in overtopping its banks or perhaps there's a smaller central channel set into a larger flood channel, similar to the LA river of the movies. The way forward is 1) recheck your data or 2) recheck your data and use a two-piece model with one piece for stage <= 6.32 and a second piece for stage > 6.32. For this second approach to work, you'll need more data than you have given us here. BTW, nls() should work fine if the model/data combination meet the requirements of 1) the model 'fits' the data, 2) the residuals are NIID(0,sigma^2), the parameters C, a, and n are identifiable from the data (should be if the last point is excluded). As always, you'll need good starting values for the parameters (get them from a log-log plot). You may find, based on the residuals, that linear regression (lm, glm) are most appropriate so that the errors meet the criteria of constant variance. If none of this makes sense, buy and study the book Nonlinear regression analysis: Its applications, D. M. Bates and D. G. Watts, Wiley, New York, 1988. ISBN 0471-816434. The nls() application is the easy part. Good luck David Stevens On 5/4/2017 4:58 PM, Zachary Shadomy wrote: > I am having some errors come up in my first section of code. I have no > issue in plotting the points. Is there an easier method for creating a > non-linear regression using C*(x+a)^n. The .txt file is named > stage_discharge with the two variables being stage and discharge. > The data is a relatively small file listed below: > > stage discharge > 6.53 2592.05 > 6.32 559.5782 > 5.96 484.2151 > 4.99 494.7527 > 3.66 456.0778 > 0.51 291.13 > > > > > >> power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n, > data=stage_discharge, start=list(C=4, a=0, n=1)) >> C<-coef(power.nls)["C"] >> a<-coef(power.nls)["a"] >> n<-coef(power.nls)["n"] >> plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25, > ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n > St age-Discharge Curve') >> curve(C*(x+a)^n, add=TRUE, col="red") > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > 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. -- David K Stevens, P.E., Ph.D. Professor and Head, Environmental Engineering Civil and Environmental Engineering Utah Water Research Laboratory 8200 Old Main Hill Logan, UT 84322-8200 435 797 3229 - voice 435 797 1363 - fax david.stev...@usu.edu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Non-Linear Regression Help
If you insist on using nls() for anything that you don't understand
extremely well, you will end up with frustration. nls() uses the same
method K F Gauss used (with good understanding of the details) over
200 years ago. The Gauss-Newton approach inside works very well and
efficiently for problems where the assumptions are met, and terribly
most other times. But nls() does have some nice "extras", and rather
than rewrite all the code, we have a wrapnls() function for the codes
in the much more modern nlsr package. It tries (and mostly succeeds) in
getting analytic derivatives in cases like this. Note that nls(), when
you output the diagnostic, tells you that it is having trouble with
the numeric derivative.
I did the following:
1) made a csv file from the data in the posting (Shadomy17.csv)
2) edited the nls() call and added trace and try()
3) ran nlxb() from nlsr. Note that it uses a lot of iterations -- the
problem is quite close to singular. The singular values have NOTHING
to do with the individual parameters. Their display position is just
convenient. Together they show that the ratio of largest / smallest sv
(a measure of the condition number) is very large -- an ill-conditioned
problem. Now we know this -- there's no guessing and hand-waving.
Best, JN
Here's the (rather verbose) output
> library(nlsr)
> shadomy <- read.csv("./Shadomy17.csv")
> power.nls<-try(nls(discharge~C*(stage+a)^n, data=shadomy, start=list(C=4,
a=0, n=1), trace=TRUE))
7585285 : 4 0 1
> print(power.nls)
[1] "Error in numericDeriv(form[[3L]], names(ind), env) : \n Missing value or an infinity produced when evaluating the
model\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
> tmp <- readline("Try a better approach")
> p.nlxb <- nlxb(discharge~C*(stage+a)^n, data=shadomy, start=list(C=4, a=0,
n=1), trace=TRUE)
formula: discharge ~ C * (stage + a)^n
lower:[1] -Inf -Inf -Inf
upper:[1] Inf Inf Inf
$watch
[1] FALSE
$phi
[1] 1
$lamda
[1] 1e-04
$offset
[1] 100
$laminc
[1] 10
$lamdec
[1] 4
$femax
[1] 1
$jemax
[1] 5000
$rofftest
[1] TRUE
$smallsstest
[1] TRUE
vn:[1] "discharge" "C" "stage" "a" "n"
Finished masks check
datvar:[1] "discharge" "stage"
Data variable discharge :[1] 2592.05 559.58 484.22 494.75 456.08 291.13
Data variable stage :[1] 6.53 6.32 5.96 4.99 3.66 0.51
trjfn:
function (prm)
{
if (is.null(names(prm)))
names(prm) <- names(pvec)
localdata <- list2env(as.list(prm), parent = data)
eval(residexpr, envir = localdata)
}
no weights
lower:[1] -Inf -Inf -Inf
upper:[1] Inf Inf Inf
Start:lamda: 1e-04 SS= 7585285 at C = 4 a = 0 n = 1 1 / 0
lamda: 0.001 SS= Inf at C = -843.56 a = 228.63 n = 123.33 2 / 1
lamda: 0.01 SS= Inf at C = -515.08 a = 148.03 n = 84.714 3 / 1
lamda: 0.1 SS= 9.0129e+100 at C = -50.129 a = 37.016 n = 29.648 4 / 1
lamda: 1 SS= 8.5013e+47 at C = 58.103 a = 30.986 n = 13.954 5 / 1
lamda: 10 SS= 4.2141e+31 at C = 49.209 a = 45.025 n = 8.0734 6 / 1
lamda: 100 SS= 9.4088e+10 at C = 17.369 a = 15.101 n = 2.9571 7 / 1
< print(p.nlxb)
nlsr object: x
residual sumsquares = 2468891 on 6 observations
after 31Jacobian and 51 function evaluations
namecoeff SE tstat pval gradient
JSingval
C17.2692 1404 0.0123 0.991 -733.5
4112
a -0.509779 60.67 -0.008403 0.9938 35772
188.7
n2.44601 31.750.07704 0.9434 -126335
0.6456
> sink()
On 2017-05-04 06:58 PM, Zachary Shadomy wrote:
I am having some errors come up in my first section of code. I have no
issue in plotting the points. Is there an easier method for creating a
non-linear regression using C*(x+a)^n. The .txt file is named
stage_discharge with the two variables being stage and discharge.
The data is a relatively small file listed below:
stage discharge
6.53 2592.05
6.32 559.5782
5.96 484.2151
4.99 494.7527
3.66 456.0778
0.51 291.13
power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n,
data=stage_discharge, start=list(C=4, a=0, n=1))
C<-coef(power.nls)["C"]
a<-coef(power.nls)["a"]
n<-coef(power.nls)["n"]
plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25,
ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n
St age-Discharge Curve')
curve(C*(x+a)^n, add=TRUE, col="red")
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the pos
Re: [R] Non-Linear Regression Help
Hi I am not an expert in nonlinear regression, but your data seems to be rather weird. Last five points has almost linear relationship, the first one is several times higher. If there is no error in your data, I doubt that you can model it by power equation. Cheers Petr > -Original Message- > From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Zachary > Shadomy > Sent: Friday, May 5, 2017 12:58 AM > To: r-help@r-project.org > Subject: [R] Non-Linear Regression Help > > I am having some errors come up in my first section of code. I have no > issue in plotting the points. Is there an easier method for creating a > non-linear regression using C*(x+a)^n. The .txt file is named > stage_discharge with the two variables being stage and discharge. > The data is a relatively small file listed below: > > stage discharge > 6.53 2592.05 > 6.32 559.5782 > 5.96 484.2151 > 4.99 494.7527 > 3.66 456.0778 > 0.51 291.13 > > > > > > > power.nls<- > nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n, > data=stage_discharge, start=list(C=4, a=0, n=1)) > > C<-coef(power.nls)["C"] > > a<-coef(power.nls)["a"] > > n<-coef(power.nls)["n"] > > plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25, > ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n > St age-Discharge Curve') > > curve(C*(x+a)^n, add=TRUE, col="red") > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > 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. Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny pouze jeho adresátům. Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze svého systému. Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat. Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či zpožděním přenosu e-mailu. V případě, že je tento e-mail součástí obchodního jednání: - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a to z jakéhokoliv důvodu i bez uvedení důvodu. - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce s dodatkem či odchylkou. - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným dosažením shody na všech jejích náležitostech. - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi či osobě jím zastoupené známá. This e-mail and any documents attached to it may be confidential and are intended only for its intended recipients. If you received this e-mail by mistake, please immediately inform its sender. Delete the contents of this e-mail with all attachments and its copies from your system. If you are not the intended recipient of this e-mail, you are not authorized to use, disseminate, copy or disclose this e-mail in any manner. The sender of this e-mail shall not be liable for any possible damage caused by modifications of the e-mail or by delay with transfer of the email. In case that this e-mail forms part of business dealings: - the sender reserves the right to end negotiations about entering into a contract in any time, for any reason, and without stating any reasoning. - if the e-mail contains an offer, the recipient is entitled to immediately accept such offer; The sender of this e-mail (offer) excludes any acceptance of the offer on the part of the recipient containing any amendment or variation. - the sender insists on that the respective contract is concluded only upon an express mutual agreement on all its aspects. - the sender of this e-mail informs that he/she is not authorized to enter into any contracts on behalf of the company except for cases in which he/she is expressly authorized to do so in writing, and such authorization or power of attorney is submitted to the recipient or the person represented by the recipient, or the existence of such authorization is known to the recipient of the person represented by the recipient. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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, r
Re: [R] Non linear regression - Von Bertalanffy Growth Function - "singular gradient matrix at initial parameter estimates"
Thank you for the tip. Indeed, nlxb in nlmrt works and results are not crazy. I would like however to assess goodness-of-fit (gof) and ultimately to compare it with gof from linear regression (fitted with same variables). Before I used AICc to compare the nls() and lm() fit, however I get now an error message concerning the method loglike and its non compatibility with nlmrt class object. I guess it is because we use now Marquardt method to minimise sum-of square instead of Gauss-Newton? I am right? Or this is just an incompatibility coming between AICc function and nlmrt objects? Is there an R function to do that? Best, Xochitl C. <>< <>< <>< <>< Xochitl CORMON +33 (0)3 21 99 56 84 Doctorante en écologie marine et science halieutique PhD student in marine ecology and fishery science <>< <>< <>< <>< IFREMER Centre Manche Mer du Nord 150 quai Gambetta 62200 Boulogne-sur-Mer <>< <>< <>< <>< Le 19/08/2015 15:11, ProfJCNash a écrit : Packages nlmrt or minpack.lm use a Marquardt method. minpack.lm won't proceed if the Jacobian singularity is at the starting point as far as I'm aware, but nlxb in nlmrt can sometimes get going. It has a policy that is aggressive in trying to improve the sum of squares, so will use more effort than nls when both work. JN On 15-08-18 12:08 PM, Xochitl CORMON wrote: Dear all, I am trying to estimate VBGF parameters K and Linf using non linear regression and nls(). First I used a classic approach where I estimate both parameters together as below with "alkdyr" being a subset per year of my age-length-key database and running in a loop. vbgf.par <- nls(Lgtcm ~ Linf *(1 - exp(-K * (Age - tzero))), start = c(K= 0.07, Linf = 177.1), data=alkdyr) I obtain an estimation of both parameters that are strongly correlated. Indeed after plotting Linf ~ K and fitting a linear regression I obtain a function (Linf = a + b*K) with R2= 0.8 and a = 215, b = -763. In this context, to take into account explicitly correlation between parameters, I decided to fit a new non linear regression derivate from VBGF but where Linf is expressed depending on K (I am most interested in K). To do so, I tried this model: vbgf.par <- nls(Lgtcm ~ (a + (b*k)) *(1 - exp(-k * (Age - tzero))), start = c(k= 0.07, a= 215, b=-763), data=alkdyr) Unfortunately at this point I cannot go further as I get the error message "singular gradient matrix at initial parameter estimates". I tried to use alg= plinear (which I am not sure I understand properly yet). If I give a starting value for a and b only, I have an error message stating "step factor below minFactor" (even when minFactor is set to 1000). Any help will be more than welcome as this is quite urgent Best, Xochitl C. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Non linear regression - Von Bertalanffy Growth Function - "singular gradient matrix at initial parameter estimates"
Packages nlmrt or minpack.lm use a Marquardt method. minpack.lm won't proceed if the Jacobian singularity is at the starting point as far as I'm aware, but nlxb in nlmrt can sometimes get going. It has a policy that is aggressive in trying to improve the sum of squares, so will use more effort than nls when both work. JN On 15-08-18 12:08 PM, Xochitl CORMON wrote: Dear all, I am trying to estimate VBGF parameters K and Linf using non linear regression and nls(). First I used a classic approach where I estimate both parameters together as below with "alkdyr" being a subset per year of my age-length-key database and running in a loop. vbgf.par <- nls(Lgtcm ~ Linf *(1 - exp(-K * (Age - tzero))), start = c(K= 0.07, Linf = 177.1), data=alkdyr) I obtain an estimation of both parameters that are strongly correlated. Indeed after plotting Linf ~ K and fitting a linear regression I obtain a function (Linf = a + b*K) with R2= 0.8 and a = 215, b = -763. In this context, to take into account explicitly correlation between parameters, I decided to fit a new non linear regression derivate from VBGF but where Linf is expressed depending on K (I am most interested in K). To do so, I tried this model: vbgf.par <- nls(Lgtcm ~ (a + (b*k)) *(1 - exp(-k * (Age - tzero))), start = c(k= 0.07, a= 215, b=-763), data=alkdyr) Unfortunately at this point I cannot go further as I get the error message "singular gradient matrix at initial parameter estimates". I tried to use alg= plinear (which I am not sure I understand properly yet). If I give a starting value for a and b only, I have an error message stating "step factor below minFactor" (even when minFactor is set to 1000). Any help will be more than welcome as this is quite urgent Best, Xochitl C. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Non linear regression - Von Bertalanffy Growth Function - "singular gradient matrix at initial parameter estimates"
These appear to be primarily statistics/nonlinear optimization issues that are off topic here, which is about R programming. Post on a statistics list like stats.stackexchange.com instead. Cheers, Bert Bert Gunter "Data is not information. Information is not knowledge. And knowledge is certainly not wisdom." -- Clifford Stoll On Tue, Aug 18, 2015 at 9:08 AM, Xochitl CORMON wrote: > Dear all, > > I am trying to estimate VBGF parameters K and Linf using non linear > regression and nls(). First I used a classic approach where I estimate both > parameters together as below with "alkdyr" being a subset per year of my > age-length-key database and running in a loop. > > vbgf.par <- nls(Lgtcm ~ Linf *(1 - exp(-K * (Age - tzero))), start = c(K= > 0.07, Linf = 177.1), data=alkdyr) > > I obtain an estimation of both parameters that are strongly correlated. > Indeed after plotting Linf ~ K and fitting a linear regression I obtain a > function (Linf = a + b*K) with R2= 0.8 and a = 215, b = -763. > > In this context, to take into account explicitly correlation between > parameters, I decided to fit a new non linear regression derivate from VBGF > but where Linf is expressed depending on K (I am most interested in K). To > do so, I tried this model: > vbgf.par <- nls(Lgtcm ~ (a + (b*k)) *(1 - exp(-k * (Age - tzero))), start = > c(k= 0.07, a= 215, b=-763), data=alkdyr) > > Unfortunately at this point I cannot go further as I get the error message > "singular gradient matrix at initial parameter estimates". > > I tried to use alg= plinear (which I am not sure I understand properly yet). > If I give a starting value for a and b only, I have an error message stating > "step factor below minFactor" (even when minFactor is set to 1000). > > Any help will be more than welcome as this is quite urgent > > Best, > > Xochitl C. > > > > > -- > > <>< <>< <>< <>< > > Xochitl CORMON > +33 (0)3 21 99 56 84 > > Doctorante en écologie marine et science halieutique > PhD student in marine ecology and fishery science > > <>< <>< <>< <>< > > IFREMER > Centre Manche Mer du Nord > 150 quai Gambetta > 62200 Boulogne-sur-Mer > > <>< <>< <>< <>< > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > 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 -- To UNSUBSCRIBE and more, see 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.
Re: [R] Non-linear regression analysis in R
Jeremy: Don't be silly. The model is overdetermined -- a and b "tradeoff" with each other. e.g. for any solution (a,b), (a/k^m,b+m) for any m is also a solution, where k = const (assuming I have correctly interpreted the model, of course). -- Bert On Wed, Dec 19, 2012 at 11:59 AM, Jeremy Miles wrote: > Could you provide the code that you're running, so we can see what > you're trying to do? Even better would be a repeatable example. > > Jeremy > > On 19 December 2012 09:42, Yann Labou wrote: > > Hey all, > > > > I'm trying to fit a non-linear model y ~ a * constant ^ b * x ^ c and > > estimates the paramaters a, b and c. > > > > Using the nls function, I'm getting following error message: > > > > Error in nlsModel(formula, mf, start, wts) : > > singular gradient matrix at initial parameter estimates > > > > If I logarithmize the whole equation log(y) ~ log (a) + b * > log(constant) + > > c * log(x) and fit the equation with lm, I get only NAs estimates for the > > second term on the right side: > > > > Coefficients: (1 not defined because of singularities) > > > > Do you have any hints on how to fit this equation or any alternatives to > > nls? > > > > Thanks, > > Yann > > > > __ > > 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. > -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm [[alternative HTML version deleted]] __ 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.
Re: [R] Non-linear regression analysis in R
Could you provide the code that you're running, so we can see what you're trying to do? Even better would be a repeatable example. Jeremy On 19 December 2012 09:42, Yann Labou wrote: > Hey all, > > I'm trying to fit a non-linear model y ~ a * constant ^ b * x ^ c and > estimates the paramaters a, b and c. > > Using the nls function, I'm getting following error message: > > Error in nlsModel(formula, mf, start, wts) : > singular gradient matrix at initial parameter estimates > > If I logarithmize the whole equation log(y) ~ log (a) + b * log(constant) + > c * log(x) and fit the equation with lm, I get only NAs estimates for the > second term on the right side: > > Coefficients: (1 not defined because of singularities) > > Do you have any hints on how to fit this equation or any alternatives to > nls? > > Thanks, > Yann > > __ > 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.
Re: [R] Non linear regression with complex equation
Apologies for the phrasing of the question.
I've sorted the problem (thanks Bert Gunter) by using the curly brackets {}
as below (using a simplified version of my real model). I hope this helps
someone else!
Jeff
---
> data
Alpha ip X
1 0.7106967 0.3616727 0.006027879
2 2.1678517 5.0615917 0.084359861
3 4.4066250 11.2282945 0.187138242
4 9.8495694 18.0534974 0.300891624
527.7247098 29.2064434 0.486774057
670.6931430 35.3946092 0.589910153
7 133.1240255 46.0347288 0.767245480
8 214.7851844 49.3811149 0.823018582
9 359.5511036 58.5069583 0.975115972
10 748.1840127 57.3744477 0.956240795
11 2129.9844080 60.000 1.0
> c<-1.83e-9
> cFe=c
> model<-nls({Fe1<-cFe/(Alpha+1+k*c)
+ X~Alpha/(Alpha+1+k*c/(1+k*Fe1))},start=list(k=1e10))
> summary(model)
Formula: X ~ Alpha/(Alpha + 1 + k * c/(1 + k * Fe1))
Parameters:
Estimate Std. Error t value Pr(>|t|)
k 3.491e+10 7.190e+09 4.856 0.000665 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05589 on 10 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 2.393e-06
--
View this message in context:
http://r.789695.n4.nabble.com/Non-linear-regression-with-complex-equation-tp4425942p4427617.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.
Re: [R] Non linear regression with complex equation
Your question is (completely) ill-posed. What is your actual model? What you have said makes no sense at all as it stands. (Minimal self-contained example .) cheers, Rolf Turner On 28/02/12 09:25, jeff_hawkes wrote: Hi all, Is it possible to model a function where the unknown parameter appears both in the fitted equation AND in the determination of other parameters? E.g. y = a^2 + b/2 + k where a = 2/k and b = k^2 and the model needs to determine k? I know this is a very simple equation (its just an example), the one I am modelling is much more complicated! k appears in the equation which the n.l.r model fits, but it also affects other parameters in the equation. Please let me know if you know a way of achieving this. I realise it is possible to set up a loop where the modelled value for k is fed back in to a and b, and the model is run again - but it seems like there should be a more elegant way within one run of the model. __ 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.
Re: [R] Non-linear Regression best-fit line
On Jun 17, 2011, at 23:14 , Sean Bignami wrote:
> I am trying to fit a curve to a cumulative mortality curve (logistic) where y
> is the cumulative proportion of mortalities, and t is the time in hours (see
> below). Asym. at 0 and 1
>> y
> [1] 0. 0.04853859 0.08303777 0.15201970 0.40995074 0.46444992
> 0.62862069 0.95885057 1.
> [10] 1. 1.
>> t
> [1] 0 13 20 24 37 42 48 61 72 86 90
>
> I tried to find starting values for a Gompertz non-linear regression by
> converting the equation (y~1*exp(-beta*exp(-gamma*t)) to a linear form per
> Dalgaard "Introductory Statistics with R" pg.279-280. But got this Error
> message:
>
>> lm(log(0-log(y))~t)
> Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
> NA/NaN/Inf in foreign function call (arg 4)
>
> I tried to change all by 0 and 1 values to non-zero and non-one values (yy
> and tt below), and was able to get starting estimates.
>
>> yy<-c(0.1,0.04853859, 0.08303777, 0.15201970, 0.40995074,
>> 0.46444992, 0.62862069, 0.95885057,
>> 0.99,0.,0.)
>> tt<-c(0.01,13,20,24,37,42,48,61,72,86,90)
>
>> lm(log(0-log(yy))~tt)
>
> Call:
> lm(formula = log(0 - log(yy)) ~ tt)
>
> Coefficients:
> (Intercept) tt
> 9.5029 -0.3681
>
> However, when I plug those values into the nls function, I get an error
> message about the "getInitial" method
>
>> nlsout<-nls(y~1*exp(-beta*exp(-gamma*t),start=c(beta=exp(9.5),gamma=.368)));summary(nlsout)
> Error in getInitial.default(func, data, mCall = as.list(match.call(func, :
> no 'getInitial' method found for "function" objects
Not really, but getting your parentheses right on the nls call should at least
get you closer. There are 3 "(" before "start" but only 1 ")", so the "start"
bit is part of the formula, etc.
Trying that, I got
> nlsout<-nls(yy~1*exp(-beta*exp(-gamma*tt)),start=list(beta=exp(9.5),gamma=.368))
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
so you still have some way to go. (I've been using the yy/tt variables because
they were easier to paste in from your post. Shouldn't matter much.)
Offhand, I'd say that your procedure for finding starting values is still too
sensitive to the zeros and ones. If you do
plot(tt, log(-log(yy)) )
you will see that the last three points (the 1s in your original data) fall way
off the linear trend. So how about omitting them rather than putting in an
arbitrary value? And get rid of the zero too.
> lm(log(-log(yy))~tt,subset=2:8)
Call:
lm(formula = log(-log(yy)) ~ tt, subset = 2:8)
Coefficients:
(Intercept) tt
2.5870 -0.0807
> nlsout<-nls(yy~1*exp(-beta*exp(-gamma*tt)),start=list(beta=exp(2.5),gamma=.08))
> summary(nlsout)
Formula: yy ~ 1 * exp(-beta * exp(-gamma * tt))
Parameters:
Estimate Std. Error t value Pr(>|t|)
beta 13.146344.617672.850.019 *
gamma 0.072840.008698.38 1.5e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0568 on 9 degrees of freedom
Number of iterations to convergence: 10
Achieved convergence tolerance: 4.29e-06
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@gmail.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.
Re: [R] Non-linear Regression best-fit line
Hi: Perhaps the self-starting functions may be helpful. See ?selfStart. There are self-starting functions for both the logistic and Gompertz models (SSlogis and SSgompertz, respectively). Go through the examples to see how they work. HTH, Dennis On Fri, Jun 17, 2011 at 2:14 PM, Sean Bignami wrote: > I am trying to fit a curve to a cumulative mortality curve (logistic) where y > is the cumulative proportion of mortalities, and t is the time in hours (see > below). Asym. at 0 and 1 >> y > [1] 0. 0.04853859 0.08303777 0.15201970 0.40995074 0.46444992 > 0.62862069 0.95885057 1. > [10] 1. 1. >> t > [1] 0 13 20 24 37 42 48 61 72 86 90 > > I tried to find starting values for a Gompertz non-linear regression by > converting the equation (y~1*exp(-beta*exp(-gamma*t)) to a linear form per > Dalgaard "Introductory Statistics with R" pg.279-280. But got this Error > message: > >> lm(log(0-log(y))~t) > Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : > NA/NaN/Inf in foreign function call (arg 4) > > I tried to change all by 0 and 1 values to non-zero and non-one values (yy > and tt below), and was able to get starting estimates. > >> yy<-c(0.1,0.04853859, 0.08303777, 0.15201970, 0.40995074, >> 0.46444992, 0.62862069, 0.95885057, >> 0.99,0.,0.) >> tt<-c(0.01,13,20,24,37,42,48,61,72,86,90) > >> lm(log(0-log(yy))~tt) > > Call: > lm(formula = log(0 - log(yy)) ~ tt) > > Coefficients: > (Intercept) tt > 9.5029 -0.3681 > > However, when I plug those values into the nls function, I get an error > message about the "getInitial" method > >> nlsout<-nls(y~1*exp(-beta*exp(-gamma*t),start=c(beta=exp(9.5),gamma=.368)));summary(nlsout) > Error in getInitial.default(func, data, mCall = as.list(match.call(func, : > no 'getInitial' method found for "function" objects > > Can anyone help clarify how I can find the parameters for a best-fit curve > for this data? Thanks!! > > Sean > [[alternative HTML version deleted]] > > __ > 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.
Re: [R] Non linear Regression: "singular gradient matrix at initial parameter estimates"
Hi Mario, yes works great. Thanks! 2011/4/12 Mario Valle > Use a more realistic starting point instead of the default one: > > fit <- nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x), > start=list(p1=410,p2=18,p4=-.03)) > > This works for me: > > fit > Nonlinear regression model > model: yeps ~ p1/(1 + exp(p2 - x)) * exp(p4 * x) > data: parent.frame() > p1p2p4 > 199.48276 16.28664 -0.01987 > residual sum-of-squares: 560.6 > > Number of iterations to convergence: 5 > Achieved convergence tolerance: 5.637e-07 > > Ciao! >mario > > > On 12-Apr-11 18:01, Felix Nensa wrote: > >> fit = nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x)) >> >> > -- > Ing. Mario Valle > Data Analysis and Visualization Group| > http://www.cscs.ch/~mvalle > Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60 > v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82 > > -- Felix Nensa Luisenstr. 15-17 44787 Bochum Germany mail: felix.ne...@googlemail.com mobile: +49 171 958 51 40 [[alternative HTML version deleted]] __ 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.
Re: [R] Non linear Regression: "singular gradient matrix at initial parameter estimates"
Use a more realistic starting point instead of the default one: fit <- nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x), start=list(p1=410,p2=18,p4=-.03)) This works for me: > fit Nonlinear regression model model: yeps ~ p1/(1 + exp(p2 - x)) * exp(p4 * x) data: parent.frame() p1p2p4 199.48276 16.28664 -0.01987 residual sum-of-squares: 560.6 Number of iterations to convergence: 5 Achieved convergence tolerance: 5.637e-07 Ciao! mario On 12-Apr-11 18:01, Felix Nensa wrote: fit = nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x)) -- Ing. Mario Valle Data Analysis and Visualization Group| http://www.cscs.ch/~mvalle Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60 v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82 __ 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.
Re: [R] Non linear Regression: "singular gradient matrix at initial parameter estimates"
Hi Peter, thank you for your reply. Now I see, that P3 is indeed redundand. But with the simplified model... fit = nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x)) ...nls still produces the same error. Any ideas? Felix 2011/4/12 Peter Ehlers > On 2011-04-11 13:29, Felix Nensa wrote: > >> Hi, >> >> I am using nls to fit a non linear function to some data but R keeps >> giving >> me "singular gradient matrix at initial parameter estimates" errors. >> For testing purposes I am doing this: >> >> ### R code ### >> >> x<- 0:140 >> y<- 200 / (1 + exp(17 - x)/2) * exp(-0.02*x) # creating 'perfect' samples >> with fitting model >> yeps<- y + rnorm(length(y), sd = 2) # adding noise >> >> # results in above error >> fit = nls(yeps ~ p1 / (1 + exp(p2 - x) / p3) * exp(p4 * x)) >> >> ### >> >> From what I've found in this list I think that my model is >>> over-parameterized. >>> >> How can I work around that? >> > > Take out p3; it's redundant. > > Peter Ehlers > > Thanks, >> >> Felix >> >>[[alternative HTML version deleted]] >> >> __ >> 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. >> > > [[alternative HTML version deleted]] __ 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.
Re: [R] Non linear Regression: "singular gradient matrix at initial parameter estimates"
On 2011-04-11 13:29, Felix Nensa wrote: Hi, I am using nls to fit a non linear function to some data but R keeps giving me "singular gradient matrix at initial parameter estimates" errors. For testing purposes I am doing this: ### R code ### x<- 0:140 y<- 200 / (1 + exp(17 - x)/2) * exp(-0.02*x) # creating 'perfect' samples with fitting model yeps<- y + rnorm(length(y), sd = 2) # adding noise # results in above error fit = nls(yeps ~ p1 / (1 + exp(p2 - x) / p3) * exp(p4 * x)) ### From what I've found in this list I think that my model is over-parameterized. How can I work around that? Take out p3; it's redundant. Peter Ehlers Thanks, Felix [[alternative HTML version deleted]] __ 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.
Re: [R] non-linear regression for 3D data
Thanks a lot, this I(xx^2) ... worked. I guess, I should learn more abot the function poly itself. (so will I... :) ) Thanks again! -- View this message in context: http://r.789695.n4.nabble.com/non-linear-regression-for-3D-data-tp2320982p2325911.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.
Re: [R] non-linear regression for 3D data
On 12/08/2010 10:35 AM, szisziszilvi wrote: I've tried lm, but something is wrong. I've made a test dataset of 599 data points, my original equation is zz = 1 +0.5*xx -3.2*xx*xx -1*yy +4.2*yy*yy but the R gives this result: --- > mp <- read.csv(file="sample.csv",sep=";",header=TRUE) > lm(zz ~ poly(xx,2) + poly(yy,2), data=mp) Call: lm(formula = zz ~ poly(xx, 2) + poly(yy, 2), data = mp) Coefficients: (Intercept) poly(xx, 2)1 poly(xx, 2)2 poly(yy, 2)1 poly(yy, 2)2 25.86 -2239.86 -595.01 2875.54776.84 --- which is definitely not the original. :( I don't think you are interpreting the coefficients properly. The basis functions are orthogonal polynomials, not xx and xx^2, so the coefficients won't match the ones you used in your definition. You should compare the predictions of the model, e.g. by looking at range(predict(lm(zz ~ poly(xx,2) + poly(yy,2), data=mp)) - zz) If you insist on the power basis, just fit the model as lm(zz ~ xx + I(xx^2) + yy + I(yy^2), data=mp) but you might get less accurate predictions due to increased collinearity. Duncan Murdoch __ 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.
Re: [R] non-linear regression for 3D data
right. How does it come that if I devide the result vector with 10*interception, I get a much better result? > zz2 <- 25.86 -2239.86*mp$xx -595.01*mp$xx*mp$xx + 2875.54*mp$yy + > 776.84*mp$yy*mp$yy > mp$zz2 <- zz2 > library(lattice) > cloud(zz2/258.6 + zz ~ xx * yy, data=mp) looks quite pretty. http://r.789695.n4.nabble.com/file/n2322812/output.jpeg -- View this message in context: http://r.789695.n4.nabble.com/non-linear-regression-for-3D-data-tp2320982p2322812.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.
Re: [R] non-linear regression for 3D data
I've tried lm, but something is wrong. I've made a test dataset of 599 data points, my original equation is zz = 1 +0.5*xx -3.2*xx*xx -1*yy +4.2*yy*yy but the R gives this result: --- > mp <- read.csv(file="sample.csv",sep=";",header=TRUE) > lm(zz ~ poly(xx,2) + poly(yy,2), data=mp) Call: lm(formula = zz ~ poly(xx, 2) + poly(yy, 2), data = mp) Coefficients: (Intercept) poly(xx, 2)1 poly(xx, 2)2 poly(yy, 2)1 poly(yy, 2)2 25.86 -2239.86 -595.01 2875.54776.84 --- which is definitely not the original. :( (In case of interest the test dataset is available here: szisziszilvi.lima-city.de/r) -- View this message in context: http://r.789695.n4.nabble.com/non-linear-regression-for-3D-data-tp2320982p2322779.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.
Re: [R] non-linear regression for 3D data
oh, god, please don't tell anybody... -- View this message in context: http://r.789695.n4.nabble.com/non-linear-regression-for-3D-data-tp2320982p2321082.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.
Re: [R] non-linear regression for 3D data
On 11/08/2010 6:15 AM, szisziszilvi wrote: Hello! Is there a simplier way in R to get a nonlinear regression (like nls) for a surface? I have 3D data, and it is definitely not a linear surface with which it would fit the best. Rather sg like z = a + f(x) + g(y) where probably both f and g are polinomes (hopefully quadratic). That's a linear model if f and g are polynomials. You can fit it with lm(z ~ poly(x, 2) + poly(y, 2)) The 2's are the degrees of each polynomial; the intercept a is implied. Duncan Murdoch __ 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.
Re: [R] Non-linear regression
It appears my suspicions about this being homework were unfounded. Given the additional problems with excess zeroes, you may want to examine the extremely informative material on analysis of such problems written by Zeileis, Kleiber and Jackman: (easily found in case you have misplaced it, as I had, with a Google search for: "r-project" zero-inflated hurdle models "Regression Models for Count Data in R" http://cran.cnr.berkeley.edu/web/packages/pscl/vignettes/countreg.pdf -- David. On Feb 6, 2010, at 10:56 PM, kupz wrote: Agreed, it would be simple to propose the relationship, however the regression is necessary to model the data properly. Unfortunately a simple decay based on those two points does not have the proper shape necessary. This is due to an extreme amount of zero inflation with this fisheries data. On another note, I have a working solution for the problem, I am excluding a portion of the zero data based on some other apriori assumptions.. Thanks for your help though. -- View this message in context: http://n4.nabble.com/Non-linear-regression-tp1471736p1471749.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. David Winsemius, MD Heritage Laboratories West Hartford, CT __ 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.
Re: [R] Non-linear regression
Agreed, it would be simple to propose the relationship, however the regression is necessary to model the data properly. Unfortunately a simple decay based on those two points does not have the proper shape necessary. This is due to an extreme amount of zero inflation with this fisheries data. On another note, I have a working solution for the problem, I am excluding a portion of the zero data based on some other apriori assumptions.. Thanks for your help though. -- View this message in context: http://n4.nabble.com/Non-linear-regression-tp1471736p1471749.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.
Re: [R] Non-linear regression
On Feb 6, 2010, at 10:33 PM, kupz wrote: So I have a data set I would like to model using a non-linear method. I know it should be an exponential decay. However I know what the first derivative of the equation should be at two points, x=0 and x=100. Is there anyway to establish this when inputing the model or how would one go about this before the nls statement Given the rather simple relationship between the exponential function and its derivative, why would you need the regression if you already have those two points for dy/dx ( as well as the value of the function at x=0)? Is this homework? -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ 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.
Re: [R] non-linear regression
>
> The problem of estimation of parameters in R is that you have to know the
> value of the initial estimates very accurately, otherwise it does not
> converge.
This was discussed on R-help in the last 2 weeks; see the thread on
'Starting estimates for nls Exponential Fit'.
>
> The example below could be resolved in Excel, however in does not converge.
> How to solve the problem?
Can you consider a model with a different functional form or do you need
to fit this exact function?
If you want to minimize log(data) log(model) RSS, then:
tx.br <- read.table('tx.br.H.txt',header=F,dec=',')
tx.br <- tx.br[,1]
id <- 1:100
qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
newD <- log(tx.br)
HP <- nls(newD~log(qx.suav(id,A,B,C,D,E,F,G,H,K)),
data=data.frame(id=id,newD=newD),
trace=TRUE,
nls.control(maxiter=5,warnOnly=TRUE),
algorithm='port',
start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
E=5.949082737,F=24.526811,G=0.46733960,H=1.0970550987,
K=0.771722501657),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
HP
par(mfrow=c(2,1))
matplot(cbind(fitted(HP), newD),type="l", main="model and fit")
matplot(cbind(exp(fitted(HP)), exp(newD)), type="l",
main="transformed back to original space")
> I made the chart on a logarithmic scale to better visualize the differences.
>
> Send the data file attached.
>
> The commands are below:
>
> tx.br <- read.table('c:/tx.br.H.txt',header=F,dec=',')
> tx.br <-tx.br[,1]
> id<-1:100
>
> qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
> (A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
>
> HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
> data=data.frame(id=id,tx.br=tx.br),
> trace=TRUE,nls.control(maxiter=5,warnOnly=TRUE,minFactor =
> 0.1),
> algorithm='port',
> start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
>
> E=5.949082737,F=24.526811,G=0.46733960,H=1.0970550987,K=0.771722501657),
> lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
>
> HP
>
> matplot(cbind(log(fitted(HP)), log(tx.br)),type="l")
>
>
>
> - Original Message -
> From: "Katharine Mullen"
> To: "AneSR"
> Cc:
> Sent: Thursday, December 10, 2009 9:55 PM
> Subject: Re: [R] non-linear regression
>
>
> > You did not provide the data or a way of generating it.
> >
> > I would guess that Excel finds the same solution (the same residual sum-of
> > squares) as nls but that it uses a weak convergence criterion and/or does
> > not give you information regarding why it terminates.
> >
> > Regarding the step size: you can set the minimum step size factor via the
> > minFactor argument of control.
> >
> > qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
> > (A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
> >
> > ## make noisy data from model
> > id <- 1:1000
> > tx.br <- qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
> > E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
> > K=0.382070638)
> > set.seed(1)
> > tx.br <- tx.br + rnorm(length(tx.br),0,.0001)
> >
> > HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
> > data=data.frame(id=id,tx.br=tx.br),
> > trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
> > algorithm='port',
> > start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
> >E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
> >K=0.382070638),
> > lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
> > matplot(cbind(fitted(HP), tx.br),type="l")
> >
> > On Thu, 10 Dec 2009, AneSR wrote:
> >
> >>
> >> I have a non-linear regression with 8 parameters to solve however it
> >> does not converge ... easily solves the excel ... including the initial
> >> estimates used in the R were found in the excel ... Another question is
> >> how
> >> to establish the increments of R by the parameters in the search ..
> >>
> >>
> >> qx.suav<-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
> >> HP<-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
> >> trace=TRUE,nls.control(maxiter=5000),algorithm='
Re: [R] non-linear regression
Katharine,
The problem of estimation of parameters in R is that you have to know the
value of the initial estimates very accurately, otherwise it does not
converge.
The example below could be resolved in Excel, however in does not converge.
How to solve the problem?
I made the chart on a logarithmic scale to better visualize the differences.
Send the data file attached.
The commands are below:
tx.br <- read.table('c:/tx.br.H.txt',header=F,dec=',')
tx.br <-tx.br[,1]
id<-1:100
qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5,warnOnly=TRUE,minFactor =
0.1),
algorithm='port',
start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
E=5.949082737,F=24.526811,G=0.46733960,H=1.0970550987,K=0.771722501657),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
HP
matplot(cbind(log(fitted(HP)), log(tx.br)),type="l")
- Original Message -
From: "Katharine Mullen"
To: "AneSR"
Cc:
Sent: Thursday, December 10, 2009 9:55 PM
Subject: Re: [R] non-linear regression
You did not provide the data or a way of generating it.
I would guess that Excel finds the same solution (the same residual sum-of
squares) as nls but that it uses a weak convergence criterion and/or does
not give you information regarding why it terminates.
Regarding the step size: you can set the minimum step size factor via the
minFactor argument of control.
qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
## make noisy data from model
id <- 1:1000
tx.br <- qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638)
set.seed(1)
tx.br <- tx.br + rnorm(length(tx.br),0,.0001)
HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
algorithm='port',
start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
matplot(cbind(fitted(HP), tx.br),type="l")
On Thu, 10 Dec 2009, AneSR wrote:
I have a non-linear regression with 8 parameters to solve however it
does not converge ... easily solves the excel ... including the initial
estimates used in the R were found in the excel ... Another question is
how
to establish the increments of R by the parameters in the search ..
qx.suav<-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
HP<-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000),algorithm='port',start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,K=0.382070638),lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
summary(HP)
How to solve this problem in R?
Ane
--
View this message in context:
http://n4.nabble.com/non-linear-regression-tp959977p959977.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
__
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.
No virus found in this incoming message.
Checked by AVG - www.avg.com
07:36:00
0,000757
0,000356
0,000290
0,000232
0,000199
0,000283
0,000128
0,000235
0,000230
0,000179
0,000197
0,000172
0,000294
0,000307
0,000411
0,000619
0,000812
0,000949
0,001116
0,001134
0,001159
0,001151
0,001321
0,001235
0,001296
0,001287
0,001241
0,001469
0,001278
0,001466
0,001440
0,001881
0,001457
0,001843
0,001754
0,001796
0,001915
0,002199
0,002120
0,002471
0,002213
0,002478
0,002975
0,002583
0,002751
0,003388
0,003956
0,004369
0,004062
0,004410
0,004939
0,005492
0,005961
0,006311
0,006846
0,007637
0,008968
0,010612
0,010907
0,011986
0,013174
0,014316
0,015557
0,017824
0,018740
0,020338
0,022034
0,022339
0,025929
0,030967
0,033520
0,039556
0,041168
0,044811
0,049932
0,047401
0,060318
0,063559
0,064084
0,077614
0,080425
0,094073
0,102383
0,116288
0,112727
0,114192
0,137330
0,159065
0,138378
0,135844
0,148244
0,165927
0,180803
0,217185
0,241570
0,249358
0,309202
0,269479
0,325142
0,461620
__
R-help@r-project.org mailing list
https://stat.ethz
Re: [R] non-linear regression
You did not provide the data or a way of generating it.
I would guess that Excel finds the same solution (the same residual sum-of
squares) as nls but that it uses a weak convergence criterion and/or does
not give you information regarding why it terminates.
Regarding the step size: you can set the minimum step size factor via the
minFactor argument of control.
qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
## make noisy data from model
id <- 1:1000
tx.br <- qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638)
set.seed(1)
tx.br <- tx.br + rnorm(length(tx.br),0,.0001)
HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
algorithm='port',
start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
matplot(cbind(fitted(HP), tx.br),type="l")
On Thu, 10 Dec 2009, AneSR wrote:
>
> I have a non-linear regression with 8 parameters to solve however it
> does not converge ... easily solves the excel ... including the initial
> estimates used in the R were found in the excel ... Another question is how
> to establish the increments of R by the parameters in the search ..
>
>
> qx.suav<-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
> HP<-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
> trace=TRUE,nls.control(maxiter=5000),algorithm='port',start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,K=0.382070638),lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
> summary(HP)
>
> How to solve this problem in R?
>
> Ane
> --
> View this message in context:
> http://n4.nabble.com/non-linear-regression-tp959977p959977.html
> Sent from the R help mailing list archive at Nabble.com.
>
> [[alternative HTML version deleted]]
>
> __
> 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.
Re: [R] Non-Linear Regression with two Predictors
Hi and thank you for your reply, in my new regression formular the parameter delta is inserted: fit <- nls(dataset$V2~(( alpha / ( 1 + exp( beta - gamma* dataset$V1 ) ) ) + (dataset$V6*delta)),data=dataset,start=startparam) The sense is, that dataset$V6 is a dummy variable that represents the german reunion. I expect, that the delta in the regression is about 1.000.000 because of the reunion. The logistic function thus has a jump at this point. But I would like to get the exact paramter value for delta (the jump) as the parameters for the logistic function of growth (alpha to gamma). The partial derivative to delta would be like a stair-function. It is 0 until 1990 and 1 there after. Any idea? Thank you! Regards Moshe Olshansky schrieb: > Hi, > > I believe that since delta does not appear in the function you are > optimizing, it's partial derivative with respect to delta is always 0 and so > the gradient is singular. > Why do you need delta at all? > > --- On Mon, 27/7/09, Berlinerfee wrote: > > >> From: Berlinerfee >> Subject: [R] Non-Linear Regression with two Predictors >> To: r-h...@stat.math.ethz.ch >> Received: Monday, 27 July, 2009, 2:52 AM >> Hello there, >> >> I am using nls the first time for a non-linear regression >> with a logistic growth function: >> startparam <- c(alpha=3e+07,beta=4000,gamma=2) >> fit <- nls(dataset$V2~(( alpha / ( 1 + exp( beta - gamma >> * dataset$V1 ) ) ) ),data=dataset,start=startparam) >> >> Everything works fine and i get good results. Now I would >> like to improve the results using my DUMMY Variable >> (dataset$V6) the runs half of the time 0 and then 1. This is >> my new nls: >> startparam <- >> c(alpha=3e+07,beta=4000,gamma=2,delta=100) >> fit <- nls(dataset$V2~(( alpha / ( 1 + exp( beta - gamma >> * dataset$V1 ) ) ) + (dataset$V6*dataset$V1*delta) >> ),data=dataset,start=startparam) >> >> I get "Singular Gradient Matrice". May anyone give me the >> right nls function for this problem?? >> >> Regards >> >> __ >> 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. >> >> > > [[alternative HTML version deleted]] __ 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.
Re: [R] Non-linear regression/Quantile regression
The coefficients are different but the predictions are the same. On Tue, Jun 9, 2009 at 12:41 PM, Greg Snow wrote: > poly by default uses orthogonal polynomials which work better mathematically > but are harder to interpret. See ?poly > > -- > Gregory (Greg) L. Snow Ph.D. > Statistical Data Center > Intermountain Healthcare > greg.s...@imail.org > 801.408.8111 > > >> -Original Message- >> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- >> project.org] On Behalf Of despaired >> Sent: Tuesday, June 09, 2009 9:59 AM >> To: r-help@r-project.org >> Subject: Re: [R] Non-linear regression/Quantile regression >> >> >> Hi, >> >> thanks, it works :-) >> But where is the difference between demand ~ Time + I(Time^2) and >> demand ~ >> poly(Time, 2) ? >> Or: How do I have to interpret the results? (I get different results >> for the >> two methods) >> >> Thank you again! >> >> >> Gabor Grothendieck wrote: >> > >> > Those are linear in the coefficients so try these: >> > >> > library(quantreg) >> > >> > rq1 <- rq(demand ~ Time + I(Time^2), data = BOD, tau= 1:3/4); rq1 >> > >> > # or >> > rq2 <- rq(demand ~ poly(Time, 2), data = BOD, tau = 1:3/4); rq2 >> > >> > >> > On Tue, Jun 9, 2009 at 10:55 AM, despaired >> wrote: >> >> >> >> Hi, >> >> >> >> I'm relatively new to R and need to do a quantile regression. Linear >> >> quantile regression works, but for my data I need some quadratic >> >> function. >> >> So I guess, I have to use a nonlinear quantile regression. I tried >> the >> >> example on the help page for nlrq with my data and it worked. But >> the >> >> example there was with a SSlogis model. Trying to write >> >> >> >> dat.nlrq <- nlrq(BM ~ I(Regen100^2), data=dat, tau=0.25, trace=TRUE) >> >> >> >> or >> >> >> >> dat.nlrq <- nlrq(BM ~ poly(Regen100^2), data=dat, tau=0.25, >> trace=TRUE) >> >> >> >> (I don't know the difference) both gave me the following error >> message: >> >> >> >> error in getInitial.default(func, data, mCall = >> as.list(match.call(func, >> >> : >> >> no 'getInitial' method found for "function" objects >> >> >> >> Looking in getInitial, it must have to do something with the >> starting >> >> parameters or selfStart model. But I have no idea, what this is and >> how I >> >> handle this problem. Can anyone please help? >> >> >> >> Thanks a lot in advance! >> >> -- >> >> View this message in context: >> >> http://www.nabble.com/Non-linear-regression-Quantile-regression- >> tp23944530p23944530.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. >> > >> > >> >> -- >> View this message in context: http://www.nabble.com/Non-linear- >> regression-Quantile-regression-tp23944530p23945900.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. > __ 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.
Re: [R] Non-linear regression/Quantile regression
poly by default uses orthogonal polynomials which work better mathematically but are harder to interpret. See ?poly -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 > -Original Message- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > project.org] On Behalf Of despaired > Sent: Tuesday, June 09, 2009 9:59 AM > To: r-help@r-project.org > Subject: Re: [R] Non-linear regression/Quantile regression > > > Hi, > > thanks, it works :-) > But where is the difference between demand ~ Time + I(Time^2) and > demand ~ > poly(Time, 2) ? > Or: How do I have to interpret the results? (I get different results > for the > two methods) > > Thank you again! > > > Gabor Grothendieck wrote: > > > > Those are linear in the coefficients so try these: > > > > library(quantreg) > > > > rq1 <- rq(demand ~ Time + I(Time^2), data = BOD, tau= 1:3/4); rq1 > > > > # or > > rq2 <- rq(demand ~ poly(Time, 2), data = BOD, tau = 1:3/4); rq2 > > > > > > On Tue, Jun 9, 2009 at 10:55 AM, despaired > wrote: > >> > >> Hi, > >> > >> I'm relatively new to R and need to do a quantile regression. Linear > >> quantile regression works, but for my data I need some quadratic > >> function. > >> So I guess, I have to use a nonlinear quantile regression. I tried > the > >> example on the help page for nlrq with my data and it worked. But > the > >> example there was with a SSlogis model. Trying to write > >> > >> dat.nlrq <- nlrq(BM ~ I(Regen100^2), data=dat, tau=0.25, trace=TRUE) > >> > >> or > >> > >> dat.nlrq <- nlrq(BM ~ poly(Regen100^2), data=dat, tau=0.25, > trace=TRUE) > >> > >> (I don't know the difference) both gave me the following error > message: > >> > >> error in getInitial.default(func, data, mCall = > as.list(match.call(func, > >> : > >> no 'getInitial' method found for "function" objects > >> > >> Looking in getInitial, it must have to do something with the > starting > >> parameters or selfStart model. But I have no idea, what this is and > how I > >> handle this problem. Can anyone please help? > >> > >> Thanks a lot in advance! > >> -- > >> View this message in context: > >> http://www.nabble.com/Non-linear-regression-Quantile-regression- > tp23944530p23944530.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. > > > > > > -- > View this message in context: http://www.nabble.com/Non-linear- > regression-Quantile-regression-tp23944530p23945900.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.
Re: [R] Non-linear regression/Quantile regression
Try `poly(Time, 2, raw=TRUE)' Here is an example: Time <- runif(100) demand <- 1 + 0.5 * Time - 1.2 * Time^2 + rt(100, df=4) rq1 <- rq(demand ~ Time + I(Time^2), tau = 1:3/4) rq2 <- rq(demand ~ poly(Time, 2, raw=TRUE), tau = 1:3/4) all.equal(c(rq1$coef), c(rq2$coef)) Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of despaired Sent: Tuesday, June 09, 2009 11:59 AM To: r-help@r-project.org Subject: Re: [R] Non-linear regression/Quantile regression Hi, thanks, it works :-) But where is the difference between demand ~ Time + I(Time^2) and demand ~ poly(Time, 2) ? Or: How do I have to interpret the results? (I get different results for the two methods) Thank you again! Gabor Grothendieck wrote: > > Those are linear in the coefficients so try these: > > library(quantreg) > > rq1 <- rq(demand ~ Time + I(Time^2), data = BOD, tau= 1:3/4); rq1 > > # or > rq2 <- rq(demand ~ poly(Time, 2), data = BOD, tau = 1:3/4); rq2 > > > On Tue, Jun 9, 2009 at 10:55 AM, despaired wrote: >> >> Hi, >> >> I'm relatively new to R and need to do a quantile regression. Linear >> quantile regression works, but for my data I need some quadratic >> function. >> So I guess, I have to use a nonlinear quantile regression. I tried >> the example on the help page for nlrq with my data and it worked. But >> the example there was with a SSlogis model. Trying to write >> >> dat.nlrq <- nlrq(BM ~ I(Regen100^2), data=dat, tau=0.25, trace=TRUE) >> >> or >> >> dat.nlrq <- nlrq(BM ~ poly(Regen100^2), data=dat, tau=0.25, >> trace=TRUE) >> >> (I don't know the difference) both gave me the following error message: >> >> error in getInitial.default(func, data, mCall = >> as.list(match.call(func, >> : >> no 'getInitial' method found for "function" objects >> >> Looking in getInitial, it must have to do something with the starting >> parameters or selfStart model. But I have no idea, what this is and >> how I handle this problem. Can anyone please help? >> >> Thanks a lot in advance! >> -- >> View this message in context: >> http://www.nabble.com/Non-linear-regression-Quantile-regression-tp239 >> 44530p23944530.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. > > -- View this message in context: http://www.nabble.com/Non-linear-regression-Quantile-regression-tp23944530p2 3945900.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.
Re: [R] Non-linear regression/Quantile regression
Hi, thanks, it works :-) But where is the difference between demand ~ Time + I(Time^2) and demand ~ poly(Time, 2) ? Or: How do I have to interpret the results? (I get different results for the two methods) Thank you again! Gabor Grothendieck wrote: > > Those are linear in the coefficients so try these: > > library(quantreg) > > rq1 <- rq(demand ~ Time + I(Time^2), data = BOD, tau= 1:3/4); rq1 > > # or > rq2 <- rq(demand ~ poly(Time, 2), data = BOD, tau = 1:3/4); rq2 > > > On Tue, Jun 9, 2009 at 10:55 AM, despaired wrote: >> >> Hi, >> >> I'm relatively new to R and need to do a quantile regression. Linear >> quantile regression works, but for my data I need some quadratic >> function. >> So I guess, I have to use a nonlinear quantile regression. I tried the >> example on the help page for nlrq with my data and it worked. But the >> example there was with a SSlogis model. Trying to write >> >> dat.nlrq <- nlrq(BM ~ I(Regen100^2), data=dat, tau=0.25, trace=TRUE) >> >> or >> >> dat.nlrq <- nlrq(BM ~ poly(Regen100^2), data=dat, tau=0.25, trace=TRUE) >> >> (I don't know the difference) both gave me the following error message: >> >> error in getInitial.default(func, data, mCall = as.list(match.call(func, >> : >> no 'getInitial' method found for "function" objects >> >> Looking in getInitial, it must have to do something with the starting >> parameters or selfStart model. But I have no idea, what this is and how I >> handle this problem. Can anyone please help? >> >> Thanks a lot in advance! >> -- >> View this message in context: >> http://www.nabble.com/Non-linear-regression-Quantile-regression-tp23944530p23944530.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. > > -- View this message in context: http://www.nabble.com/Non-linear-regression-Quantile-regression-tp23944530p23945900.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.
Re: [R] Non-linear regression/Quantile regression
Those are linear in the coefficients so try these: library(quantreg) rq1 <- rq(demand ~ Time + I(Time^2), data = BOD, tau= 1:3/4); rq1 # or rq2 <- rq(demand ~ poly(Time, 2), data = BOD, tau = 1:3/4); rq2 On Tue, Jun 9, 2009 at 10:55 AM, despaired wrote: > > Hi, > > I'm relatively new to R and need to do a quantile regression. Linear > quantile regression works, but for my data I need some quadratic function. > So I guess, I have to use a nonlinear quantile regression. I tried the > example on the help page for nlrq with my data and it worked. But the > example there was with a SSlogis model. Trying to write > > dat.nlrq <- nlrq(BM ~ I(Regen100^2), data=dat, tau=0.25, trace=TRUE) > > or > > dat.nlrq <- nlrq(BM ~ poly(Regen100^2), data=dat, tau=0.25, trace=TRUE) > > (I don't know the difference) both gave me the following error message: > > error in getInitial.default(func, data, mCall = as.list(match.call(func, : > no 'getInitial' method found for "function" objects > > Looking in getInitial, it must have to do something with the starting > parameters or selfStart model. But I have no idea, what this is and how I > handle this problem. Can anyone please help? > > Thanks a lot in advance! > -- > View this message in context: > http://www.nabble.com/Non-linear-regression-Quantile-regression-tp23944530p23944530.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.
Re: [R] Non-linear regression with latent variable
If you have the "RSiteSearch" package installed, you can do the
following:
library(RSiteSearch)
nrow(nll <- RSiteSearch.function("nonlinear regression with latent"))
HTML(nll)
This just produced 8 hits for me. If this doesn't solve your
problem, you might try other search terms. If all those fail, you
could write your own likelihood function and use the "maxLik" package.
Hope this helps.
Spencer
Samiul Hasan wrote:
Hi
Can anyone please suggest me a package where I can estimate a non-linear
regression model? One of the independent variables is latent or unobserved.
I have an indicator variable for this unobserved variable; however the
relationship is known to be non-linear also. In terms of equations my
problem is
y=f(latent, fixed)
q=g(latent) where q is the indicator variable
For me both f and g are non-linear.
Thanks
Samiul Hasan
__
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.
Re: [R] non linear regression with nls
Thank you! It worked perfectly, also for the other variables! Messaggio originale Da: r...@life.ku.dk Data: 06.02.2009 13.29 A: Oggetto: Re: [R] non linear regression with nls Hi Laura, I think you have to make a list formulas: formList <- list(NT.N ~ fz1(Portata, a, b), NT.N ~ fz2(Portata, a, b), NT.N ~ fz3(Portata, a, b, d, e), NT.N ~ fz4(Portata, a, b), NT.N ~ fz5(Portata, a, b, d)) and then in the loop: resultList[[i]] <- nls(formList[[i]], ... Christian __ 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.
Re: [R] non linear regression with nls
Hello, thanks for the advice of nlsList!
I tried to look at the help page of nlsList, but I didnt understand how to use
the subset argument of the function
and it's not clear to
me if this only allows you to choose one subset or if it run the regression for
every given subset, in this case how
can someone specify the groups/subset?
Thanks a lot!
Laura
Messaggio originale
Da: kfr...@wisc.edu
Data: 03.02.2009 15.36
A:
Oggetto: Re: [R] non linear regression with nls
Hi, Laura-
You might have a look at ?nlsList().
Ken
- Original Message -
From: "lauramorg...@bluewin.ch"
Date: Tuesday, February 3, 2009 4:38 am
Subject: [R] non linear regression with nls
To: r-help@r-project.org
> Hello,
> I'm a beginner with R and it's the first time I'm using the R-help
> list... I hope I'm in the right place, if not:
> Sorry!!
>
> I need to do non linear regressions on a data set which columns are:
> "river.name""Portata" "PTG.P" "PO4.P" "NT.N""NH4.N"
> "NO3.N" "BOD5""SiO2"
> I need to predict every variable (PTG, PO4, NT, ..., which are
> concentration of substances in water) starting from
> the "Portata" variable (which is the water flow)
> The functions that I'm using are:
> > fz1<-function(Portata, a, b){a+b/Portata}
> > fz2<-function(Portata, a, b){a*exp(b*Portata)}
> > fz3<-function(Portata, a, b, d, e){a+b/Portata+d*(Portata^e)}
> > fz4<-function(Portata, b, d){b*Portata^d}
> > fz5<-function(Portata, a, b, d){a+b*(Portata^d)}
> I've made a list of the functions with list(fz1, fz2, fz3, fz4, fz5)
>
> and the starting , lower and upper parameters for each function are:
> fz1: start=list(a=10, b=10), lower=list(a=0, b=0), upper=list(a=1000,
> b=1000)
> fz2: start=list(a=10, b=1), lower=list(a=0, b=0), upper=list(a=1000, b=1000)
> fz3: start=list(a=10, b=10, d=10, e=10), lower=list(a=0, b=0, d=0,
> e=-50), upper=list(a=1000, b=1000, d=1000, e=50)
> fz4: start=list (a=10, b=1), lower=list(a=0, b=-50),
> upper=list(a=1000, b=50)
> fz5: start=list(a=10, b=10, d=1), lower=list(a=0, b=0, d=-50),
> upper=list(a=1000, b=1000, d=50)
>
> so far i manage to do non linear regression one at a time that is,
> using one function for one river) using nls(), for
> example:
> r.NT.lav<-nls(NT.N~ fz1(Portata, a,b),
> data=subset(dati,river.name=="Laveggio"), start=list(a=10, b=10), nls.control
> (maxiter=200), algorithm='port', trace=TRUE, na.action=na.omit,
> lower=list(a=0, b=0), upper=list(a=1000, b=1000))
> and then I get the results with summary() and the graph using plotfit()
>
> but this will get extremly long since I have 12 rivers to analize for
> every variable and then compare the results, so
> I'd like to use a loop (cycle for??) but I can't figure out how it
> works. I've tried to look on the help page on ?
> Control (control flow) but I didn't understand it...
> Can somebody help me (give me a hint or an example of a loop) to
> automize the regression and save the results
> Please consider that my knowledge of computer programming is
> practically non-existent!!
> Thanks a lot!
> Laura F.
>
> __
> 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.
Re: [R] non linear regression with nls
Thank you a lot Mr. Ritz!
I've tried the loop you suggested and I added a list for lower and upper limits
of the parameters, but there is still
a problem...
I have a list of functions, which works...
fz1<-function(Portata, a, b){a+(b/Portata)}
fz2<-function(Portata, a, b){a*exp(b*Portata)}
fz3<-function(Portata, a, b, d, e){a+(b/Portata)+d*(Portata^e)}
fz4<-function(Portata, a, b){a*Portata^b}
fz5<-function(Portata, a, b, d){a+b*(Portata^d)}
fctList <- list(fz1, fz2, fz3, fz4, fz5)
as well as lists for starting values, upper and lower values, which work as
well:
startList <- list(list(a=10, b=10), list(a=10, b=1), list(a=10, b=10, d=10,
e=1), list(a=10, b=1), list(a=10, b=10,
d=1))
lowerList<-list(list(a=0,b=0),list(a=0,b=0), list(a=0,b=0,d=0,e=-50),list(a=0,
b=-50), list(a=0, b=0, d=-50))
upperList<-list(list(a=1000, b=1000), list(a=1000, b=1000),
list(a=1000,b=1000,d=1000,e=50), list(a=1000,b=50), list
(a=1000, b=1000, d=50))
but if I try to run this for loop
resultList <- list()
for (i in 1:5)
{
resultList[[i]] <- nls(NT.N ~ fctList[[i]](Portata, a,b), data=subset(dati,
Fiume=="Laveggio"), start=startList[[i]],
nls.control(maxiter=200), algorithm='port', trace=TRUE, na.action=na.omit,
upper=upperList[[i]], lower=lowerList[[i]])
}
I get the following error message:
"Error in fctList[[i]](Portata, a, b) : element 1 is empty;
the part of the args list of '*' being evaluated was:
(d, (Portata^e))"
I realized that the problem is the element after the function, i.e. (Portata,
a, b), since fct 3 and 5 have more
parameters (Portata,a,b,d,e).
So I tried to make a list (parList) for this too, i tried 2 versions:
#version1
list("Portata","a", "b")->pf1.2.4
list("Portata","a","b","d","e")->pf3
list("Portata","a","b","d")->pf5
parList<-list(pf1.2.4, pf1.2.4,pf3,pf1.2.4,pf5)
#version 2
parList<-list("Portata,a,b","Portata,a,b","Portata,a,b,d,e",
"Portata,a,b","Portata,a,b,d")
and then I tried them (one at a time) in the loop:
resultList <- list()
for (i in 1:5)
{
resultList[[i]] <- nls(NT.N ~ fctList[[i]](parList[[i]]), data=subset(dati,
Fiume=="Laveggio"), start=startList
[[i]], nls.control(maxiter=200), algorithm='port', trace=TRUE,
na.action=na.omit,upper=upperList[[i]], lower=lowerList
[[i]])
}
but I got this error message: " Error in fctList[[i]](parList[i]) : element 1
is empty;
the part of the args list of '+' being evaluated was:
(a, (b/Portata))"
What can I do to fix it?
I'm also wondering which kind of function (maybe another loop?) I could use to
automize the regression not only for
the variable NT.N but for every variable (PTG.P, PO4.P,. ..)
My dataframe look like this (a sample):
Fiume giorno mese anno Portata PTG.P PO4.P NT.N NH4.N NO3.N
BOD5 SiO2 data
1Vedeggio 101 1995 0.981 218.40 118.000 9.196 6.5700 2.06
6.080 4.33 34709
2Vedeggio 72 1995 0.965 125.84 54.000 8.701 5.2600 2.31
16.480 4.43 34737
3Vedeggio 73 1995 1.536 37.44 12.000 7.271 5.5600 1.88
5.240 4.15 34765
...
190 Cassarate 299 2008 1.240 26.00 20.000 2.480 0.1200 1.79
1.700 4.03 39720
191 Cassarate 13 10 2008 0.860 23.00 16.000 2.720 0.0200 2.13
1.780 3.71 39734
192 Cassarate 10 11 2008 8.840 26.00 14.000 2.900 0.0500NA
1.400 3.62 39762
193 Cassarate 9 12 2008 2.030 35.00 23.000 2.190 0.0700 1.79
1.950 3.74 39791
...
279 Laveggio 151 2002 0.347 77.00 30.000 9.690 0.4300 7.23
1.950 4.17 37271
280 Laveggio 112 2002 0.527 54.00 17.000 7.520 0.8800 5.87
2.410 3.58 37298
281 Laveggio 133 2002 0.900 34.00 15.000 7.520 0.7100 6.17
6.550 3.03 37328
...
Thanks to anyone that could give me any hint!!
Laura
## a for loop
resultList <- list()
for (i in 1:5)
{
## storing the result as the i'th list component
## notice that the i'th list components in fctList and startList
## are used for the i'th fit
resultList[[i]] <- nls(NT.N ~ fctList[[i]](parList[[i]]), data=subset(dati,
Fiume=="Laveggio"), start=startList
[[i]], nls.control(maxiter=200), algorithm='port', trace=TRUE,
na.action=na.omit)
}
Messaggio originale
Da: r...@life.ku.dk
Data: 03.02.2009 19.00
A:
Oggetto: Re: [R] non linear regression with nls
Hi Laura,
I've the following suggestion for you using several lists and a "for" loop:
fz1<-function(Portata, a, b){a+b/Portata}
fz2<-function(Portata, a, b){a*exp(b*Portata)}
fz3<-function(Portata, a, b, d, e){a+b/Portata+d*(Portata^e)}
fz4<-function(Portata, b, d){b*Portata^d}
fz5<-function(Portata, a,
Re: [R] Non linear regression with 2 explanatory variables
Dear Rolf,
One thing that sometimes makes nls easier to apply is using the 'formula'
argument like you would use the 'fn' argument of optim. That is, if you
have a residual function that has arguments x, y, a, b and you need to
optimize a and b, you would make a call like
nls(~resid(x,y,a=astart, b=bstart), control = nls.control(warnOnly =
TRUE, printEval = TRUE), start = list(a=astart, b=bstart))
This did not work easily before R-2.6.0, but does now. The Puromycin
analysis from the help files is an example of this useage and below is
another.
Or do you already use nls this way and still have problems?
# get data as a sum of exponentials
dataSumOfExp <- function(rates = seq(.05, .005, length=3),
times = 1:100,
amps = rep(1, length(rates))) {
tfun <- function(t,r) exp(-r*t)
## get C with tfun
C <- mapply(tfun, r=rates, MoreArgs=list(t=times))
## add the columns of C with relative amplitudes 1, and add noise
C %*% amps + rnorm( nrow(C) ) * max(C) * .1
}
# residual function
resFun <- function(rates, amps, measured, times = 1:100) {
tfun <- function(t,r) exp(-r*t)
CEst <- mapply(tfun, r=rates, MoreArgs=list(t=times))
measured - CEst %*% amps
}
# get data
measured <- dataSumOfExp()
# optimize rates of exponentials and their relative amplitudes
res <- nls(~resFun(rates = rates, measured = measured, amps = amps),
control = nls.control(warnOnly = TRUE, printEval = TRUE),
start = list(rates = c(.04, .1, .001),
amps = rep(1,3)), trace = TRUE)
summary(res)
On Thu, 17 Jan 2008, Rolf Turner wrote:
>
> I have never had much success in using nls(). If you scan the archives
> you will find one or two postings from me on this topic. I have
> received
> no useful responses to these postings.
>
> I have found that anything that I tried (and failed) to do using nls()
> could be done quite easily using optim().
>
> cheers,
>
> Rolf Turner
>
>
> On 17/01/2008, at 3:56 AM, Gavin Simpson wrote:
>
> > hits=-2.6 tests=BAYES_00
> > X-USF-Spam-Flag: NO
> >
> > On Wed, 2008-01-16 at 11:02 +0100, Janice Kielbassa wrote:
> >> Hello!
> >>
> >> I want to do a non-linear regression with 2 explanatory variables
> >> (something like : length ~ a * time * exp( b* temperature)), having a
> >> data set (length, time, temperature). Which function could I use (I
> >> tried nls but I think it doesn't work)
> >
> > Janice, I'll start by saying I can't help you as I have never used
> > nls()
> > myself and I am not familiar with this type of analysis.
> >
> > Why do you think that nls() "doesn't work"? It is a widely used
> > part of
> > R and thus probably very well tested.
> >
> > My understanding of these things is that nls is a sophisticated tool
> > that requires some effort on the part of the user, such as selecting
> > appropriate starting values.
>
>
> ##
> Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
>
> __
> 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.
Re: [R] Non linear regression with 2 explanatory variables
I have never had much success in using nls(). If you scan the archives
you will find one or two postings from me on this topic. I have
received
no useful responses to these postings.
I have found that anything that I tried (and failed) to do using nls()
could be done quite easily using optim().
cheers,
Rolf Turner
On 17/01/2008, at 3:56 AM, Gavin Simpson wrote:
> hits=-2.6 tests=BAYES_00
> X-USF-Spam-Flag: NO
>
> On Wed, 2008-01-16 at 11:02 +0100, Janice Kielbassa wrote:
>> Hello!
>>
>> I want to do a non-linear regression with 2 explanatory variables
>> (something like : length ~ a * time * exp( b* temperature)), having a
>> data set (length, time, temperature). Which function could I use (I
>> tried nls but I think it doesn't work)
>
> Janice, I'll start by saying I can't help you as I have never used
> nls()
> myself and I am not familiar with this type of analysis.
>
> Why do you think that nls() "doesn't work"? It is a widely used
> part of
> R and thus probably very well tested.
>
> My understanding of these things is that nls is a sophisticated tool
> that requires some effort on the part of the user, such as selecting
> appropriate starting values.
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
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.
Re: [R] Non linear regression with 2 explanatory variables
Gavin Simpson wrote: > hits=-2.6 tests=BAYES_00 > X-USF-Spam-Flag: NO > > On Wed, 2008-01-16 at 11:02 +0100, Janice Kielbassa wrote: >> Hello! >> >> I want to do a non-linear regression with 2 explanatory variables >> (something like : length ~ a * time * exp( b* temperature)), having a >> data set (length, time, temperature). Which function could I use (I >> tried nls but I think it doesn't work) > > Janice, I'll start by saying I can't help you as I have never used nls() > myself and I am not familiar with this type of analysis. > maybe it helps if you have a look at Chapter 10 "Nonlinear Models" by Douglas M. Bates and John M. Chambers in: John M. Chambers, Trevor J. Hastie (Eds.): Statistical Models in S. Chapman & Hall/CRC , 1992 Best, Roland __ 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.
Re: [R] Non linear regression with 2 explanatory variables
hits=-2.6 tests=BAYES_00 X-USF-Spam-Flag: NO On Wed, 2008-01-16 at 11:02 +0100, Janice Kielbassa wrote: > Hello! > > I want to do a non-linear regression with 2 explanatory variables > (something like : length ~ a * time * exp( b* temperature)), having a > data set (length, time, temperature). Which function could I use (I > tried nls but I think it doesn't work) Janice, I'll start by saying I can't help you as I have never used nls() myself and I am not familiar with this type of analysis. Why do you think that nls() "doesn't work"? It is a widely used part of R and thus probably very well tested. My understanding of these things is that nls is a sophisticated tool that requires some effort on the part of the user, such as selecting appropriate starting values. You are unlikely to get any further assistance from the list unless you read the posting guide and post an example of what you did (preferably with the actual data or dummy data with the same properties if not) and the exact error message or output from R that lead you to believe that nls() did not work. HTH G > Thanks a lot! > > Janice > > __ > 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. -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr. Gavin Simpson [t] +44 (0)20 7679 0522 ECRC, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ 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.
