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.
poly(..., raw=TRUE)
to generate the usual polynomial terms $X^2$, $X^3$, etc. By default, R'spoly()
uses the argumentraw=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. $\endgroup$