Suppose I have a set of data ${(x_i,y_i)}$ in which the uncertainty in the measurements ${(\Delta x_i,\Delta y_i)}$ (which come from the propagation of systematic errors from the measurement apparatus) is different for each point. If I do a linear regression on the set of data how do I calculate the uncertainty in the slope? I would like an explicit procedure or formula.
-
2$\begingroup$ Do you have a sense about which measurement error is bigger? $\endgroup$– dimitriyCommented Jan 18, 2018 at 1:21
-
1$\begingroup$ By the deltas, do you mean difference between sequential measurements? Is your data sequential? Do you expect the error to be correlated? Do you expect the correlation to decay? Do you have any sort of independent replications? More information is necessary to provide a concrete answer. $\endgroup$– user3903581Commented Jan 18, 2018 at 1:24
-
$\begingroup$ The term you are looking for is error propagation. You have errors on the input side and compute (more specifically: estimate) two parameters from it. Unfortunately the name "linear regression" only describes a (popular) model, but not the method by which you estimate the parameters. For the most used methods you can probably look up the solution (e.g. least squares). If not, you can either calculate it analytically or approximate it by numeric evaluation. $\endgroup$– cherubCommented Jan 22, 2018 at 15:10
4 Answers
We can model the experiment as $$x_i=x_i^*+\tilde u_i$$ $$y_i=y_i^*+\tilde v_i$$ $$\tilde u_i=\bar u + v_i$$ $$\tilde v_i=\bar v + u_i$$ where $x_i^*, y_i^*$ denote true values, $\tilde u_i,\tilde v_i $ are measurement errors, $\bar u,\bar v $ are their "fixed" components independent from observation (which could arise from wrong calibration of the sensors) and $u,v$ vary from observation to observation and correspond to many possible factors which we treat as random.
Simple linear regression is $$y_i^*=\alpha+\beta x_i^*+e_i$$ and OLS estimate of the slope is $$\hat\beta=\frac{Cov(x^*,y^*)}{Var(x^*)}$$ What we obtain is however $$\tilde\beta=\frac{Cov(x,y)}{Var(x)}=\frac{Cov(x^* + u,y^*+ v)}{Var(x^* + u)}=\frac{Cov(x^*,y^*)+Cov(x^*,v)+Cov(y^*,u)+Cov(u,v)}{Var(x^*) + Var(u) + 2Cov(x,u)}$$
Now let's assume that $v,u$ are uncorrelated with $x^*,y^*$ and each other (a rather strong assumption that can be improved if we have more inferences about the nature of errors). Then our estimate is $$\tilde\beta=\beta\frac{\sigma^2_{x^*}}{\sigma^2_{x^*}+\sigma^2_{u}}\approx\beta\frac{\hat\sigma^2_x-\hat\sigma^2_u}{\hat\sigma^2_x}=\beta\hat\lambda$$ We can estimate $\hat\sigma^2_x$ as sample variation of $x_i$. We also need to estimate $\sigma^2_u$. If we have an experiment when we can observe $x^*_i$ multiple times, then one simple approach is to estimate $\sigma^2_u=E[\sigma^2_x|x^*_i$].
Now we can use our $\hat\sigma^2_{\tilde\beta}$ calculated with, for example, bootstrap method, and correct it for $\hat\beta =\tilde\beta /\hat\lambda$ so that $$\hat\sigma^2_{\hat\beta}=\frac{\hat\sigma^2_{\tilde\beta}}{\hat\lambda^2}$$.
I think the answer given by @yshilov is definitely awesome by considering the measurement error into the error term and significantly, deduces the result $$\tilde \beta = \beta \frac{\sigma_x^2}{\sigma_x^2 + \sigma_u^2}$$
To elaborate, this beta has special properties that it is an biased estimator, but biased towards 0. Specifically, for linear regression, $E(\hat \beta_1)=\beta_1 \cdot\Big[\frac{\sigma_x^2+\sigma_{x\delta}}{\sigma_x^2+2\sigma_{x\delta}+\sigma_{\delta}^2}\Big]$
The proof is as follows: in simple linear regression, recall $$\hat \beta_1 = \frac{\sum_{i=1}^n(x_i-\bar x)y_i}{\sum_{i=1}^n(x_i-\bar x)^2}$$ In the case of measurement error, we have $x_i^O=x_i^A=\delta_i$, $y_i^O=y_i^A+\epsilon_i$, and $y_i^A=\beta_0 +\beta_1 x_i^A$, so we get $$y_i^O=\beta_0+\beta_1(x_i^O-\delta_i)+\epsilon_i=\beta_0+\beta_1x_i^O+(\epsilon_i-\beta_1 \delta_i)$$ Assuming that $E(\epsilon_i)=E(\delta_i)=0$, $var(\epsilon_i)=\sigma_{\epsilon}^2$, $var(\delta_i)=\sigma_{\delta}^2 = \frac{1}{n}\sum_{i=1}^n(\delta_i-\bar \delta)^2$ and the variance of true predictor value $\sigma_{x}^2=\frac{\sum(x_i^A-\bar {x^A})^2}{n}$ and correlation of true predictor and error $\sigma_{x \delta}=cov(x^A,\delta)= \frac{1}{n}\sum_{i=1}^n(x_i^A-\bar {x_i^A})(\delta_i- \bar \delta)$, then
$$cov(x_i^O,\delta)=E(x_i^O\delta)-E(x_i^O)\cdot E(\delta)=E(x_i^O\delta)=E[(x_i^A+\delta)\delta]=E(x_i^A \delta)+E(\delta^2)$$ $$=\big[E(x_i^A \delta)-E(x_i^A)\cdot E(\delta)\big]+\big[var(\delta)+[E(\delta)]^2\big]=cov(x_i^A,\delta)+\sigma_{\delta}^2=\sigma_{x\delta}+\sigma_{\delta}^2$$ Then, by $\bar x = E(x_i)$ and bilinearity property in covariance, the expectation of $\hat \beta_1$ is $$E(\hat \beta_1)=E\Big[\frac{\sum_{i=1}^n(x_i^O-\bar x^O)y_i^O}{\sum_{i=1}^n(x_i^O-\bar x^O)^2}\Big]=\frac{E(\sum_{i=1}^nx^O_iy_i^O)-E(\sum_{i=1}^n \bar x^Oy_i^O)}{\sum_{i=1}^n E\big[(x_i^O-E(x_i^O))^2\big]}=\frac{E(\sum_{i=1}^nx_i^Oy_i^O)-E(x_i^O)\cdot E(\sum_{i=1}^n y_i^O)}{\sum_{i=1}^nvar(x_i^O)}$$ $$=\frac{\sum_{i=1}^ncov(y_i^O,x_i^O)}{\sum_{i=1}^nvar(x_i^O)}=\frac{\sum_{i=1}^ncov(\beta_0+\beta_1x_i^O+\epsilon_i-\beta_1\delta_i,~x_i^O)}{\sum_{i=1}^nvar(x_i^O)}=\frac{\beta_1\cdot \sum_{i=1}^nvar(x_i^O)-\beta_1\cdot \sum_{i=1}^ncov(x_i^O, \delta_i)}{\sum_{i=1}^nvar(x_i^O)}$$ $$=\beta_1 \cdot \Big[ 1-\frac{{\sum_{i=1}^ncov(x_i^O,\delta_i)}/{n}}{\sum_{i=1}^nvar(x_i^A+\delta_i)/n}\Big]=\beta_1 \cdot\Big[1-\frac{\sigma_{x\delta}+\sigma_{\delta}^2}{\sigma_x^2+2cov(x_i^A,\delta_i)+\sigma_{\delta}^2}\Big] =\beta_1 \cdot\Big[\frac{\sigma_x^2+\sigma_{x\delta}}{\sigma_x^2+2\sigma_{x\delta}+\sigma_{\delta}^2}\Big]$$ , as desired. Hence, the result $E(\hat \beta_1)=\beta_1 \cdot\Big[\frac{\sigma_x^2+\sigma_{x\delta}}{\sigma_x^2+2\sigma_{x\delta}+\sigma_{\delta}^2}\Big]$ is well-established.
I have a similar problem - posted here - and no certain answer still. What I did for the moment is simply gather a set of very similar Xs and check if there's a big variation for Y within those lines. Another kind of approach could be some a simulation: you use a single X from your dataset, but replicate the lines following the predictors systematic error (something like rnorm(...,0,0.3)). The confidence interval for slope may be something similar to the systematic error span.
I would recommend a parametric bootstrap on the data. That means generating new datasets that are similar to the real dataset, but are different to the extent implied by your uncertainty in each observation.
Here's some pseudo-code for that. Notice I'm using vector inputs to rnorm
, as is normal in the R language. Also I'm assuming that what you are calling $\Delta$ are standard errors.
For each b in 1...B:
x_PB = rnorm(x, x_se)
y_PB = rnorm(y, y_se)
r[b] = cor(x_PB, y_PB)
Then look at the distribution of the values in r.