How to Calculate Total Least Squares in R? (Orthogonal Regression)

ggplot2: How to plot an orthogonal regression line?

It turns out that you can extract the slope and intercept from principal components analysis on (x,y), as shown here. This is just a little simpler, runs in base R, and gives the identical result to using Deming(...) in MethComp.

# same `x and `y` as @user20650's answer
df  <- data.frame(y, x)
pca <- prcomp(~x+y, df)
slp <- with(pca, rotation[2,1] / rotation[1,1])
int <- with(pca, center[2] - slp*center[1])

ggplot(df, aes(x,y)) + 
  geom_point() + 
  stat_smooth(method=lm, color="green", se=FALSE) +
  geom_abline(slope=slp, intercept=int, color="blue")

Nonlinear total least squares/Deming regression

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")

Adding orthogonal regression line in ggplot

Borrowed from this blog post & this answer. Basically, you will need Deming function from MethComp or prcomp from stats packages together with a custom function perp.segment.coord. Below is an example taken from above mentioned blog post.

library(ggplot2)
library(MethComp)

data(airquality)
airquality <- na.exclude(airquality)

# Orthogonal, total least squares or Deming regression
deming <- Deming(y=airquality$Wind, x=airquality$Temp)[1:2]
deming  
#>  Intercept      Slope 
#> 24.8083259 -0.1906826

# Check with prcomp {stats}
r <- prcomp( ~ airquality$Temp + airquality$Wind )
slope <- r$rotation[2,1] / r$rotation[1,1]
slope   
#> [1] -0.1906826

intercept <- r$center[2] - slope*r$center[1]
intercept
#> airquality$Wind 
#>        24.80833

# https://stackoverflow.com/a/30399576/786542
perp.segment.coord <- function(x0, y0, ortho){
  # finds endpoint for a perpendicular segment from the point (x0,y0) to the line
  # defined by ortho as y = a + b*x
  a <- ortho[1]  # intercept
  b <- ortho[2]  # slope
  x1 <- (x0 + b*y0 - a*b)/(1 + b^2)
  y1 <- a + b*x1
  list(x0=x0, y0=y0, x1=x1, y1=y1)
}

perp.segment <- perp.segment.coord(airquality$Temp, airquality$Wind, deming)
perp.segment <- as.data.frame(perp.segment)
    
# plot
plot.y <- ggplot(data = airquality, aes(x = Temp, y = Wind)) + 
  geom_point() +
  geom_abline(intercept = deming[1],
              slope = deming[2]) +
  geom_segment(data = perp.segment, 
               aes(x = x0, y = y0, xend = x1, yend = y1), 
               colour = "blue") +
  theme_bw()

Created on 2018-03-19 by the reprex package (v0.2.0).

{Methcomp} – Deming / orthogonal regression – goodness of fit + confidence intervals

There are many proposed methods to calculate goodness of fit and tolerance intervals for Deming Regression but none of them widely accepted. The conventional methods we use for OLS regression may not make sense. This is an area of active research. I don't think there many R-packages which will help you compute that since not many mathematicians agree on any particular method. Most methods for calculating intervals are based on Resampling techniques.

However you can check out the 'mcr' package for intervals...
https://cran.r-project.org/web/packages/mcr/

Orthogonal Distance Regression in 3d

Try http://tolstoy.newcastle.edu.au/R/help/05/07/8039.html

Apparently there's a package, prcomp, which may meet your needs.

---

Related Topics