6
$\begingroup$

I am looking for someone to just confirm / double-check something for me with regards to errors on measurements.

Let's say I am trying to determine the slope of a relationship by varying one quantity and measuring another, and then I plot the graph and do a least-squares fit straight line to the data (graph on the left). Then I repeat this procedure twice more, to get the middle and right-most graphs. enter image description here

Each fit routune will typically give me back a slope and the corresponding 95% confidence interval, so that I obtain $(m_1\pm\Delta m_1), (m_2\pm\Delta m_2)$ and $(m_3\pm\Delta m_3)$. Now I know that the underlying quantity which determines $m$ in each case is the same, so I should be able to quote a best estimate for the slope as their mean

$$ \bar{m} = \frac{m_1+m_2+m_3}{3}. \tag{1} $$

My question is about the appropriate way to quote the error. We know that for a function $f(x,y)$ with errors in $x$ and $y$ given by $\Delta x$ and $\Delta y$, respectively, the error on $f$ is given by

$$ \Delta f = \sqrt{ (\Delta x)^2 \bigg(\frac{\partial f}{\partial x}\bigg)^2 + (\Delta y)^2 \bigg(\frac{\partial f}{\partial y}\bigg)^2 } \tag{2} $$

So I would think I can determine the error in $\bar{m}$ to be

$$ \begin{align} \Delta \bar{m} &= \sqrt{ (\Delta m_1)^2 \bigg(\frac{\partial \bar{m}}{\partial m_1}\bigg)^2 + (\Delta m_2)^2 \bigg(\frac{\partial \bar{m}}{\partial m_2}\bigg)^2 + (\Delta m_3)^2 \bigg(\frac{\partial \bar{m}}{\partial m_3}\bigg)^2} \tag{3} \\ &= \frac{1}{3} \sqrt{ (\Delta m_1)^2 + (\Delta m_2)^2 + (\Delta m_3)^2 } \tag{4} \end{align} $$

First question, is this correct?

Second question, is it okay to propagate 95% confidence intervals in this way? Should I simply quote now the result as $\bar{m} \pm \Delta \bar{m}$ and just explain that $\Delta \bar{m}$ is the combined 95% confidence interval, or should I convert the 95% number from the fits into standard errors (through the factor of 1.96)?

Thanks in advance,

(I am for now assuming Gaussian errors everywhere.)


EDIT

It was suggested in the comments that I first implement weighting in the averaging step before worrying about the errors. This should help to give more weight to slopes which have tighter confidence intervals (and vice versa).

According to this link, the weighted version of the mean would be given by $$ \bar{m}_\textrm{w} = \frac{\sum_i w_i m_i}{\sum_iw_i}, \hspace{1cm} \textrm{where} \hspace{0.5cm} w_i = \frac{1}{\sigma_i^2}\tag{5} $$ and $\sigma_i$ is the variance of each slope. Therefore, in my case with the three example slopes, it should be $$ \bar{m}_\textrm{w} = \frac{m_1/\sigma_1^2 + m_2/\sigma_2^2 + m_3/\sigma_3^2}{1/\sigma_1^2 + 1/\sigma_2^2 + 1/\sigma_3^2}. \tag{6} $$

The variance on the weighted mean slope is given at the above link again by $$ \begin{align} \textrm{Var}(\bar{m}_\textrm{w}) &= \frac{\sum_iw_i^2\sigma_i^2}{\big( \sum_iw_i\big)^2}\tag{7}\\ &= \frac{1/\sigma_1^2 + 1/\sigma_2^2 + 1/\sigma_3^2}{\big(1/\sigma_1^2 + 1/\sigma_2^2 + 1/\sigma_3^2\big)^2}\tag{8}\\ &= \big(1/\sigma_1^2 + 1/\sigma_2^2 + 1/\sigma_3^2\big)^{-1}.\tag{9} \end{align} $$

So now my main question remains - these are variances, so should we convert the 95% confidence intervals $\Delta m_i$ returned by a fitting algorthm somehow into a variance?

Maybe for a concrete example we could imagine the following values were returned from the fitting routine: $$ \begin{align} m_1 &= 5.5\; (4.9, 6.1)\rightarrow \Delta m_1 = 0.6\\ m_2 &= 5.5\; (5.3, 5.7)\rightarrow \Delta m_2 = 0.2\\ m_3 &= 5.2\; (4.5, 5.9)\rightarrow \Delta m_3 = 0.7 \end{align} $$ where the values in brackets represent the 95% confidence intervals. How should the estimate of the slope be reported, including errors? Let's imagine I only have access to these values (and not the underlying data that was used for fitting to obtain these slopes).

$\endgroup$
20
  • 1
    $\begingroup$ You can improve things using en.wikipedia.org/wiki/Inverse-variance_weighting $\endgroup$ – Jarle Tufto Jul 7 '20 at 10:23
  • 2
    $\begingroup$ To first order, you handle the confidence intervals by using them to weight your estimate, as indicated by @Jarle. Only then can you construct a reasonable confidence interval for that estimate. $\endgroup$ – whuber Jul 7 '20 at 13:34
  • 1
    $\begingroup$ Thanks both, I have edited the question to try to incorporate your advice. $\endgroup$ – teeeeee Jul 7 '20 at 14:26
  • 2
    $\begingroup$ Yes exactly, I only have access to $m_i \pm \Delta m_i$ where $\Delta m_i$ is the 95% interval. $\endgroup$ – teeeeee Jul 9 '20 at 11:56
  • 1
    $\begingroup$ Are your experiments/runs supposed to have different variance in the noise? Or can we can we consider the process to be the same in which case we can use pooled variance. $\endgroup$ – Sextus Empiricus Jul 10 '20 at 12:23
2
+50
$\begingroup$

I imagine the 95% confidence intervals come from some assumptions on normality of data. Otherwise, please state how you got these CI. This implies you believe the mean of each slope (viewed as a RV) is $m_i$ with some variance $\sigma_i$ In this case you can average the slopes as you did and get the new variance of the averaged estimator (assuming independent errors). From said variance you can get 95% CI (using 1.96 standard deviations).

So, to summarize (assuming $m_i$ are independent is crucial):

  1. Let $m := \frac{ \sum_{i=1}^N m_i \sigma_i^{-2}}{\sum_{i=1}^N \sigma_i^{-2}}$
  2. Let $\sigma^2 := Var(m) = Var(\frac{ \sum_{i=1}^N m_i \sigma_i^{-2}}{\sum_{i=1}^N \sigma_i^{-2}}) = (\frac{1}{\sum_{i=1}^N \sigma_i^{-2}})^2 \sum_{i=1}^N \sigma^{-4}Var(m_i) = \frac{1}{\sum_{i=1}^N \sigma_i^{-2}}$.
  3. Note that this is the harmonic mean! Finally you see it in the wild after learning that inequality in your first calculus class!!
  4. A 95% CI for the true value of the slope is $[m- 1.96\sigma, m+ 1.96\sigma]$ (see 1.96)
$\endgroup$
2
  • $\begingroup$ Just to check - these are exactly the same as the formulae I quoted in the question, right? $\endgroup$ – teeeeee Jul 15 '20 at 17:06
  • 1
    $\begingroup$ @teeeeee Yes they are (you wrote so much text I got lazy reading and decided it's easier to just post an answer). Note that you are relying on $m_i$'s being independent - you did not mention that anywhere and that is an assumption you are implicitly making. $\endgroup$ – Yair Daon Jul 16 '20 at 7:23
1
$\begingroup$

I don't have much to offer in terms of methods, I think the ones presented here (esp inverse variance weighted approaches) are good ones. What I can add is a small simulation study to prove that under the assumption of Gaussian errors in the regression, this process has good enough coverage


set.seed(0)

library(tidyverse)

simulate_data<-function(n){
  x = rnorm(n)
  y = 2*x + 1 + rnorm(n, 0, 0.5)
  model = lm(y~x)
  results = tibble(beta = coef(model)['x'],
                   w = 1/vcov(model)['x','x'])
}




simulate_procedure<-function(iter){
  
  n =  rnbinom(3,200,0.9)
  
  results = map_dfr(n, simulate_data)
  
  m = sum(results$beta*results$w)/sum(results$w)
  
  sig = sqrt(1/sum(results$w))
  
  interval = tibble(lower = m - 1.96*sig, 
                    est = m, 
                    upper = m + 1.96*sig)
  
  
  interval
}


map_dfr(1:10000, simulate_procedure, .id = 'iter')  %>% 
  mutate(contains = (2<upper)&(2>lower)) %>% 
  summarise(mean(contains))

>>>0.922

So what does this mean? It means that were I to repeat this procedure to construct a 95% interval for the slope, the resulting interval would capture the true slope (here 2) only 92% of the time. So barring I didn't make a mistake (entirely possible) that seems to be good enough.

How should the estimate of the slope be reported, including errors? Let's imagine I only have access to these values (and not the underlying data that was used for fitting to obtain these slopes).

So I would compute $m$ and $\sigma^2$ as mentioned by Yair Daon. You don't need to access the data in order to do these. In your example, the $m$ would be 5.5, 5.5, 5.2. The variances are found by doing a little algebra on the confidence interval. Remember, confidence intervals look like

$$m \pm 1.96 se $$

Here, $se$ is the standard error (or the standard deviation of the sampling distribution). You can find the variance by taking the difference between interval endpoints and then dividing by $3.92 = 2\times 1.96$. Your sigmas (not squared) would then be 0.306, 0.102, 0.357.

So your best estimate for $m$ from the example you've provided is 5.47, with an accompanying interval of 5.29 to 5.66. These were computed using the formulae provided by Yair.

$\endgroup$
2
  • $\begingroup$ Can you explain a little what the code is doing? I am not familiar with the programming language you have used. Thanks! $\endgroup$ – teeeeee Jul 14 '20 at 18:01
  • $\begingroup$ @teeeeee . I'm basically simulating your example 10,000 times and seeing how many times the confidence interval captures the true parameter of interest. $\endgroup$ – Demetri Pananos Jul 14 '20 at 18:32

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.