3
$\begingroup$

Disclaimer: Statistics is not my strong side, so if my question is nonsense I apologize. I'm a beginner, but really wanting to understand this.

My question is: why do I get so widely different parameter estimates when using different transformations on my data in a non-linear regression ?

I'm trying to do a nonlinear regression and to estimate the uncertainty of the fit (confidence interval) using linear approximation. From my understanding the more linear-like the shape of the nonlinear function, the more accurate will the confidence interval calculation by linear approximation be. I therefore want to transform the data to make it as linear as possible. The errors in $y$ can be assumed to be log-normal. My data is monotonic and assumed to follow a power function in most cases.

$$ y = a*(x-x_0)^b $$

where $y$ is river discharge, $x$ is an arbitrary water level in the river and $x_0$ is the water level where where discharge $y$ is 0. This can be rewritten as log transformed, and nice and linear $$ log(y) = a + b \times log(x-x_0) $$.

I need to estimate the parameters $a$, $b$ and $x_0$, so to do so simultaneously I use nonlinear regression. I also have some data that follows quadratic functions, so I would like to set up (and understand) a non-linear method.

I use r and nlsLM() from minpack.lm to carry out the non-linear regression. Here is some example code:

library(minpack.lm)

xdata <- c(19,  21,  24,    25, 29, 34, 35, 40, 40, 46, 48, 48, 52, 56, 57, 65, 65, 68)
ydata <- c(10,  11, 14, 20, 24, 50, 42, 96, 89, 134,    135,    161,    171,    218,    261,    371,    347,    393)
df<-data.frame(x=xdata, y=ydata)

#weights applied in the case of no transformation (relative error assumed to be the same for all y data)
W<-1/ydata

# NLS regression with weights, no transformation
nlsmodel1<-nlsLM(y ~ a*(x-x0)^b,data=df,start=list(a=0.1, b=2.5,x0=0))

# log transformed
nlsmodel2<-nlsLM(log(y) ~ a+b*(log(x-x0)),data=df,start=list(a=0.1, b=2.5,x0=0))
> coef(nlsmodel1)
          a           b          x0 
0.005158377 2.719693093 4.896772931 
> coef(nlsmodel2)
        a         b        x0 
-8.683758  3.445699 -4.139127 

> exp(-8.683758)
[1] 0.0001693136

I understand that the weights are very important, and can have a say in the differences here, but not by this much? My judgement of the two parameter sets is that nlsmodel1 performs "better", and that the b coefficient is too high in the fit from nlsmodel2. nlsmodel2 does a poor job in the upper end of the data, with large residuals there. But why are they so different? I feel like I'm doing something very silly here, and is unable to see the error. I have tried some other transformations, for example only transforming LHS as log(y), but the problem remains.

I appreciate any tips that can help me improve, and not the least understand, the transformed fit.

Cheers

Related post #1 and post #2

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ The differences in $b$ and $x_0$ do not appear large: check their standard errors. The difference in $a$ is due to a basic mistake with logs: $\log(y) = \log(a) + b\log(x-x_0)$ (which I guess is compensated for in the last line), whence the "a" in the log model should be close to the log of the "a" in the original model. When this change is made (by replacing "a" in the first model with "exp(a)") the changes in all three coefficients are within the ranges indicated by their standard errors. In short, all appears fine: is there really a problem to be resolved here? $\endgroup$
    – whuber ♦
    Commented Dec 3, 2013 at 20:04
  • 1
    $\begingroup$ Thanks @whuber! Cheers for seeing the log(a) error. I do optimize it as a nontransformed parameter, as it is prone to creating NaN's. And after thinking some more about this, I think you are right and there is no problem. I was surprised it gave such different exponent values, but when looking closer at the $x_0$ differences and how the two solutions optimize it, it does make sense. Glad you made me see it this way! $\endgroup$
    – rhkarls
    Commented Dec 4, 2013 at 15:32
  • 2
    $\begingroup$ I have a couple of practical suggestions then for achieving stable and reliable solutions. (1) Use $\exp(a)$ in model 1 and $a$ in model 2 to obtain comparable results without NaN problems. (2) Run model 2 first and feed its estimates as starting values to model 1. (3) Break model 2 into two parts: only $x_0$ is nonlinear; $a$ and $b$ are efficiently determined by least squares (e.g., use lm). You can then use a one variable application of nls to estimate $x_0$. If you want the full output of nlsLM, feed it those estimates as starting values. $\endgroup$
    – whuber ♦
    Commented Dec 4, 2013 at 15:41
  • $\begingroup$ Thats a great tip on the starting values. The parameters are very strongly correlated, and I guess the solution is prone to local minima as well. The reason I apply a nonlinear solution is to estimate the 3 parameters simultaneously to quantify the correlation and hence get a more realistic uncertainty estimate. Now I just need to come up with a good way to do the same with a 3 degree polynomial with zero intercept. Any good tips on that @whuber ? $\endgroup$
    – rhkarls
    Commented Dec 5, 2013 at 14:58
  • 2
    $\begingroup$ Offhand I would consider orthonormalizing the basis of polynomial functions $x^i,$ $i=1,2,3$ (as evaluated at the data) and using ordinary least squares. $\endgroup$
    – whuber ♦
    Commented Dec 5, 2013 at 16:28

2 Answers 2

Reset to default
2
$\begingroup$

One note about weights: It is better to weight by the predicted value of Y (for given X) rather than the observed value of Y. The observed values can vary for random reasons, and this can (in some cases) throw the weights way off, and lead to suboptimal results.

I don't know enough about R to know if it can handle weighting by the reciprocal of the predicted Y value squared (rather than the reciprocal of the observed Y value squared). You could fit with the free demo of GraphPad Prism (which weights by predicted Y) and see how different the results are. It depends, of course, on how scattered your data are -- how far apart the predicted and observed Y values tend to be.

Improve this answer
$\endgroup$
0
$\begingroup$

I am not a statistician either but, as a physicist, I use to face quite often the same kind of problems. So, I shall give you here my experience.

If your original model is y = a (x - x0)^b and you want to weight using w =1/y, what I suggest is to run the linearized model without any weight first. From there, you have the parameters from log(a) and b; now, run the nonlinear regression with the weights.

If your weights were w = 1/y^2, corresponding to the minimization of the sum of squares of relkative errors, you should notice, via a simple Taylor expansion, that, if errors are not large, the sum of squares for the linearized model (no weight) is "similar" to the sum of squares of the non linear model.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.