# How to recover the formula of the polynomial regression model?

I am trying to recover the formula of my regression model. I build the polynomial regression model using glmer(optionval ~ nt1 + nt2 + nt3 + section:(nt1+ nt2+ nt3)+(nt1-1|participant.id),...). "nts" are time based natural polynomials generated using poly(..., raw=TRUE) and "section" is a categorical factor consisting of 2 levels.

I got the model summary including estimates like below:

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.04591    0.03217  -1.427    0.154
nt1            2.96992    2.70780   1.097    0.273
nt2           61.03648    1.78963  34.106  < 2e-16 ***
nt3          -29.26076    1.60500 -18.231  < 2e-16 ***
nt1:section1  -6.31923    2.67106  -2.366    0.018 *
nt2:section1   7.57970    1.46541   5.172 2.31e-07 ***
nt3:section1 -10.13101    1.41280  -7.171 7.45e-13 ***


Form what I know, I think these estimates are the βs of the model formula, but I am n̶o̶t̶ ̶s̶u̶r̶e̶ now sure if value for "section" is associated with the coding defined by "contrast()"(in my case is (-1,1)). Anyway, based on my assumption, the formula I recovered (61.03x^2-29.26x^3+6.32x-7.58x^2+10.13x^3; 61.03x^2-29.26x^3-6.32x-7.58x^2-10.13x^3) is different from the one visualized by averaging model predictions from predict().

Stats pros on the platform please give me some ideas!

Update 1:

I managed to recover the formula for the model with no random effect using the procedure I described above (bold sentence; see Fig.1). However, when there is a random effect, the formula recovered directly using the summarized parameters does not produce a similar curve as the one produced by averaging individual predict() values (see Fig.2). Fig.1 Fixed effect model. dashed lines are produced by recovered model formula. Solid lines are produced by averaging predict(). Fig.2 Mixed effect model. dashed lines are produced by recovered model formula. Solid lines are produced by averaging predict().

I looked the coef(model) and found different nt1 for each sample, so I realized that the predict() values for each sample must be generated using individual coefficients instead of the united model coefficients. Furthermore, the averaged predict() curve even has 1 more inflection point (3 inflection points) than the model (3rd order so 2 inflection points) could possibly produce.

In this case, I doubt the possibility of recovering the formula for a mixed effect model.

• In R, try using poly(..., raw=TRUE) to generate the usual polynomial terms $X^2$, $X^3$, etc. By default, R's poly() uses the argument raw=FALSE which uses "orthogonal polynomials." These give equivalent predictions $hat Y$, and they can be a bit nicer for computational reasons, but also a bit harder to interpret directly. Mar 13 at 0:45
• @civilstat Hi, I am clear about natural or orthogonal polynomials and I did use poly(..., raw=TRUE), so this is not really my confusion. How about the way to recover the formula?
– Jin
Mar 13 at 1:16
• @Jin What are the “nt” polynomials? Can you write down expressions like $x$ or $x^2$ (maybe more complicated than that)? // You seem to be changing some of the coefficient signs. Also, you seem to have missed the intercept, the “nt” term on its own, and the interaction between the categorical variable and your polynomial terms. It might be as simple as that.
– Dave
Mar 13 at 5:18
• @Dave Hi. Yeah, nt1 refer to x, nt2 - x², nt3 - x³. Sorry I just forget the signs (just edited). I thought Intercept is added by default and I leave it out because it is not significant. Other than those, do you suggest the rest of the formula is just right?
– Jin
Mar 13 at 10:06
• As I explain at stats.stackexchange.com/a/249202/919, you can always find the coefficients of polynomials or splines even when their methods are not documented. An example with splines is given at stats.stackexchange.com/a/101484/919.
– whuber
Mar 13 at 15:35