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.
