# How to find the best fit line in this case?

Suppose that I have some data like this: There are $$n$$ data points $$(x_i,y_i)$$ and associated with each point are standard errors, $$\sigma_{xi}$$ and $$\sigma_{yi}$$ each with confidence level of $$\sim 68\%$$.

Now I want to find the line of best fit that approximates the model best, assuming that the relation between $$x_i$$’s and $$y_i$$’s is linear. I also want to find the best estimate of the errors (with confidence level) in the parameters (slope and intercept) of the best fit line and also the best estimate for the error (again with confidence level) in a $$y$$ value predicted by the best fit line at a given $$x$$.

Of all the sources that I’ve gone through so far (Squires, Bevington and Robinson), none mentions what to do in this case. At most they only mention least squares method with wheighting. But this is not what I’m after, as this approach doesn’t take into account the absolute values of the associated errors — only their ratios to weight each $$(x_i, y_i)$$.

And clearly that is not adequate for if I scale up all the errors by a constant, this approach, though will give me the same line, but will also give me the same values for error in the parameters and $$y$$ values, which is obviously wrong.

So please guide me what to do in such crisis, or point me to a resource which does answer this.

• Hi, where do do you know the standard errors from? If you already know them, shouldn't you also know the mean then? Oct 14, 2019 at 13:11
• @jottbe Standard errors after taking each $(x_i, y_i)$ observation large number of times. Yes I know mean of each of them, but that's not what I ask. I want the best fit parameters with confidence levels.
– Atom
Oct 14, 2019 at 13:14
• Weighted reduced major axis regression? eprints.usq.edu.au/24210 The major axis regression deals with the error in both x and y, the weighting with that you actually know that error. Never used it, only a possible hint. Oct 14, 2019 at 13:26

I wanted to suggest something similar to the approach pointed by @Carsten in the comment (I also never used it, nor heard about it ealier). He refers to the paper by Saqr and Khan (2012) who describe weighted reduced major axis regression. The method assumes that there is a regression model defined in terms of latent variables $$\eta_i$$ and $$\xi_i$$

\begin{align} \eta_i &= \alpha + \beta \,\xi_i + \varepsilon_i \\ \varepsilon &\sim \mathcal{N}(0, \sigma) \end{align}

where what we observe are the draws from the random variables $$\eta_i$$ and $$\xi_i$$ that are measured with noise:

$$x_i \sim \mathcal{N}(\xi_i, \sigma_{x_i}) \\ y_i \sim \mathcal{N}(\eta_i, \sigma_{y_i})$$

where the measurement error levels $$\sigma_{x_i}$$ and $$\sigma_{y_i}$$ are known.

I didn't go into details of the referred paper, but such model can be easily estimated using probabilistic programming framework like Stan. The model definition in Stan would be something like the following code (I'm using flat priors for simplicity, but you could use more reasonable ones if possible).

data {
int<lower=0> n;
vector[n] x;
vector[n] x_sd;
vector[n] y;
vector[n] y_sd;
}
parameters {
vector[n] x_mean;
vector[n] y_mean;
real alpha;
real beta;
real<lower=0> sigma;
}
model {
x ~ normal(x_mean, x_sd);
y ~ normal(y_mean, y_sd);
y_mean ~ normal(alpha + beta * x_mean, sigma);
}


where it can be easily extended to more complicated cases.

You can verify the solution by simulating some data from the model and comparing the results from linear regression and the suggested approach.

set.seed(42)
n <- 100

alpha <- 3
beta <- 0.32
sigma <- 0.5

x_mean <- 1:10 + rnorm(n)
x_sd <- abs(rnorm(n, sd=1.5))
x <- rnorm(n, x_mean, x_sd)

y_mean <- alpha + beta * x_mean + rnorm(n, sd=sigma)
y_sd <- abs(rnorm(n, sd=2))
y <- rnorm(n, y_mean, y_sd)


The linear regression model gives us the following results:

summary(fit_lm <- lm(y ~ x))
## Call:
## lm(formula = y ~ x)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -10.7904  -0.6837  -0.0704   1.1214   6.8163
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.90409    0.52236    5.56 2.34e-07 ***
## x            0.30848    0.08075    3.82 0.000234 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 2.51 on 98 degrees of freedom
## Multiple R-squared:  0.1296, Adjusted R-squared:  0.1207
## F-statistic: 14.59 on 1 and 98 DF,  p-value: 0.0002339

sd(resid(fit_lm))
## 2.49763398072135


While the model estimated using the Stan code returns:

## Inference for Stan model: 0040187cdf195ac63987c34c5b701b58.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##       mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## alpha 3.49    0.01 0.23 3.05 3.33 3.48 3.64  3.96   901 1.00
## beta  0.25    0.00 0.04 0.17 0.22 0.25 0.27  0.32   781 1.00
## sigma 0.50    0.01 0.11 0.30 0.42 0.49 0.56  0.75   245 1.02
##
## Samples were drawn using NUTS(diag_e) at Mon Oct 14 16:20:49 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).


As you can see, the estimated parameters are close, but only the proposed model estimated the correct error standard deviation $$\sigma$$.

Below you can see the results plotted, where true and estimated regression lines are shown overlayed over the observed points (red points) and surrounding 95% confidence intervals over the true latent variables (gray intervals). Saqr, A. and Khan, S. (2012). Weighted reduced major axis method for regression model. In: 12th Islamic Countries Conference on Statistical Sciences (ICCS 2012): Statistics for Everyone and Everywhere, 19-22.