R How to construct polynom with form c0 + c1 * x + .. + cn * x^n from orthogonal lm (with poly(raw=F))

Question

How can I extract the coefficients of an orthogonal polynomial regression in R?

It works like a charm with a raw regression:

#create datas
set.seed(120)
x= 1:20
y=x^2 + rnorm(length(x),0,10)
datas = data.frame(x,y)

#compute raw model with function poly
mod= lm(data=datas, y~poly(x,2,raw=T))

#get coefficients with function coef()
coefficients = coef(mod)

#construct polynom and check fitted values
fitted_values = mod$fitted.values
x0 = datas$x[1]
solution = coefficients[1]+ coefficients[2]*x0^1 + coefficients[3]*x0^2
print(solution)
# 1.001596 
print(fitted_values[1])
# 1.001596
# 1.001596  == 1.001596

But the coefficient obtained with the function coef on an orthogonal lm does not work:

#create datas
set.seed(120)
x= 1:20
y=x^2 + rnorm(length(x),0,10)
datas = data.frame(x,y)

#compute raw model with function poly
mod = lm(data=datas, y~poly(x,2,raw=F))

#get coefficients with function coef()
coefficients = coef(mod)
fitted_values = mod$fitted.values

#construct polynom and check fitted values
x0 = datas$x[1]
solution = coefficients[1]+ coefficients[2]*x0^1 + coefficients[3]*x0^2
print(solution)
# 805.8476 
print(fitted_values[1])
# 1.001596
# 1.001596 != 805.8476 

Is there another way to get the right parameters to construct a polynom of the form c0 + c1 * x + .. + cn * x^n and use it to solve or predict?

I need to solve the equation, meaning getting the x given a y with the function base::solve.

Thank you

Answer 1

The issue is that you don't only need the coefficients but also the orthogonal polynomials (instead of the raw polynomials you are trying to use). The latter are constructed by model.matrix:

newdata <- model.matrix(~ poly(x, 2), data = datas)
solution <- newdata %*% coefficients
print(solution[1])
# [1] 1.001596 
print(fitted_values[1])
#1.001596

I don't understand the connection to solve but trust you can take it from here.