# GLM - inconsistency with R on ISLR dataset?

**URL:** <https://discourse.julialang.org/t/glm-inconsistency-with-r-on-islr-dataset/75427>\
**Category:** Performance\
**Tags:** question, glm\
**Created:** [January 29, 2022, 6:38pm UTC](https://discourse.julialang.org/t/glm-inconsistency-with-r-on-islr-dataset/75427 "2022-01-29T18:38:43Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![compleat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/compleat/32/8958_2.png) [@compleat](https://discourse.julialang.org/u/compleat)\
**Post date:** [January 29, 2022, 6:38pm UTC](https://discourse.julialang.org/t/glm-inconsistency-with-r-on-islr-dataset/75427/1 "2022-01-29T18:38:43Z")

</div>

I am teaching a class using the (well-known) book Intro to Statistical Learning in R (ISLR). For the Lab example in Chapter 3, there is an example using polynomial regression in GLM.

According to documentation, the way to do this in Julia should be:

```julia
lm_fit=lm(@formula(MedV ~ LStat+LStat^2+LStat^3+LStat^4+LStat^5),Boston)

```

This give output

```julia
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}

MedV ~ 1 + LStat + :(LStat ^ 2) + :(LStat ^ 3) + :(LStat ^ 4) + :(LStat ^ 5)

Coefficients:
─────────────────────────────────────────────────────────────────────────────────────
                   Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────────────────
(Intercept) 0.0 NaN NaN NaN NaN NaN
LStat 15.8973 0.458112 34.70 <1e-99 14.9972 16.7973
LStat ^ 2 -2.60236 0.111004 -23.44 <1e-81 -2.82045 -2.38427
LStat ^ 3 0.167498 0.0091596 18.29 <1e-56 0.149502 0.185494
LStat ^ 4 -0.00472568 0.000307712 -15.36 <1e-43 -0.00533024 -0.00412112
LStat ^ 5 4.85095e-5 3.60338e-6 13.46 <1e-34 4.14299e-5 5.55891e-5
─────────────────────────────────────────────────────────────────────────────────────

```

Which seems buggy and anyway does not correspond to corresponding result run in R.

Am I understanding the usage correctly?

Thanks for any help.

---

<div class="post-metadata">

**Author:** ![jbrea](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbrea/32/3879_2.png) [@jbrea](https://discourse.julialang.org/u/jbrea)\
**Post date:** [January 29, 2022, 7:29pm UTC](https://discourse.julialang.org/t/glm-inconsistency-with-r-on-islr-dataset/75427/2 "2022-01-29T19:29:20Z")

</div>

I think you are running into an [open GLM issue](https://github.com/JuliaStats/GLM.jl/issues/426). With `res = lm(@formula(MedV ~ LStat + LStat^2 + LStat^3 + LStat^4 + LStat^5), Boston, dropcollinear=false)` you get the correct result. Note that [R’s `poly` function](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/poly) constructs orthogonal polynomials by default. Use `poly(LStat, 5, raw = T)` in R, if you want to compare to the solution with raw polynomial in R.

---

<div class="post-metadata">

**Author:** ![compleat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/compleat/32/8958_2.png) [@compleat](https://discourse.julialang.org/u/compleat)\
**Post date:** [January 30, 2022, 5:31pm UTC](https://discourse.julialang.org/t/glm-inconsistency-with-r-on-islr-dataset/75427/3 "2022-01-30T17:31:11Z")

</div>

Thank you so much, that is exactly what I needed. My class will appreciate. Presumably the overall fit (i.e., R^2 and F) will be the same. I was trying to figure out an easy alternative to poly(), but I couldn’t come up with (or find) one.