12
$\begingroup$

I made $n$ measurements of two variables $x$ and $y$. They both have known uncertainties $\sigma_x$ and $\sigma_y$ associated with them. I want to find the relation between $x$ and $y$. How can I do it?

EDIT: each $x_i$ has a different $\sigma_{x,i}$ associated with it, and the same with the $y_i$.


Reproducible R example:

## pick some real x and y values 
true_x <- 1:100
true_y <- 2*true_x+1

## pick the uncertainty on them
sigma_x <- runif(length(true_x), 1, 10) # 10
sigma_y <- runif(length(true_y), 1, 15) # 15

## perturb both x and y with noise 
noisy_x <- rnorm(length(true_x), true_x, sigma_x)
noisy_y <- rnorm(length(true_y), true_y, sigma_y)

## make a plot 
plot(NA, xlab="x", ylab="y",
    xlim=range(noisy_x-sigma_x, noisy_x+sigma_x), 
    ylim=range(noisy_y-sigma_y, noisy_y+sigma_y))
arrows(noisy_x, noisy_y-sigma_y, 
       noisy_x, noisy_y+sigma_y, 
       length=0, angle=90, code=3, col="darkgray")
arrows(noisy_x-sigma_x, noisy_y,
       noisy_x+sigma_x, noisy_y,
       length=0, angle=90, code=3, col="darkgray")
points(noisy_y ~ noisy_x)

## fit a line 
mdl <- lm(noisy_y ~ noisy_x)
abline(mdl)

## show confidence interval around line 
newXs <- seq(-100, 200, 1)
prd <- predict(mdl, newdata=data.frame(noisy_x=newXs), 
    interval=c('confidence'), level=0.99, type='response')
lines(newXs, prd[,2], col='black', lty=3)
lines(newXs, prd[,3], col='black', lty=3)

linear regression without considering errors in variables

The problem with this example is that I think it assumes that there are no uncertainties in $x$. How can I fix this?

$\endgroup$
  • $\begingroup$ True, lm fits a linear regression model, that is: a model of the expectation of $Y$ with respect to $P(Y | X)$, in which clearly $Y$ is as random and $X$ is considered known. To deal with uncertainty in $X$ you will need a different model. $\endgroup$ – conjugateprior Mar 15 '16 at 17:59
  • 1
    $\begingroup$ For your rather special case (univariate with a known ratio of noise levels for X and Y) Deming regression will do the trick, e.g. the Deming function in R package MethComp. $\endgroup$ – conjugateprior Mar 15 '16 at 18:08
  • 1
    $\begingroup$ @conjugateprior Thanks, this looks promising. I'm wondering: does Deming regression still work if I have a different (but still known) variance on each individual x and y? i.e. if the x's are lengths, and I used rulers with different precisions to obtain each x $\endgroup$ – rhombidodecahedron Mar 15 '16 at 20:57
  • $\begingroup$ I think perhaps the way to solve it when there are different variances for each measurement is using York's method. Does anyone happen to know if there's an R implementation of this method? $\endgroup$ – rhombidodecahedron Mar 15 '16 at 21:50
  • 1
    $\begingroup$ @rhombidodecahedron See the "with measured errors" fit in my answer there: stats.stackexchange.com/questions/174533/… (which was taken from the documentation of package deming). $\endgroup$ – Roland Mar 16 '16 at 9:47
9
$\begingroup$

Let the true line $L$, given by an angle $\theta$ and a value $\gamma$, be the set

$$(x,y): \cos(\theta) x + \sin(\theta) y = \gamma.$$

The signed distance between any point $(x,y)$ and this line is

$$d(x,y;L) = \cos(\theta) x + \sin(\theta) y - \gamma.$$

Letting the variance of $x_i$ be $\sigma_i^2$ and that of $y_i$ be $\tau_i^2$, independence of $x_i$ and $y_i$ implies the variance of this distance is

$$\operatorname{Var}(d(x_i,y_i;L)) = \cos^2(\theta)\sigma_i^2 + \sin^2(\theta)\tau_i^2.$$

Let us therefore find $\theta$ and $\gamma$ for which the inverse variance weighted sum of squared distances is as small as possible: it will be the maximum likelihood solution if we assume the errors have bivariate normal distributions. This requires a numerical solution, but it's straightforward to find a with a few Newton-Raphson steps beginning with a value suggested by an ordinary least-squares fit.

Simulations suggest this solution is good even with small amounts of data and relatively large values of $\sigma_i$ and $\tau_i$. You can, of course, obtain standard errors for the parameters in the usual ways. If you're interested in the standard error of the position of the line, as well as the slope, then you might wish first to center both variables at $0$: that should eliminate almost all the correlation between the estimates of the two parameters.


The method works so well with the example of the question that the fitted line is almost distinguishable from the true line in the plot: they are within one unit or so of each other everywhere. Instead, in this example the $\tau_i$ are drawn iid from an exponential distribution and the $\sigma_i$ are drawn iid from an exponential distribution with twice the scale (so that most of the error tends to occur in the $x$ coordinate). There are only $n=8$ points, a small number. The true points are equally spaced along the line with unit spacing. This is a fairly severe test, because the potential errors are noticeable compared to the range of the points.

Figure

The true line is shown in dotted blue. Along it the original points are plotted as hollow circles. Gray arrows connect them to the observed points, plotted as solid black disks. The solution is drawn as a solid red line. Despite the presence of large deviations between observed and actual values, the solution is remarkably close to the correct line within this region.

#
# Generate data.
#
theta <- c(1, -2, 3) # The line is theta %*% c(x,y,-1) == 0
theta[-3] <- theta[-3]/sqrt(crossprod(theta[-3]))
n <- 8
set.seed(17)
sigma <- rexp(n, 1/2)
tau <- rexp(n, 1)
u <- 1:n
xy.0 <- t(outer(c(-theta[2], theta[1]), 0:(n-1)) + c(theta[3]/theta[1], 0))
xy <- xy.0 + cbind(rnorm(n, sd=sigma), rnorm(n, sd=tau))
#
# Fit a line.
#
x <- xy[, 1]
y <- xy[, 2]
f <- function(phi) { # Negative log likelihood, up to an additive constant
  a <- phi[1]
  gamma <- phi[2]
  sum((x*cos(a) + y*sin(a) - gamma)^2 / ((sigma*cos(a))^2 + (tau*sin(a))^2))/2
}
fit <- lm(y ~ x) # Yields starting estimates
slope <- coef(fit)[2]
theta.0 <- atan2(1, -slope)
gamma.0 <- coef(fit)[1] / sqrt(1 + slope^2)
sol <- nlm(f,c(theta.0, gamma.0))
#
# Plot the data and the fit.
#
theta.hat <- sol$estimate[1] %% (2*pi)
gamma.hat <- sol$estimate[2]
plot(rbind(xy.0, xy), type="n", xlab="x", ylab="y")
invisible(sapply(1:n, function(i) 
  arrows(xy.0[i,1], xy.0[i,2], xy[i,1], xy[i,2], 
         length=0.15, angle=20, col="Gray")))
points(xy.0)
points(xy, pch=16)
abline(c(theta[3] / theta[2], -theta[1]/theta[2]), col="Blue", lwd=2, lty=3)
abline(c(gamma.hat / sin(theta.hat), -1/tan(theta.hat)), col="Red", lwd=2)
$\endgroup$
  • $\begingroup$ +1. As far as I understand, this answers this older Q too: stats.stackexchange.com/questions/178727? We should close it as a duplicate then. $\endgroup$ – amoeba says Reinstate Monica Mar 15 '16 at 23:22
  • $\begingroup$ Also, as per my comment to the answer in that thread, it looks like deming function can handle variable errors too. It should probably yield a fit very similar to yours. $\endgroup$ – amoeba says Reinstate Monica Mar 15 '16 at 23:26
  • $\begingroup$ I wonder if the flow of the discussion makes more sense if you switch the places of the 2 paragraphs above & below the figure? $\endgroup$ – gung - Reinstate Monica Mar 16 '16 at 0:42
  • 3
    $\begingroup$ I was reminded this morning (by a voter) that this question had been asked and answered in multiple ways, with working code, several years ago on the Mathematica SE site. $\endgroup$ – whuber Mar 16 '16 at 14:20
  • $\begingroup$ Does this solution have a name? and possibly a resource for further reading (besides the Mathematica SE site i mean)? $\endgroup$ – JustGettinStarted Sep 17 '18 at 14:14
0
$\begingroup$

The maximum likelihood optimization for the case of uncertainties in x and y has been addressed by York (2004). Here is R code for his function.

"YorkFit", written by Rick Wehr, 2011, translated into R by Rachel Chang

Universal routine for finding the best straight line fit to data with variable, correlated errors, including error and goodness of fit estimates, following Eq. (13) of York 2004, American Journal of Physics, which was based in turn on York 1969, Earth and Planetary Sciences Letters

YorkFit <- function(X,Y, Xstd, Ystd, Ri=0, b0=0, printCoefs=0, makeLine=0,eps=1e-7)

X, Y, Xstd, Ystd: waves containing X points, Y points, and their standard deviations

WARNING: Xstd and Ystd cannot be zero as this will cause Xw or Yw to be NaN. Use a very small value instead.

Ri: correlation coefficients for X and Y errors -- length 1 or length of X and Y

b0: rough initial guess for the slope (can be gotten from a standard least-squares fit without errors)

printCoefs: set equal to 1 to display results in the command window

makeLine: set equal to 1 to generate a Y wave for the fit line

Returns a matrix with the intercept and slope plus their uncertainties

If no initial guess for b0 is provided, then just use OLS if (b0 == 0) {b0 = lm(Y~X)$coefficients[2]}

tol = abs(b0)*eps #the fit will stop iterating when the slope converges to within this value

a,b: final intercept and slope a.err, b.err: estimated uncertainties in intercept and slope

# WAVE DEFINITIONS #

Xw = 1/(Xstd^2) #X weights
Yw = 1/(Ystd^2) #Y weights


# ITERATIVE CALCULATION OF SLOPE AND INTERCEPT #

b = b0
b.diff = tol + 1
while(b.diff>tol)
{
    b.old = b
    alpha.i = sqrt(Xw*Yw)
    Wi = (Xw*Yw)/((b^2)*Yw + Xw - 2*b*Ri*alpha.i)
    WiX = Wi*X
    WiY = Wi*Y
    sumWiX = sum(WiX, na.rm = TRUE)
    sumWiY = sum(WiY, na.rm = TRUE)
    sumWi = sum(Wi, na.rm = TRUE)
    Xbar = sumWiX/sumWi
    Ybar = sumWiY/sumWi
    Ui = X - Xbar
    Vi = Y - Ybar

    Bi = Wi*((Ui/Yw) + (b*Vi/Xw) - (b*Ui+Vi)*Ri/alpha.i)
    wTOPint = Bi*Wi*Vi
    wBOTint = Bi*Wi*Ui
    sumTOP = sum(wTOPint, na.rm=TRUE)
    sumBOT = sum(wBOTint, na.rm=TRUE)
    b = sumTOP/sumBOT

    b.diff = abs(b-b.old)
  }     

   a = Ybar - b*Xbar
   wYorkFitCoefs = c(a,b)

# ERROR CALCULATION #

Xadj = Xbar + Bi
WiXadj = Wi*Xadj
sumWiXadj = sum(WiXadj, na.rm=TRUE)
Xadjbar = sumWiXadj/sumWi
Uadj = Xadj - Xadjbar
wErrorTerm = Wi*Uadj*Uadj
errorSum = sum(wErrorTerm, na.rm=TRUE)
b.err = sqrt(1/errorSum)
a.err = sqrt((1/sumWi) + (Xadjbar^2)*(b.err^2))
wYorkFitErrors = c(a.err,b.err)

# GOODNESS OF FIT CALCULATION #
lgth = length(X)
wSint = Wi*(Y - b*X - a)^2
sumSint = sum(wSint, na.rm=TRUE)
wYorkGOF = c(sumSint/(lgth-2),sqrt(2/(lgth-2))) #GOF (should equal 1 if assumptions are valid), #standard error in GOF

# OPTIONAL OUTPUTS #

if(printCoefs==1)
 {
    print(paste("intercept = ", a, " +/- ", a.err, sep=""))
    print(paste("slope = ", b, " +/- ", b.err, sep=""))
  }
if(makeLine==1)
 {
    wYorkFitLine = a + b*X
  }
 ans=rbind(c(a,a.err),c(b, b.err)); dimnames(ans)=list(c("Int","Slope"),c("Value","Sigma"))
return(ans)
 }
$\endgroup$
  • $\begingroup$ Also note, the R package "IsoplotR" includes the york() function, giving the same results as the YorkFit code here. $\endgroup$ – Steven Wofsy Sep 23 at 13:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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