# Plot the confidence interval for a model fit

**URL:** <https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767>\
**Category:** Statistics\
**Tags:** plotting, fit, glm\
**Created:** [April 17, 2020, 3:52pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767 "2020-04-17T15:52:05Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [April 17, 2020, 3:52pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/1 "2020-04-17T15:52:05Z")

</div>

I’m no data scientist, so it’s likely that my question does not make much sense… But here goes:

I recently tried to do some data analysis that involved fitting models. I’m able to fit the model, compute represent a prediction that it produces, but I’d like to also (compute and) represent confidence intervals for the model.

Here is my minimal (but incomplete) example:

```julia
using DataFrames
data = DataFrame(x = rand(100));
data.y = 1 .+ 2*data.x .+ 0.1*rand(100);

using GLM
model = lm(@formula(y ~ x), data)
pred = DataFrame(x = 0:0.01:1);
pred.y = predict(model, pred);

using Plots
plot(xlabel="x", ylabel="y", legend=:bottomright)
plot!(data.x, data.y, label="data", seriestype=:scatter)
plot!(pred.x, pred.y, label="model", linewidth=3)
savefig("/tmp/plot.png")

```

![plot](https://global.discourse-cdn.com/julialang/original/3X/8/4/84b7b7376ce81d5240c2b9bd85d1af15cb62c32e.png)

To rephrase (because I’m not even sure to use the correct words here): I’d like to represent the uncertainty about the model using a shaded area around the “model” curve.

Is there a (not too complicated) way to do that? Ideally, the solution should be as independent of the model as possible, because my real use-case involves more complex models…

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 17, 2020, 8:20pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/2 "2020-04-17T20:20:12Z")

</div>

If you replace your last `plot!` statement (before `savefig`) with:

```julia
plot!(pred.x, pred.y, ribbon=(fill(0.3,length(pred.x)),0.5.+0.1*cos.(10pi*pred.x)),fc=:orange,fa=0.3,label="model", linewidth=3)

```

you get:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/8/6885d702940f93173c3c8d367a9660004a7aa224.png)  
OK – my example of `ribbon` argument is not realistic, but then I don’t know how to use GLM… If GLM lets you produce prediction confidence intervals, my example should give an idea of how you can do it.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 9:06pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/3 "2020-04-17T21:06:36Z")

</div>

Could you try `]add StatsPlots#mkb/statsmodels`? It’s the branch of this PR [https://github.com/JuliaPlots/StatsPlots.jl/pull/293](https://github.com/JuliaPlots/StatsPlots.jl/pull/293)  
Check the example in the first comment for usage. If you comment on the PR with your experiences we can merge it and make this available to everyone.

---

<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:** [April 17, 2020, 9:06pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/4 "2020-04-17T21:06:43Z")

</div>

EDIT: I accidentally led this thread down a rabbit hole by plotting a _prediction_ interval rather than a _confidence_ interval below, and also bungling the use of `ribbon`. Long story short the ggplot2 plot shown below can also be obtained in Plots.jl/GLM.jl if one calls `predict(model, pred, interval = :confidence, level = 0.95)`.

@BLI has the right plot command, let me add the GLM command you are probably after:

```julia

using DataFrames, GLM, Plots
data = DataFrame(x = rand(100));
data.y = 1 .+ 2*data.x .+ 0.1*rand(100);

model = lm(@formula(y ~ x), data)
pred = DataFrame(x = 0:0.01:1);
pr = predict(model, pred, interval = :prediction, level = 0.95)

plot(xlabel="x", ylabel="y", legend=:bottomright)
plot!(data.x, data.y, label="data", seriestype=:scatter)
plot!(pred.x, pr.prediction, label="model", linewidth=3, 
        ribbon = (pr.prediction .- pr.lower, pr.upper .- pr.prediction))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/c/ccd5e3396801b969a79277239b618fc7685b7371.png)

The relevant docstring is:

```julia
?help> predict

  predict(mm::LinearModel, newx::AbstractMatrix;
          interval::Union{Symbol,Nothing} = nothing, level::Real = 0.95)

  If interval is nothing (the default), return a vector with the predicted values for model mm and new data newx. Otherwise, return a 3-column matrix with the prediction and the lower and upper confidence bounds for a given level (0.95 equates alpha = 0.05).
  Valid values of interval are :confidence delimiting the uncertainty of the predicted relationship, and :prediction delimiting estimated bounds for new data points.

```

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 9:07pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/5 "2020-04-17T21:07:48Z")

</div>

@nilshg I believe the ribbon is called wrongly in your example

---

<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:** [April 17, 2020, 9:11pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/6 "2020-04-17T21:11:43Z")

</div>

Ah yes sorry as usual I got confused with what `ribbon` does - let me fix that

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 9:17pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/7 "2020-04-17T21:17:38Z")

</div>

They are in fact so small you almost can’t see the ribbon. This is with your code using the StatsPlots PR:

```julia
using DataFrames, StatsPlots, GLM
data = DataFrame(x = rand(100));
data.y = 1 .+ 2*data.x .+ 0.1*rand(100);
model = lm(@formula(y ~ x), data)

@df data scatter(:x, :y)
plot!(model)

```

 ![Screenshot 2020-04-17 at 23.17.21](https://global.discourse-cdn.com/julialang/original/3X/c/c/cc2428c98ba4743a4ebb8060bd8a8c4fc192830b.png)

---

<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:** [April 17, 2020, 9:21pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/8 "2020-04-17T21:21:16Z")

</div>

Now I’m confused - for my first prediction I get [0.98, 1.10] for upper and lower confidence interval, which looks about right for the ribbon my plot above. What is StatsPlots plotting?

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 9:28pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/9 "2020-04-17T21:28:15Z")

</div>

This is how it looks in ggplot2 (same data with RCall)

 ![Screenshot 2020-04-17 at 23.27.06](https://global.discourse-cdn.com/julialang/original/3X/a/7/a7d293a2607661c9e707eb4a17de9180d58c6178.png)  
I think you swapped the upper and the lower confidence interval perhaps? The `ribbon` keyword really sucks for this kind of thing in fact.

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [April 17, 2020, 10:01pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/10 "2020-04-17T22:01:37Z")

</div>

Many thanks to both of you!

I have been able to successfully use both @nilshg’s “manual” method and @mkborregaard’s `StatsPlot` branch to get the same plots in my simple example.

In my real case, I’m only able to use the “manual” method, which I guess is because my model is slightly mode complicated: `@formula(y ~ x + x^2)`. It looks like I’m hitting the limitations mentioned in the [PR](https://github.com/JuliaPlots/StatsPlots.jl/pull/293) : “it does only address the bivariate case”.

  

As a follow-up, would there be a way to extend this to Generalized Linear Models? It looks like `predict` does not support the `interval` keyword for such models…

For example, the previous case, with

```julia
model = glm(@formula(y ~ x), data, Normal(), IdentityLink())

```

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 10:13pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/11 "2020-04-17T22:13:55Z")

</div>

You’re right, it’s not yet implemented on GLM (it isn’t in R either if that’s any help). There’s a PR on GLM.jl for it but its bugged (IMHO).  
Could you comment the polynomial case on the PR in StatsPlots? That seems like the first thing to add.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [April 17, 2020, 10:15pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/12 "2020-04-17T22:15:58Z")

</div>

Hilariously I contributed the `predict` with `intervals` to GLM like 3 years ago, exactly because I wanted to create support for this - and I still haven’t really come around to it. So nice with a bit of a push.

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [April 17, 2020, 10:26pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/13 "2020-04-17T22:26:46Z")

</div>

> [@mkborregaard](#):
>
> Could you comment the polynomial case on the PR in StatsPlots? That seems like the first thing to add.

Will do. Thanks!

---

<div class="post-metadata">

**Author:** ![hafez\_ahmad](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hafez_ahmad/32/14423_2.png) [@hafez\_ahmad](https://discourse.julialang.org/u/hafez_ahmad)\
**Post date:** [August 15, 2020, 9:32am UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/14 "2020-08-15T09:32:53Z")

</div>

I am getting error (julia 1.5)  
plot!(model)

ERROR: Cannot convert StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}} to series data for plotting  
Stacktrace:  
[1] error(::String) at .\error.jl:33  
[2] \_prepare\_series\_data(::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}) at C:\Users\hafez.julia\packages\RecipesPipeline\5RD7m\src\series.jl:8  
[3] \_series\_data\_vector(::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}, ::Dict{Symbol,Any}) at C:\Users\hafez.julia\packages\RecipesPipeline\5RD7m\src\series.jl:27  
[4] macro expansion at C:\Users\hafez.julia\packages\RecipesPipeline\5RD7m\src\series.jl:139 [inlined]  
[5] apply\_recipe(::Dict{Symbol,Any}, ::Type{RecipesPipeline.SliceIt}, ::Nothing, ::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}, ::Nothing) at C:\Users\hafez.julia\packages\RecipesBase\AN696\src\RecipesBase.jl:282  
[6] \_process\_userrecipes!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}}) at C:\Users\hafez.julia\packages\RecipesPipeline\5RD7m\src\user\_recipe.jl:35  
[7] recipe\_pipeline!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}}) at C:\Users\hafez.julia\packages\RecipesPipeline\5RD7m\src\RecipesPipeline.jl:68  
[8] \_plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}}) at C:\Users\hafez.julia\packages\Plots\ViMfq\src\plot.jl:167  
[9] #plot!#127 at C:\Users\hafez.julia\packages\Plots\ViMfq\src\plot.jl:158 [inlined]  
[10] plot!(::Plots.Plot{Plots.GRBackend}, ::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}) at C:\Users\hafez.julia\packages\Plots\ViMfq\src\plot.jl:155  
[11] plot!(::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}; kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}) at C:\Users\hafez.julia\packages\Plots\ViMfq\src\plot.jl:150  
[12] plot!(::StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}) at C:\Users\hafez.julia\packages\Plots\ViMfq\src\plot.jl:144  
[13] top-level scope at REPL[12]:1  
[14] include\_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [August 16, 2020, 8:44pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/15 "2020-08-16T20:44:47Z")

</div>

did yuo check out the pr?

---

<div class="post-metadata">

**Author:** ![floswald](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floswald/32/195_2.png) [@floswald](https://discourse.julialang.org/u/floswald)\
**Post date:** [September 23, 2021, 12:59pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/16 "2021-09-23T12:59:47Z")

</div>

to revive this, I currently get

```julia
data = DataFrame(x = rand(100));
data.y = 1 .+ 2*data.x .+ rand(100);
r = lm( y ~ x, data)
pred = GLM.predict(r, data, interval = :confidence, level = 0.95)
p = @df d scatter(:x, :y, leg = false, 
xlab = "Model", 
ylab = "Data")
plot!(p, data[!,x], pred.prediction, linewidth = 2,
ribbon = (pred.prediction .- pred.lower, pred.upper .- pred.prediction))

```

 ![Screenshot 2021-09-23 at 14.59.19](https://global.discourse-cdn.com/julialang/original/3X/8/3/83ad5462904c0ce1fdceaa7576c942ec00a3aade.png)

---

<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:** [September 23, 2021, 1:40pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/17 "2021-09-23T13:40:03Z")

</div>

FWIW, say for QC, some shameless plug:

```julia
using LinearFitXYerrors, DataFrames
data = DataFrame(x = rand(100));
data.y = 1 .+ 2*data.x .+ rand(100);
st = linearfitxy(data.x, data.y, isplot=true, ratio=:auto)

```

 ![Plot_linearfitxy_confidence_interval](https://global.discourse-cdn.com/julialang/original/3X/f/8/f83177f2f13b12c98d144e8d90bd80407bfef0c3.png)

---

<div class="post-metadata">

**Author:** ![floswald](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floswald/32/195_2.png) [@floswald](https://discourse.julialang.org/u/floswald)\
**Post date:** [September 23, 2021, 2:12pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/18 "2021-09-23T14:12:04Z")

</div>

Nice! Is this a standalone plotting package?

---

<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:** [September 23, 2021, 2:32pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/19 "2021-09-23T14:32:27Z")

</div>

@floswald, [LinearFitXYerrors.jl](https://github.com/rafael-guerra-www/LinearFitXYerrors.jl) is a small Julia package to perform 1D linear fitting of experimental data with possibly correlated errors in X and Y variables, based in York et al. (2004).

`linearfitxy()` returns the regression results and uncertainties in a structure. A plot is displayed by setting optional parameter `isplot = true`

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [September 24, 2021, 12:16pm UTC](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767/20 "2021-09-24T12:16:22Z")

</div>

@floswald for that to work the data you predict on need to be sorted

[Next page](https://discourse.julialang.org/t/plot-the-confidence-interval-for-a-model-fit/37767.md?page=2)