Nonlinear total least squares/Deming regression

Question

I've been using nls() to fit a custom model to my data, but I don't like how the model is fitting and I would like to use an approach that minimizes residuals in both x and y axes.

I've done a lot of searching, and have found solutions for fitting linear models such as the deming package (http://cran.r-project.org/web/packages/deming/index.html) and these stackoverflow posts: Total Least Square method using R, How to calculate Total least squares in R? (Orthogonal regression). I've also found matlab solutions (https://stats.stackexchange.com/questions/110772/total-least-squares-curve-fit-problem), but these fit a second order polynomial and not a custom, user-defined model.

What I would like is something similar to nls() that does the x and y residual minimization. This would allow me to enter my custom model. Is anyone aware of any solution in R?

Many thanks!

EDIT

Here's an example, but please note that I'm seeking suggestions on a general solution for nonlinear total least squares regression, and not something specific to this dataset (this is just an example data based on Modifying a curve to prevent singular gradient matrix at initial parameter estimates):

df <- structure(list(x = c(3, 4, 5, 6, 7, 8, 9, 10, 11), y = c(1.0385, 
1.0195, 1.0176, 1.01, 1.009, 1.0079, 1.0068, 1.0099, 1.0038)), .Names = c("x", 
"y"), row.names = c(NA, -9L), class = "data.frame")

(nlsfit <- nls(y ~ a^b^x, data = df, start = c(a=0.9, b=0.6)))

library(ggplot2)
ggplot(df, aes(x=x, y=y)) + 
    geom_point() + 
    geom_smooth(method="nls", formula = y ~ a^b^x, se=F, start = list(a=0.9, b=0.6))

It wasn't me that down voted it but the question should at least show the code and the input data for what you have tried. — G. Grothendieck
thanks G. Grothendieck! I've added an example for nonlinear regression. Still wondering if anyone knows of a nonlinear TLS analogue. thnaks! — Thomas
maybe PCA regression? principal components minimize the orthogonal distance from a fitted line to any given point. cerebralmastication.com/wp-content/uploads/2010/09/pca.png . i think you were probably downvoted because your question is more appropriate for cross validated than stackoverflow (stats.stackexchange.com) — bjoseph
Thanks bjoseph, I appreciate the feedback. will post my question there. Also, I did came across PCA regression during my search but wonder if that is only for linear case (not nonlinear?). thank you! — Thomas
Its not completely clear how to know if a particular answer is what you want but it does seem that the first point has high leverage and so perhaps that is distorting the fit relative to what you want. If its true in general that your curves have high leverage points you could use a robust fitting method: library(robustbase); nlrob(y ~ a^b^x, df, start = coef(nlsfit), method = "tau", lower = c(a = 0, b = 0), upper = 2) . For orthogonal least squares you could try the onls package but I am not so sure that would give a fit much different than that shown in the question. — G. Grothendieck

Thomas Thomas · Accepted Answer · 2015-03-15T02:21:45

G. Grothendieck and Brian Borchers at CrossValidated suggested the onls package, which is exactly what I was looking for. Thanks everyone, the help is much appreciated. see here for more info: http://www.r-bloggers.com/introducing-orthogonal-nonlinear-least-squares-regression-in-r/

Here's an example using the same data from above - this gives the same fitted parameters as the regular nls(), however it does make a difference in my real data. But this at least shows how carry out the operations.

df <- structure(list(x = c(3, 4, 5, 6, 7, 8, 9, 10, 11), y = c(1.0385, 
1.0195, 1.0176, 1.01, 1.009, 1.0079, 1.0068, 1.0099, 1.0038)), .Names = c("x", 
"y"), row.names = c(NA, -9L), class = "data.frame")

library(onls)
(onlsfit <- onls(y ~ a^b^x, data = df, start = c(a=0.9, b=0.6)))

# define function containing onls fitted parameters for plotting
fit.model <- function(x) {1.0934^0.7242^x}

library(ggplot2)
ggplot(df, aes(x=x, y=y)) + 
    geom_point() + 
    stat_function(fun=fit.model, color="black")

Nonlinear total least squares/Deming regression

1 Answers