GLM - inconsistency with R on ISLR dataset?

General Usage Performance
,
1

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:

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

This give output

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.

2

I think you are running into an open GLM issue. 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 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.

1 Like
3

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.