# How can I tell Julia what are the parameters in a linear model?

**URL:** <https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692>\
**Category:** Statistics\
**Tags:** regression, fit, curve-fitting, glm, linear-regression\
**Created:** [September 23, 2022, 11:27am UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692 "2022-09-23T11:27:06Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 23, 2022, 11:27am UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/1 "2022-09-23T11:27:06Z")

</div>

Hi all!

So, I’m experimenting with the linear regression ecosystem (`GLM`, `Plots`, `DataFrames`, `StatsPlots` and `CSV`). I understood how to produce a simple linear regression model and a multivariate one.  
I know that OLS works even when your function is not linear in the variables, provided it is in the parameters. For instance, I tried fitting this paraboloid:

```julia
# excerpt to generate it
f((x,y)) = (x, y, x^2 + 2 * y^2 + 10*randn())

# excerpt to fit it
fm = @formula(z ~ x^2 + y^2)

```

and **it works! Yay!**

However, now I’m trying shifting the vertex:

```julia
f((x,y)) = (x, y, (x - x_0)^2 + 2 * (y - y_0)^2 + 10*randn())

```

Of course, now the estimate is off. **How can I tell Julia that there are two more parameters to fit?**

I’ve tried this

```julia
fm = @formula(z ~ (x - x_0)^2 + (y - y_0)^2)

```

but it panics when I try to fit it:

```julia
julia> fm = @formula(z ~ (x - x_0)^2 + (y - y_0)^2)
FormulaTerm
Response:
  z(unknown)
Predictors:
  (x,x_0)->(x - x_0) ^ 2
  (y,y_0)->(y - y_0) ^ 2

julia> @time(model = lm(fm, df))
ERROR: ArgumentError: There isn't a variable called 'x_0' in your data; the nearest names appear to be: x
Stacktrace:
 [1] ModelFrame(f::FormulaTerm{Term, Tuple{FunctionTerm{typeof(^), var"#87#91", (:x, :x_0)}, FunctionTerm{typeof(^), var"#89#93", (:y, :y_0)}}}, data::NamedTuple{(:x, :y, :z), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}; model::Type{LinearModel}, contrasts::Dict{Symbol, Any})
   @ StatsModels ~/.julia/packages/StatsModels/6uCgU/src/modelframe.jl:76
 [2] fit(::Type{LinearModel}, f::FormulaTerm{Term, Tuple{FunctionTerm{typeof(^), var"#87#91", (:x, :x_0)}, FunctionTerm{typeof(^), var"#89#93", (:y, :y_0)}}}, data::DataFrame, args::Nothing; contrasts::Dict{Symbol, Any}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ StatsModels ~/.julia/packages/StatsModels/6uCgU/src/statsmodel.jl:85
 [3] fit
   @ ~/.julia/packages/StatsModels/6uCgU/src/statsmodel.jl:78 [inlined]
 [4] #lm#5
   @ ~/.julia/packages/GLM/P0Ris/src/lm.jl:184 [inlined]
 [5] lm (repeats 2 times)
   @ ~/.julia/packages/GLM/P0Ris/src/lm.jl:184 [inlined]
 [6] top-level scope
   @ ./timing.jl:262 [inlined]
 [7] top-level scope
   @ ./REPL[287]:0

```

Is there a way to do it?

Thanks!

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 23, 2022, 2:31pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/2 "2022-09-23T14:31:45Z")

</div>

I was reasoning about my problem and it occured to me that expanding `(x - x_0)^2` I get:

```julia
(x - x_0)^2 = x^2 - 2x_0*x + x_0^2

```

Is the last term (`x_0^2`) making the problem non-linear? Is this the meaning of the error produced by Julia? If yes, is there a way to solve it (which I presume is another function, for non linear models)?

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [September 23, 2022, 3:02pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/3 "2022-09-23T15:02:35Z")

</div>

The problem is just about parsing, not the math.

Inside `@formula`, Julia treats all symbols as columns in a the input `DataFrame`. `x_0` is not a variable in your data frame, so its throwing an error.

Just make the `x - x_0` variables outside the `@formula`. Maybe call it `recentered_x`. Then run `lm` with a `@formula` acting on the `recentered_x`.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [September 23, 2022, 4:59pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/4 "2022-09-23T16:59:44Z")

</div>

Expanding is the right idea, yet, from the perspective of the model x\_0 is a constant parameter, i.e.,

```julia
(x - x_0)^2 = x^2 - 2 * x_0 * x + x_0^2 = x^2 + b * x + c

```

Thus, you could fit the more general model with the formula `z ~ x^2 + x + 1`. I’m not quite sure, how you could impose the constraint that `b = 2*x_0` and `c = x_0^2`, i.e., `c = (b/2)^2`, in a generalized linear model.

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 23, 2022, 5:31pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/5 "2022-09-23T17:31:48Z")

</div>

Thanks to both!  
Actually, I wanted to estimate _both_ the coefficients, AND the vertex coordinates (`x_0` and `y_0`). But if I did my calcs right, this makes the problem non linear, so using `GLM.lm` is not the way to go. **Is this correct? And if so, does Julia have a quadratic equivalent?**

P.S. Interestingly enough, MATLAB _can_ evaluate it (so I suppose it recognizes it and switches to a non-linear model), but OTOH it cannot deal with more than 2 predictors (for what I understood).

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [September 23, 2022, 5:34pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/6 "2022-09-23T17:34:59Z")

</div>

So you want to solve the problem below?

\min\_{x\_0, \vec{\beta}} \ \sum\_{i = 1}^{N} \left(y\_i - \beta\_1 (x\_{1i} - x\_0) - \beta\_2(x\_{2i} - x\_0)\right)^2

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 23, 2022, 5:38pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/7 "2022-09-23T17:38:53Z")

</div>

Yes, exactly

BTW sorry for the noob question: how do you typeset LaTeX in posts?

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [September 23, 2022, 5:39pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/8 "2022-09-23T17:39:55Z")

</div>

GLM.jl can’t do that for you. Hopefully someone can chime in with an appropriate tutorial for an optimization package.

You do math in Discourse

```julia
$$
Like this
$$

```

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 23, 2022, 5:45pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/9 "2022-09-23T17:45:12Z")

</div>

Thanks (for both things)! Then it’s what I suspected. I don’t know enough about neither Julia nor regression analysis to create it myself, so I’ll check for other packages or revert to MATLAB/Octave if needed. Kudos!

---

<div class="post-metadata">

**Author:** ![ayushpatnaikgit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ayushpatnaikgit/32/25229_2.png) [@ayushpatnaikgit](https://discourse.julialang.org/u/ayushpatnaikgit)\
**Post date:** [September 24, 2022, 5:14am UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/10 "2022-09-24T05:14:04Z")

</div>

Perhaps something like this:

```julia
using Optim
using DataFrames
df = DataFrame(y = rand(10), x₁ = rand(10), x₂ = rand(10), x₀ = repeat([rand(1:10)], 10))

f(β) = sum((df.y .- β[1]*(df.x₁ .- df.x₀) - β[2]*(df.x₂ .- df.x₀)).^2) # equation written by pdeffebach

result = optimize(f, [0.0, 0.0]) # Ideally, put suitable initial values instead of [0.0,0.0]
Optim.minimizer(result)

```

If x₀ is not a constant, it can be β[3].

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 26, 2022, 10:36am UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/11 "2022-09-26T10:36:24Z")

</div>

Unfortunately, this doesn’t work. Even using the correct initial values, it overestimates a parameter. Let me explain:  
I have this code to generate a DataFrame (`generate_quadratic.jl`):

```julia
using DataFrames
using CSV

x = range(-50, 50, length=101);
y = range(-70, 50, length=51);
domain = Iterators.product(x,y) |> collect
x_0 = 10;
y_0 = 20;
f((x,y)) = (x, y, (x - x_0)^2 + 2 * (y - y_0)^2 + 10*randn())
points = vec(map(f, domain))

@time(df = DataFrame(NamedTuple{(:x, :y, :z)}.(points)))
# print(df)

CSV.write("raw_quadratic_data.csv", df);

```

Then I tweaked your suggestion like this to estimate it (`quadratic.jl`):

```julia
using Optim
using DataFrames
using StatsPlots

df = DataFrame(CSV.File("raw_quadratic_data.csv"));

f(β) = sum((df.z .- β[1]*(df.x .- β[3]).^2 - β[2]*(df.y .- β[4])).^2) # equation written by pdeffebach

@df df scatter(:x, :y, :z)
result = optimize(f, [1.0, 2.0, 10.0, 20.0]) # Ideally, put suitable initial values instead of [0.0,0.0]
coefficients = Optim.minimizer(result)

best_fit_surface(x, y) = coefficients[1] * (x - coefficients[3])^2 + coefficients[2] * (y - coefficients[4])^2;

x = range(-50, 50, length=101);
y = range(-70, 50, length=51);
graph = surface!(x, y, best_fit_surface)

savefig(graph, "tentative_quadratic_plot.png");

```

I got in the REPL

```julia
julia> coefficients = Optim.minimizer(result)
4-element Vector{Float64}:
    1.0003341599739208
 -120.00187363980439
    9.997530795236802
   25.795981453902133

```

and the resulting plot is:  
 ![tentative_quadratic_plot](https://global.discourse-cdn.com/julialang/original/3X/b/5/b59396f69079b5b32094f7ce76508f043a8a8cf4.png)

Am I doing something wrong?

---

<div class="post-metadata">

**Author:** ![ayushpatnaikgit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ayushpatnaikgit/32/25229_2.png) [@ayushpatnaikgit](https://discourse.julialang.org/u/ayushpatnaikgit)\
**Post date:** [September 26, 2022, 1:53pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/12 "2022-09-26T13:53:24Z")

</div>

I have made minor changes, please have a look:

```julia
using DataFrames
using Optim

x = range(-50, 50, length=101);
y = range(-70, 50, length=51);
domain = Iterators.product(x,y) |> collect
β = [1.0, 2.0, 10.0, 20.0]
f((x,y)) = (x, y, (β[1]*(x - β[3])^2 + β[2] * (y - β[4])^2 + 10*randn()))
points = vec(map(f, domain))

df = DataFrame(NamedTuple{(:x, :y, :z)}.(points))

function g(β) # since f is already used. 
    z = df.z; x = df.x; y = df.y
    x₀ = β[3]; y₀ = β[4]
    ẑ = β[1]*(x .- x₀).^2 + β[2]*(y .- y₀).^2 
    return sum((z - ẑ).^2)
end

result = optimize(g, [1.0, 2.0, 10.0, 20.0]) # Actual initial conditions
coefficients = Optim.minimizer(result) # looks like it's the right answer. 

result = optimize(g, [1.6, 20.0, 0.0, 0.0]) # Some random initial conditions. 
coefficients = Optim.minimizer(result) # Still the correct answer. 

```

---

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 26, 2022, 4:43pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/13 "2022-09-26T16:43:12Z")

</div>

**Indeed it works!** 🎉🎊 Thank you!

I’m not sure if I can flag two answers as solutions. In case I can’t, I’m going to credit you in a separate thread, to keep the question-solution pair more coherent.

---

<div class="post-metadata">

**Author:** ![ayushpatnaikgit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ayushpatnaikgit/32/25229_2.png) [@ayushpatnaikgit](https://discourse.julialang.org/u/ayushpatnaikgit)\
**Post date:** [September 26, 2022, 4:51pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/14 "2022-09-26T16:51:51Z")

</div>

Let me know if you need help with making it faster.

There are several algorithms for optimisation, LBFGS if often the fastest. You can also use automatic differentiation. I think adding constraints on beta will also make it faster.

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [September 26, 2022, 5:24pm UTC](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/15 "2022-09-26T17:24:11Z")

</div>

There are many performance improvements to be had before changing the algorithm. In particular, you want to “capture” the variables used in the optimization inside a `let` block.

Also, never do `df.z` in performant code, since `DataFrame`s are type-unstable. Always work on the vector directly.

```julia
julia> z = df.z; x = df.x; y = df.y;

julia> g = let x = x, y = y, z = z 
           β -> begin 
                   x₀ = β[3]; y₀ = β[4]
                   ẑ = β[1]*(x .- x₀).^2 + β[2]*(y .- y₀).^2 
                   return sum((z - ẑ).^2)
           end
       end;

```