Search code examples
rmodelingpolynomials

Model extraction with polynomial involved

I want to extract the function, which defines my model. I´m working with lme4 and the lmer function in Rstudio to get a Model. The function I´m using looks like:

model<-lmer(y~poly(x, 6)+ z +(1|d))

to show the estimates and other things:

summary(model)

I´ve assumed that the function looks like this: y~intercept+x*0.3418+x^2*-0.1.....+x^6*0.05871+z*-0.002566

if I insert values for x and z the predicted y values are far off from anything they should get to. when I draw the graph it shows realistic values.

plot(x, predict(Model))

Rplot output

summary(Modeltest)
Linear mixed model fit by REML. t-tests use Satterthwaite's method 
Formula: sum.ARM.g.cm³.ges ~ poly(Baumalter.bei.Jahrringsbildung, 6) +  
    Bonität.h100.gemessen.interpoliert + (1 | Baumnummer)

REML criterion at convergence: -8023.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.3263 -0.5221 -0.0344  0.6118  3.1398 

Random effects:
 Groups     Name        Variance  Std.Dev.
 Baumnummer (Intercept) 0.0010255 0.03202 
 Residual               0.0001138 0.01067 
Number of obs: 1307, groups:  Baumnummer, 12

Fixed effects:
                                           Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                               5.841e-01  4.858e-02  1.000e+01  12.024 2.87e-07 ***
poly(Baumalter.bei.Jahrringsbildung, 6)1  3.418e-01  1.114e-02  1.289e+03  30.675  < 2e-16 ***
poly(Baumalter.bei.Jahrringsbildung, 6)2 -1.066e-01  1.122e-02  1.289e+03  -9.506  < 2e-16 ***
poly(Baumalter.bei.Jahrringsbildung, 6)3 -1.467e-01  1.092e-02  1.289e+03 -13.441  < 2e-16 ***
poly(Baumalter.bei.Jahrringsbildung, 6)4  1.547e-01  1.081e-02  1.289e+03  14.305  < 2e-16 ***
poly(Baumalter.bei.Jahrringsbildung, 6)5 -6.797e-02  1.073e-02  1.289e+03  -6.334 3.30e-10 ***
poly(Baumalter.bei.Jahrringsbildung, 6)6  5.871e-02  1.072e-02  1.289e+03   5.474 5.27e-08 ***
Bonität.h100.gemessen.interpoliert       -2.566e-03  1.415e-03  1.000e+01  -1.814   0.0998 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) p(B..J,6)1 p(B..J,6)2 p(B..J,6)3 p(B..J,6)4 p(B..J,6)5 p(B..J,6)6
pl(B..J,6)1 -0.005                                                                  
pl(B..J,6)2 -0.003  0.079                                                           
pl(B..J,6)3 -0.001  0.040      0.057                                                
pl(B..J,6)4  0.003 -0.017      0.004      0.014                                     
pl(B..J,6)5  0.001 -0.020     -0.025     -0.011      0.002                          
pl(B..J,6)6  0.001 -0.012     -0.021     -0.020     -0.008      0.004               
Bntth100gms -0.982  0.005      0.004      0.002     -0.003     -0.001     -0.001

I hope you guys can help me getting a fitting function.

Solution

  • This is not an lme4-specific issue, it has to do with the poly() function. By default poly() creates an orthogonal polynomial; this is not necessarily an easy function to reconstruct. You might try poly(x, 6, raw=TRUE) which will give you a "raw" polynomial of the form b0+b1*x+b2*x^2+...

    this question gives the details of how to derive the coefficients for an orthogonal polynomial.