diff --git a/README.md b/README.md
index 9c1fb73..c22b7ff 100644
--- a/README.md
+++ b/README.md
@@ -372,6 +372,18 @@ plot(
 
 ![heatmap dendrogram optimal](https://user-images.githubusercontent.com/3502975/59949464-20778680-9441-11e9-8ed7-9a639b50dfb2.png)
 
+## Plotting statistical model objects
+So far there are only recipes for the bivariate case, but the plan is to make this work for all types of statistical models
+
+```julia
+using StatsPlots, GLM, RDatasets
+iris = dataset("datasets", "iris")
+mod = lm(@formula(PetalLength ~ PetalWidth), iris)
+plot(mod)
+@df iris scatter!(:PetalWidth, :PetalLength, ms = 2, c = :black, legend = :topleft)
+```
+![model_example](https://user-images.githubusercontent.com/8429802/73453356-1ff68500-436c-11ea-872c-03ae851be281.png)
+
 ## GroupedErrors.jl for population analysis
 
 Population analysis on a table-like data structures can be done using the highly recommended [GroupedErrors](https://github.com/piever/GroupedErrors.jl) package.
diff --git a/src/StatsPlots.jl b/src/StatsPlots.jl
index de28950..95201c5 100644
--- a/src/StatsPlots.jl
+++ b/src/StatsPlots.jl
@@ -39,6 +39,7 @@ include("ecdf.jl")
 include("hist.jl")
 include("marginalhist.jl")
 include("bar.jl")
+include("statsmodels.jl")
 include("dendrogram.jl")
 include("andrews.jl")
 include("ordinations.jl")
diff --git a/src/statsmodels.jl b/src/statsmodels.jl
new file mode 100644
index 0000000..fa52f07
--- /dev/null
+++ b/src/statsmodels.jl
@@ -0,0 +1,8 @@
+# Note both predict and RegressionModel are defined in StatsBase
+@recipe function f(mod::RegressionModel)
+    newx = [ones(200) range(extrema(mod.model.pp.X[:,2])..., length = 200)]
+    newy, l, u = predict(mod, newx, interval = :confidence)
+    ribbon := (vec(u)-newy, newy-vec(l))
+    label --> "model"
+    newx[:,2], newy
+end