# Help understanding GLM (the concept, not the package)

**URL:** <https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896>\
**Category:** Statistics\
**Tags:** glm\
**Created:** [March 12, 2020, 6:18pm UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896 "2020-03-12T18:18:08Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [March 12, 2020, 6:18pm UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896/1 "2020-03-12T18:18:08Z")

</div>

I’m trying to understand how generalized linear models work by coding a simple example from scratch and I’m hoping some of the stats experts here (I know there are several of you lurking around 🧐) can help. In the case of simple linear regression, the code is super simple:

```julia
using Random
using StatsPlots

X = 1:10
Y = rand(collect(-1:0.1:1), 10) .+ collect(1:2:20)
β = hcat(ones(X), X) \ Y # estimate coefficients
Yᵢ(xᵢ) = β[1] + β[2]xᵢ # function to estimate Yᵢ
scatter(X, Y, legend=false, xlabel="X", ylabel="Y")
plot!(X, Yᵢ.(X))

```

![ols](https://global.discourse-cdn.com/julialang/original/3X/b/3/b33bf0305a1ff02fb0b7b76acc97aab0c7e65d93.png)

Now, I’m trying to understand how to code a simple GLM for something like this:

```julia
using DataFrames
using GLM
using Random
using StatsPlots

X = 1:10
Y2 = Y.^4 .+ Y
scatter(X, Y2, legend=false, xlabel="X", ylabel="Y")

```

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

With the GLM package, it’s obviously very easy to fit a model to this:

```julia
model = glm(@formula(Y2 ~ X), DataFrame(X=X, Y2=Y2), Normal(), LogLink())

```

But I’m trying to understand how to code something like I’ve done above, but for this example.

Wikipedia has a decent example of how to compute the expected value of Yᵢ:

![glm](https://wikimedia.org/api/rest_v1/media/math/render/svg/f28cf2b294ec8dfda8fbccb7fabc69e3fee2639a)  
I guess my specific questions are:

1. Is `g` simply the `log()` function in this case?
2. `β = hcat(ones(X), X) \ Y` no longer holds true for this case and I need to estimate `β` via MLE. Is this correct? Since Distributions.jl uses MLE for fitting distributions, is it possible to do what I need with that package or do I need to use Optim.jl to solve this problem (or am I totally wrong about all of this 😆)?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 13, 2020, 9:01am UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896/2 "2020-03-13T09:01:38Z")

</div>

1. yes,
2. `glm` should already do the fitting for you, see

[https://juliastats.org/GLM.jl/dev/manual/#Fitting-GLM-models-1](https://juliastats.org/GLM.jl/dev/manual/#Fitting-GLM-models-1)

and the examples in the manual.

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [March 13, 2020, 1:22pm UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896/3 "2020-03-13T13:22:00Z")

</div>

> [@Tamas\_Papp](#):
>
> - yes,
> - `glm` should already do the fitting for you, see

Hi Tamas - yes, GLM makes it very easy to fit models to these examples. What I’m trying to do is to figure out how to code it myself. For example, the estimates for the coefficients that are produced by

```julia
ols = lm(@formula(Y ~ X), DataFrame(X=X,Y=Y))

```

can be found very simply with`β = hcat(ones(X), X) \ Y` and then the simple function `Yᵢ(xᵢ) = β[1] + β[2]xᵢ` will produce a point estimate.

I’m trying to figure out how to do the same, but for the exponential example. In GLM, it’s

```julia
model = glm(@formula(Y2 ~ X), DataFrame(X=X, Y2=Y2), Normal(), LogLink())

```

so I need to figure out how to estimate `β` and then I believe I can plug that, along with `log` in place of `g`, into the below function to produce point estimates:

![glm](https://wikimedia.org/api/rest_v1/media/math/render/svg/f28cf2b294ec8dfda8fbccb7fabc69e3fee2639a)

---

<div class="post-metadata">

**Author:** ![dilumaluthge](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dilumaluthge/32/29283_2.png) [@dilumaluthge](https://discourse.julialang.org/u/dilumaluthge)\
**Post date:** [March 13, 2020, 1:25pm UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896/4 "2020-03-13T13:25:12Z")

</div>

If you are interested in doing maximum likelihood estimation of your coefficient vector beta, then you could try Googling “generalized linear models Fisher’s scoring” and “generalized linear models iteratively weighted least squares”.

Unlike linear models, generalized linear models do not have closed form solutions. You need to use iterative algorithms.

Here is an example of some lecture slides that might be useful: [http://civil.colorado.edu/~balajir/CVEN6833/lectures/glm-estimation-presentation.pdf](http://civil.colorado.edu/~balajir/CVEN6833/lectures/glm-estimation-presentation.pdf)

Essentially, the algorithms that you have are:

1. Newton’s method (also known as the Newton-Raphson method)
2. Fisher’s scoring method
3. Iteratively reweighted least squares (IRLS)

IRLS is basically just a rewritten form of Fisher’s scoring, so the two are exactly equivalent,

Newton-Raphson is equivalent to Fisher’s scoring if you are using the canonical link function. If your link function is not the canonical link function, then Newton-Raphson is not necessarily equivalent to Fisher’s scoring.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 13, 2020, 3:56pm UTC](https://discourse.julialang.org/t/help-understanding-glm-the-concept-not-the-package/35896/5 "2020-03-13T15:56:23Z")

</div>

As pointed out by @dilumaluthge, there is no closed form solution in general; various iterative methods exist using the special structure of the problem. Practical implementations use some standard tricks to help with the usual numerical problems (mostly conditioning), the key algorithm is [IRLS](https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares). GLM.jl [uses this](https://github.com/JuliaStats/GLM.jl/blob/master/src/glmfit.jl) and it is a nice codebase to follow along.

If you are interested in the numerical implementation of these algorithms, another nice package to read is

> **[GitHub - JuliaStats/MixedModels.jl: A Julia package for fitting (statistical)...](https://github.com/JuliaStats/MixedModels.jl)**
>
> A Julia package for fitting (statistical) mixed-effects models - GitHub - JuliaStats/MixedModels.jl: A Julia package for fitting (statistical) mixed-effects models