Here is a Julia 1.x update for a response surface model:

```
import Cairo
using DataFrames, Gadfly, RDatasets
module t
using DataFrames, Statistics
VecOrSA = Union{Vector, SubArray}
expand_grid(xv::AbstractRange, yv::AbstractRange) =
vcat([[x y] for x in xv, y in yv]...)
function rsm(xv::VecOrSA, yv::VecOrSA, zv::VecOrSA, k::Int;
xlim::Tuple=extrema(xv), ylim::Tuple=extrema(yv))
# k is the order of the polynomial surface
y = Vector{Float64}(zv)
μx, μy = mean(xv), mean(yv)
p = reshape([1:k;],1,:)
X = hcat(ones(length(xv)), (xv.-μx).^p, (yv.-μy).^p)
B = X'X \ X'y
xp = range(xlim[1], xlim[2], length=11)
yp = range(ylim[1], ylim[2], length=11)
g = expand_grid(xp, yp)
X = hcat(ones(size(g,1)), (g[:,1].-μx).^p, (g[:,2].-μy).^p)
return DataFrame(x=g[:,1], y=g[:,2], z=X*B)
end
end
Doz = dropmissing(dataset("datasets","airquality"), disallowmissing=true)
Doz[:SolarRad] = Doz[:,2] .> 200
labelD = Dict(false=>"SolarRad ≤200", true=>"SolarRad >200")
Dpredict = by(Doz, :SolarRad, [:Wind, :Temp, :Ozone] =>
x->t.rsm(x.Wind, x.Temp, x.Ozone, 2, xlim=(0, 21), ylim=(50,100)) )
set_default_plot_size(16cm, 8cm)
p = plot(Doz, x=:Wind, y=:Temp, xgroup=:SolarRad, color=:Ozone,
Geom.subplot_grid(
layer( Geom.point),
layer(Dpredict, x=:x, y=:y, z=:z, xgroup=:SolarRad, Geom.contour),
Guide.xticks(ticks=[0:5:20;]) ),
Scale.color_continuous(minvalue=0, maxvalue=200),
Scale.xgroup(labels=x->labelD[x], levels=[false, true])
)
```