# Using a given polynomial formula in a lm() model in R

I am currently trying to fit a polynomial model to measurement data using the lm() function.

fit_poly4 <- lm(y ~ poly(x, degree = 4, raw = T), weights = w)


with x as independent,

y as dependent variable

and w = 1/variance of the measurements

summary(fit_poly4)

Call:
lm(formula = y ~ poly(x, degree = 4, raw = T), weights = w)

Weighted Residuals:
Min       1Q   Median       3Q      Max
-1.57259 -1.02934  0.00252  0.98814  1.48758

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                      39.0375     3.4460  11.328 2.94e-13 ***
poly(x, degree = 4, raw = T)1  55.6996    17.7858   3.132 0.003501 **
poly(x, degree = 4, raw = T)2 -71.8194    19.9575  -3.599 0.000979 ***
poly(x, degree = 4, raw = T)3  23.8642     7.0456   3.387 0.001759 **
poly(x, degree = 4, raw = T)4  -2.4069     0.7507  -3.206 0.002872 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.069 on 35 degrees of freedom
Multiple R-squared:  0.8903,    Adjusted R-squared:  0.8778
F-statistic: 71.04 on 4 and 35 DF,  p-value: 2.637e-16


While this gives me nice enough results on the statistical side, the derived polynomial

y = 39.0375 + 55.6996*x - 71.8194*x^2 + 23.8642*x^3 - 2.4069*x^4

specifically the second peak between 3 and 5 (see picture) does not make sense for my specific case. Therefore I want to try and use a given polynomial formula

y = -3,3583*x^4 + 43*x^3 - 191,14*x^2 + 328,2*x - 137,7

in lm().

I tried

fit_poly4 <- lm(y ~ 328.2*x-191.14*I(x^2)+43*I(x^3)-3.3583*I(x^4)-137.3, weights = w)


which just returns

Error in terms.formula(formula, data = data) :
invalid model formula in ExtractVars


How do I need to write my polynomial in lm() to get this to work?

Alternatively is there any way to force lm() to not produce a second peak but to smoothly decrease from 3 to 5?

Edit: The data represents biomass (y) over different vol% nutrient concentrations (x) and should be continuous.

I've also tested simpler polynomials which either come back to imprecise in case of degree = 2 or have their own flaws (negative biomass) in case of degree = 3 • Polynomials are not a good choice if you're trying to fit a smooth curve through discrete sets of x-values. You might consider (for example) natural regression splines with a few prespecified knots, or perhaps smoothing splines (or a variety of other options. Another possibility would be to use orthogonal polynomials with fairly strong regularization on the higher order terms, such as ridge regression or perhaps elastic net. At the very least you would need to impose the more obvious constraints on it. – Glen_b -Reinstate Monica Apr 1 '19 at 0:25
• fitting on a log-y scale will take care of the possibility of negative predictions (also deal with some heteroscedasticity). Can we have a reproducible example? (I answered the technical question here, but it looks like there's room for some statistical questions as well ... – Ben Bolker Apr 1 '19 at 0:38