3
$\begingroup$

I have a set of points and I would like to fit a linear regression model to them, where each point has its own error value, and I want to find the gradient of the regression line. How do I calculate the standard error on the gradient? In R particularly would be helpful.

Edit: I feel like I am misusing these terms horribly, so here is what my data looks like:

enter image description here

with a lm fitted line in R. The gradient of the line corresponds to a physical quantity, I need to find the value and standard error in this quantity.

Edit 2: I feel that I should point out that in this case, all the data values have the same precision. I am interested in both this and the general case though.

$\endgroup$
3
  • 1
    $\begingroup$ What do you mean by the gradient of the regression line? Do you mean the normal of the hyperplane that fits the data the best? Also, what do you mean by the standard error on the gradient? $\endgroup$ – Sobi Dec 13 '15 at 21:49
  • 1
    $\begingroup$ What do the "error values" represent, exactly? $\endgroup$ – Glen_b Dec 15 '15 at 1:20
  • $\begingroup$ Error values are our best estimate of 1 sigma precision in the measuring apparatus. $\endgroup$ – kylergs Dec 15 '15 at 12:46
1
$\begingroup$

One approach that comes to my mind looking at your data is to conduct a simulation study. You have five mean values $\bar y_1,...,\bar y_5$ and five corresponding standard deviations $s_1,...,s_5$. It seems from the plot you provided that standard errors for the groups do not differ much from each other (so sample sizes and standard deviations are probably close for the groups), but I assumed more general case when they differ.

Knowing all this you can simulate different regression slopes by sampling $r=1,...,R$ times new values for $y_1^{(r)},..., y_5^{(r)}$ groups from normal distributions (but other choice is also possible, e.g. $t$-distribution)

$$ y_i^{(r)} \sim \mathrm{Normal}(\bar y_i, s_i) $$

and then estimating regression using those values

$$ y_i^{(r)} = \beta_0^{(r)} + \beta_1^{(r)} x_i + \varepsilon_i^{(r)} $$

using $\beta_0^{(r)}$ and $\beta_1^{(r)}$ values from each of the simulation repetitions you can compute the average slope and intercept and compute confidence intervals for those values and around the regression line applying the same methods like you would do with bootstrap results (e.g. using quantiles).

Below you can see R code for such simulation.

# generating example data

set.seed(123)

N <- 5
n <- sample(20, N, replace = TRUE)
s <- runif(N, 0, 3)
m <- runif(N, 0, 5)
X <- rnorm(sum(n), rep(m, n), rep(s, n))
x <- tapply(X, rep(1:N, n), mean)
y <- -6.2 * x + 2 + rnorm(N)

# simulation

f <- function() {
  ysamp <- rnorm(y, y, s)
  fit <- lm(ysamp ~ x)
  out <- c(coef(fit),
           fitted(fit),
           ysamp)
  names(out)[-c(1:2)] <- paste0(rep(c("yhat", "ysim"), each = N), 1:5)
  out
}

sim <- replicate(1e3, f())
coef <- rowMeans(sim[1:2, ])
quant.ci <- apply(sim[3:(2+N), ], 1, quantile, c(.025, .975))

enter image description here

$\endgroup$
4
  • $\begingroup$ Thanks, this is really helpful, how are you choosing the other 2 lines to be the 1 sigma values? Could you post the whole code for me to pick through? $\endgroup$ – kylergs Dec 15 '15 at 15:44
  • $\begingroup$ Those lines are the quant.ci values, but as I said, you can use other method for computing CI's, e.g. normal or $t$-approximation etc. $\endgroup$ – Tim Dec 15 '15 at 15:48
  • $\begingroup$ By using coef.quant.ci <-apply(sim[1:2, ], 1, quantile, c(.159, .881)); abline(coef=coef.quant.ci[1,],col='blue'); abline(coef=coef.quant.ci[2,],col='blue'); I get results that visually look wrong, am I misusing apply() or am I just misjudging it? $\endgroup$ – kylergs Dec 15 '15 at 16:17
  • $\begingroup$ @kylergs abline is used for drawing straight lines, you should rather use lines(x, coef.quant.ci[1,]) with slight adjustments for nicer look. $\endgroup$ – Tim Dec 15 '15 at 18:28
3
$\begingroup$

One way to approach this is to perform a weighted regression, with the weight of each of the 5 values inversely related to its estimated error, as determined from the size of each error bar. That gives a regression fit that takes into account the relative certainty about each of the mean values.

Useful hints for choosing weights can be found on this Cross Validated page.

In the case where all points have the same precision, as in an edit to the original question, you have a standard linear regression. What you call the "gradient" is usually called the "slope" or the "regression coefficient"; the lm function in R (and, as far as I know, all standard statistical software) provides information on the standard error of the regression coefficient.

$\endgroup$
1
$\begingroup$

If by gradient you mean the gradient of the linear function $f(x) = \beta_0 + x'\beta$, then the estimation of the gradient is $\hat\beta$ and one possible confidence interval is the combined confidence area of all the coefficients, that is, the CI is

$$ \hat\beta \pm \hat{SE}(\hat\beta) \cdot Z_{\alpha/2}, $$

where $\hat\beta$ is a vector that lm outputs under the name Estimate, $\hat{SE}(\hat\beta)$ are the corresponding standard deviations (outputted by lm as well) and $Z_{\alpha/2}$ is the $\alpha/2$ quantile of a standard normal distribution. These are asymptotic and quite frankly I hardly see the use of thinking of this as the "gradient" of the linear function, but that's how you might do it.

$\endgroup$

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.