# Linear regression with a positive definite matrix in GLM.jl?

**URL:** <https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638>\
**Category:** Statistics\
**Tags:** glm\
**Created:** [September 6, 2018, 7:26pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638 "2018-09-06T19:26:26Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![sebastianskejoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebastianskejoe/32/5348_2.png) [@sebastianskejoe](https://discourse.julialang.org/u/sebastianskejoe)\
**Post date:** [September 6, 2018, 7:26pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/1 "2018-09-06T19:26:26Z")

</div>

Hi

I’m trying to perform linear regression on a small dataset (these are just sample data, the real dataset will be much larger), but I’m running into a problem with the matrix not being positive definite, which it very well might be, but I wouldn’t know since I have no idea what that means 🙂 However, if I perform linear regression in R, it does not complain. So how can I work around this in Julia?

Julia:

```julia
julia> using GLM

julia> using CSV

julia> d = CSV.read("test.csv")
10×6 DataFrame
│ Row │ x1 │ x2 │ x3 │ x4 │ x5 │ y │
├─────┼───────┼─────────┼─────────┼────┼───────┼─────────┤
│ 1 │ 187.7 │ 84.1768 │ 82.4464 │ 10 │ 76.81 │ 2714.95 │
│ 2 │ 187.7 │ 87.5767 │ 84.3167 │ 12 │ 76.81 │ 2754.11 │
│ 3 │ 187.7 │ 91.9267 │ 86.0667 │ 14 │ 76.81 │ 4480.88 │
│ 4 │ 187.7 │ 91.5833 │ 89.1167 │ 16 │ 76.81 │ 4760.54 │
│ 5 │ 187.7 │ 86.915 │ 92.35 │ 18 │ 76.81 │ 5205.43 │
│ 6 │ 178.4 │ 106.757 │ 81.1379 │ 10 │ 78.97 │ 2598.67 │
│ 7 │ 178.4 │ 108.783 │ 82.9833 │ 12 │ 78.97 │ 4295.96 │
│ 8 │ 178.4 │ 103.832 │ 86.5667 │ 14 │ 78.97 │ 4647.4 │
│ 9 │ 178.4 │ 96.5583 │ 90.2203 │ 16 │ 78.97 │ 4300.89 │
│ 10 │ 178.4 │ 93.7767 │ 92.9333 │ 18 │ 78.97 │ 4254.61 │

julia> lm(@formula(y ~ x1 + x2 + x3 + x4 + x5), d)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
...

```

R:

```julia
> d = read.csv("test.csv")
> lm(y ~ x1 + x2 + x3 + x4 + x5, d)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = d)

Coefficients:
(Intercept) x1 x2 x3 x4 x5  
   -9234.73 91.71 64.45 -206.79 590.43 NA 

```

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [September 6, 2018, 8:22pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/2 "2018-09-06T20:22:51Z")

</div>

Someone else can give you a more technical answer on the meaning of positive definite.

My guess is that the error is because your `x1` and `x5` terms vary together, so you don’t want to use them as separate terms. Note the `NA` in the R results for `x5`. Try taking one of those terms out.

---

<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 6, 2018, 9:22pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/3 "2018-09-06T21:22:31Z")

</div>

R and Stata silently omit duplicate variables from regression output, and I like this feature a lot. It’s useful for heterogeneity analysis with something in your existing vector of covariates.

It would be nice if Julia had this behavior, too.

---

<div class="post-metadata">

**Author:** ![alejandromerchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alejandromerchan/32/10500_2.png) [@alejandromerchan](https://discourse.julialang.org/u/alejandromerchan)\
**Post date:** [September 7, 2018, 12:05am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/4 "2018-09-07T00:05:38Z")

</div>

As @tshort pointed out, this is due to collinear variables X1 and X5. Is that a feature of your dataset or something you inserted to test the package? There are ways to deal with collinearity problems, such as PCA, but if you’re just testing, try removing X5, or change its values and you should obtain an answer.

---

<div class="post-metadata">

**Author:** ![tk3369](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tk3369/32/2824_2.png) [@tk3369](https://discourse.julialang.org/u/tk3369)\
**Post date:** [September 7, 2018, 6:29am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/5 "2018-09-07T06:29:20Z")

</div>

It may be good to figure out how R does it and consider making an enhancement in GLM?

---

<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:** [September 7, 2018, 7:13am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/6 "2018-09-07T07:13:48Z")

</div>

FWIW, I don’t think that automagically doing something clever (automatic removal of dependent columns, regularization, shrinkage priors) is the right solution, as it just hides the problem from the user.

OTOH, an informative error message would go a long way.

---

<div class="post-metadata">

**Author:** ![sebastianskejoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebastianskejoe/32/5348_2.png) [@sebastianskejoe](https://discourse.julialang.org/u/sebastianskejoe)\
**Post date:** [September 7, 2018, 7:35am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/8 "2018-09-07T07:35:49Z")

</div>

Thanks for the quick replies! The two columns in question are weight and height, so this might sort itself out, when I get more data, but it is probably better to just combine these into body mass index. Removing either column or using BMI instead, makes it run perfectly.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [September 7, 2018, 7:51am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/9 "2018-09-07T07:51:22Z")

</div>

At least we could improve the error message. In case of errors like that, it would be nice to identify the problematic columns and mention then in the message. Can you file an issue?

Being able to get the same behavior as R and Stata could also be useful, maybe as an option. Though an explicit error can be more user-friendly than setting coefficients to `NA` without any explanation.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [September 7, 2018, 8:18am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/10 "2018-09-07T08:18:02Z")

</div>

We have this functionality already but, unfortunately, it is currently undocumented. It also seems to remove the intercept which might be a little annoying.

```julia
julia> lm(@formula(y ~ x1 + x2 + x3 + x4 + x5), df, true)
StatsModels.DataFrameRegressionModel{LinearModel{LmResp{Array{Float64,1}},DensePredChol{Float64,LinearAlgebra.CholeskyPivoted{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + x1 + x2 + x3 + x4 + x5

Coefficients:
             Estimate Std.Error t value Pr(>|t|)
(Intercept) 0.0 NaN NaN NaN
x1 73.8966 47.6088 1.55216 0.1956
x2 64.4484 164.694 0.391322 0.7155
x3 -206.787 1058.1 -0.195432 0.8546
x4 590.431 1353.83 0.436118 0.6853
x5 -76.7016 1076.81 -0.0712303 0.9466

```

---

<div class="post-metadata">

**Author:** ![sebastianskejoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebastianskejoe/32/5348_2.png) [@sebastianskejoe](https://discourse.julialang.org/u/sebastianskejoe)\
**Post date:** [September 7, 2018, 11:37am UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/11 "2018-09-07T11:37:10Z")

</div>

I have opened an issue suggesting to add what @nalimilan and @andreasnoack wrote.

---

<div class="post-metadata">

**Author:** ![sagearbor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sagearbor/32/6871_2.png) [@sagearbor](https://discourse.julialang.org/u/sagearbor)\
**Post date:** [January 31, 2019, 6:56pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/12 "2019-01-31T18:56:23Z")

</div>

I wouldnt combine them into BMI. In that particular example, BMI is not the same as weight and height. For a 5 foot and 6 foot person could have the same BMI but if the height alone was responsible for an effect you would have lost that. Fat and muscular people can also have the same BMI.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [February 1, 2019, 7:37pm UTC](https://discourse.julialang.org/t/linear-regression-with-a-positive-definite-matrix-in-glm-jl/14638/13 "2019-02-01T19:37:47Z")

</div>

What functionality?