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