# Weighted linear regression with confidence interval fitted to error bars

**URL:** <https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743>\
**Category:** General Usage\
**Tags:** statistics, regression, fit, glm\
**Created:** [May 7, 2021, 9:10pm UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743 "2021-05-07T21:10:36Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Okarin99](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/okarin99/32/21397_2.png) [@Okarin99](https://discourse.julialang.org/u/Okarin99)\
**Post date:** [May 7, 2021, 9:10pm UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/1 "2021-05-07T21:10:36Z")

</div>

Hi,  
I’m pretty new to Data Science.  
I have some Data with measurement errors and I want to fit a linear model.

```julia
using DataFrames
using GLM
using StatsPlots

df = DataFrame(x = [0.0, 0.0669873, 0.25, 0.5, 0.75, 0.933013, 1.0],
y = [0.223, 0.291, 0.393, 0.549, 0.73, 0.85, 0.896],
u_y = [0.023, 0.024, 0.027, 0.031, 0.037, 0.041, 0.043])

# Weighted Least Squares
wlr = glm(@formula(y~x), df, Normal(), wts = 1 ./ df.u_y.^2)

x_interval = DataFrame(x = 0:0.01:1);
pred = predict(wlr, x_interval, interval = :confidence)
@df df scatter(:x, :y, yerror = :u_y)
plot!(x_interval.x, pred.prediction,
ribbon = (pred.prediction .- pred.lower, pred.upper .- pred.prediction))

```

I get this result:  
 ![plot](https://global.discourse-cdn.com/julialang/original/3X/b/d/bd58823e4efc61bf2c2092af1b965982b42e8759.png)  
The confidence interval should get plotted but it is so small that you can’t even see it.  
Shouldn’t it be as wide as the measurement errors?  
How do I fit it accordingly?

Thank your for all answers and have a nice day 🙂

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [May 8, 2021, 6:47am UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/2 "2021-05-08T06:47:33Z")

</div>

Your question is hopefully answered here [Plot the confidence interval for a model fit](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 8, 2021, 7:44am UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/3 "2021-05-08T07:44:21Z")

</div>

> [@Okarin99](#):
>
> The confidence interval should get plotted but it is so small that you can’t even see it

[GLM.jl](https://github.com/JuliaStats/GLM.jl/blob/master/src/glmfit.jl)’s `glm` does not interpret the weights `wts` as inverse variances but as prior frequencies:

```julia
    - `wts::Vector=similar(y,0)`: Prior frequency (a.k.a. case) weights of observations.
      Such weights are equivalent to repeating each observation a number of times equal
      to its weight. Do note that this interpretation gives equal point estimates but
      different standard errors from analytical (a.k.a. inverse variance) weights and
      from probability (a.k.a. sampling) weights which are the default in some other
      software.
      Can be length 0 to indicate no weighting (default).

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 8, 2021, 8:31am UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/4 "2021-05-08T08:31:17Z")

</div>

There are two different confidence intervals to consider, the confidence interval for the posterior mean, and the interval for a new measurement. You are talking about the interval for a new measurement, which contains uncertainty about both the mean and the measurement.

---

<div class="post-metadata">

**Author:** ![Okarin99](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/okarin99/32/21397_2.png) [@Okarin99](https://discourse.julialang.org/u/Okarin99)\
**Post date:** [May 8, 2021, 5:37pm UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/5 "2021-05-08T17:37:51Z")

</div>

So to get the right confidence interval I have to sum up the confidence interval I get from glm with the mean of the measurement errors?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 8, 2021, 10:12pm UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/6 "2021-05-08T22:12:52Z")

</div>

As indicated in GLM.jl’s doc above, `glm` does not handle inverse-variance weighting.

For this purpose, you may use LsqFit.jl.

```julia
using DataFrames, LsqFit, Printf, Plots; gr()

df = DataFrame(x = [0.0, 0.0669873, 0.25, 0.5, 0.75, 0.933013, 1.0],
        y = [0.223, 0.291, 0.393, 0.549, 0.73, 0.85, 0.896],
        u_y = [0.023, 0.024, 0.027, 0.031, 0.037, 0.041, 0.043])

x, y = df.x, df.y
wt = 1 ./ df.u_y .^2

p0 = [0.5, 0.5]
m(x, p) = p[1] .+ p[2] * x # p: model parameters
fit = curve_fit(m, x, y, wt, p0)

cf = coef(fit)
ci = confidence_interval(fit, 0.05) # 5% significance level

str = @sprintf("Y = (%.2f +/- %.2f) + (%.2f +/- %.2f)*X",
          cf[1],diff([ci[1]...])[1]/2, cf[2],diff([ci[2]...])[1]/2)

tl, bl = ci[1][1] .+ ci[2][2]*x, ci[1][2] .+ ci[2][1]*x
σp, σm = maximum([tl bl], dims=2) .- m(x,cf), m(x,cf) .- minimum([tl bl], dims=2)

plot(x, cf[1] .+ cf[2]*x, color=:lightblue, ribbon=(σp,σm), label=str)
plot!(x, cf[1] .+ cf[2]*x, color=:blues, lw=1, label=false, xlabel="X",ylabel="Y")
scatter!(x,y, ms=3,label=false,mc=:blue, yerror=df.u_y, legend=:topleft)

```

![LsqFit_weighted_linear_regression](https://global.discourse-cdn.com/julialang/original/3X/0/3/038fd8ca3b5ce5d8581108966b05d1d97c8cdb03.png)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2021, 11:01am UTC](https://discourse.julialang.org/t/weighted-linear-regression-with-confidence-interval-fitted-to-error-bars/60743/7 "2021-05-09T11:01:23Z")

</div>

Maybe, but probably not. It really depends on what question you want to answer. Usually, you use measurements to learn about some underlying system, and you postulate that the system behaves according to the model you have identified. You are then typically interested in learning about the posterior over the model parameters after having seen the data. The confidence interval for a new measurement is usually not very interesting since you often do not care about the measurement, you care about the thing you are trying to measure. It therefore makes sense to plot the confidence bounds like @rafael.guerra did with a shadede ribbon in the figure above, around the posterior mean of the model prediction.