# Problems with simple linear model (R and Julia comparison)

**URL:** <https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229>\
**Category:** General Usage\
**Created:** [March 21, 2025, 3:15pm UTC](https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229 "2025-03-21T15:15:45Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [March 21, 2025, 3:15pm UTC](https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229/1 "2025-03-21T15:15:45Z")

</div>

Hi all, I’m trying to reproduce some R code. Can someone help me figure out why the difference in these polynomial regressions is so different between Julia and R? Yes, a 9th order polynomial fit is not practical; I’m trying to reproduce a figure from an R-based textbook.

Julia:

```julia
using RDatasets, DataFrames
using GLM
mcycle = dataset("MASS", "mcycle")
lm = fit(LinearModel, @formula(Accel ~ Times + Times^2 + Times^3 + Times^4 + Times^5 + Times^6 + Times^7 + Times^8 + Times^9), mcycle)

```

Gives

```julia
Accel ~ 1 + Times + :(Times ^ 2) + :(Times ^ 3) + :(Times ^ 4) + :(Times ^ 5) + :(Times ^ 6) + :(Times ^ 7) + :(Times ^ 8) + :(Times ^ 9)

Coefficients:
────────────────────────────────────────────────────────────────────────────────────────
                    Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 0.0 NaN NaN NaN NaN NaN
Times 0.0 NaN NaN NaN NaN NaN
Times ^ 2 0.0 NaN NaN NaN NaN NaN
Times ^ 3 0.0 NaN NaN NaN NaN NaN
Times ^ 4 0.0 NaN NaN NaN NaN NaN
Times ^ 5 0.0 NaN NaN NaN NaN NaN
Times ^ 6 0.0 NaN NaN NaN NaN NaN
Times ^ 7 3.64577e-10 1.39406e-9 0.26 0.7941 -2.3934e-9 3.12255e-9
Times ^ 8 -1.3951e-11 5.54462e-11 -0.25 0.8017 -1.23645e-10 9.57425e-11
Times ^ 9 1.33146e-13 5.48143e-13 0.24 0.8085 -9.51289e-13 1.21758e-12
────────────────────────────────────────────────────────────────────────────────────────

```

Whereas in R we have:

```r
> library(MASS)
> data(mcycle)
> lm(accel ~ poly(times, 9), data = mcycle)

Call:
lm(formula = accel ~ poly(times, 9), data = mcycle)

Coefficients:
    (Intercept) poly(times, 9)1 poly(times, 9)2 poly(times, 9)3  
        -25.546 164.557 131.227 -239.790  
poly(times, 9)4 poly(times, 9)5 poly(times, 9)6 poly(times, 9)7  
         -6.738 245.799 -83.906 -153.596  
poly(times, 9)8 poly(times, 9)9  
        163.064 31.879  

```

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [March 21, 2025, 3:48pm UTC](https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229/2 "2025-03-21T15:48:50Z")

</div>

The model matrix is ill-conditioned (see all the issues related to QR vs Cholesky on the GLM.jl repository), but there is also another difference: `poly` in R by default gives _orthogonal_ polynomials, not raw ones, so the estimates will be different. What happens if you do

```R
lm(accel ~ poly(times, 9, raw=TRUE), data = mcycle)

```

in R?

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [March 21, 2025, 3:59pm UTC](https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229/3 "2025-03-21T15:59:27Z")

</div>

Okay, with raw polynomials in R:

```R
> summary(lm(accel ~ poly(times, 9, raw=TRUE), data=mcycle))

Call:
lm(formula = accel ~ poly(times, 9, raw = TRUE), data = mcycle)

Residuals:
    Min 1Q Median 3Q Max 
-92.250 -18.707 -0.545 19.571 54.739 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)  
(Intercept) 3.467e+02 1.395e+02 2.486 0.0143 *
poly(times, 9, raw = TRUE)1 -2.288e+02 9.266e+01 -2.469 0.0149 *
poly(times, 9, raw = TRUE)2 4.963e+01 2.226e+01 2.230 0.0276 *
poly(times, 9, raw = TRUE)3 -4.811e+00 2.691e+00 -1.788 0.0762 .
poly(times, 9, raw = TRUE)4 2.271e-01 1.870e-01 1.214 0.2269  
poly(times, 9, raw = TRUE)5 -4.851e-03 7.926e-03 -0.612 0.5417  
poly(times, 9, raw = TRUE)6 1.025e-05 2.081e-04 0.049 0.9608  
poly(times, 9, raw = TRUE)7 1.460e-06 3.300e-06 0.442 0.6591  
poly(times, 9, raw = TRUE)8 -2.473e-08 2.892e-08 -0.855 0.3942  
poly(times, 9, raw = TRUE)9 1.283e-10 1.075e-10 1.194 0.2349  

```

In Julia, if we disable rank deficiency checks:

```julia
julia> fit(LinearModel, @formula(Accel ~ Times + Times^2 + Times^3 + Times^4 + Times^5 + Times^6 + Times^7 + Times^8 + Times^9), mcycle, dropcollinear=false)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}}}, Matrix{Float64}}

Accel ~ 1 + Times + :(Times ^ 2) + :(Times ^ 3) + :(Times ^ 4) + :(Times ^ 5) + :(Times ^ 6) + :(Times ^ 7) + :(Times ^ 8) + :(Times ^ 9)

Coefficients:
─────────────────────────────────────────────────────────────────────────────────────────
                     Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 346.708 139.457 2.49 0.0143 70.6622 622.754
Times -228.771 92.6506 -2.47 0.0149 -412.167 -45.3748
Times ^ 2 49.6301 22.2542 2.23 0.0276 5.57919 93.681
Times ^ 3 -4.81176 2.69039 -1.79 0.0762 -10.1372 0.513707
Times ^ 4 0.227133 0.186976 1.21 0.2268 -0.142975 0.59724
Times ^ 5 -0.0048522 0.00792551 -0.61 0.5415 -0.0205403 0.0108359
Times ^ 6 1.02888e-5 0.000208089 0.05 0.9606 -0.00040161 0.000422188
Times ^ 7 1.45897e-6 3.2998e-6 0.44 0.6592 -5.07279e-6 7.99072e-6
Times ^ 8 -2.47269e-8 2.8922e-8 -0.85 0.3942 -8.19762e-8 3.25224e-8
Times ^ 9 1.2828e-10 1.07472e-10 1.19 0.2349 -8.44551e-11 3.41015e-10
─────────────────────────────────────────────────────────────────────────────────────────

```

which is the same answer as in R.

---

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [March 21, 2025, 4:01pm UTC](https://discourse.julialang.org/t/problems-with-simple-linear-model-r-and-julia-comparison/127229/4 "2025-03-21T16:01:32Z")

</div>

Ah thanks. Those differences are helpful to know. Because its so impractical in real use I never investigated the details of `poly` in R.