Coefficients of my polynomial model in R don't match graph

Question

Using Greg's helpful answer here, I fit a second order polynomial regression line to my dataset:

poly.fit<-lm(y~poly(x,2),df)

When I plot the line, I get the graph below:

The coefficients are:

# Coefficients:
#   (Intercept)     poly(x, 2)1     poly(x, 2)2  
#         727.1           362.4          -269.0

I then wanted to find the x-value of the peak. I assume there is an easy way to do so in R but I did not know it,* so I went to Wolfram Alpha. I entered the equation:

y=727.1+362.4x-269x^2

Wolfram Alpha returned the following:

As you can see, the function intersects the x-axis at approximately x=2.4. This is obviously different from my plot in R, which ranges from 0≤x≤80. Why are these different? Does R interpret my x-values as a fraction of some backroom variable?

*I would also appreciate answers on how to find this peak. Obviously I could take the derivative, but how do I set to zero?

Answer 1

In the case of a quadratic polynomial, you can of course use a little calculus and algebra (once you have friendly coefficients).

Somewhat more generally, you can get an estimate by evaluating your model over a range of candidate values and determining which one gives you the maximum response value.

Here is a (only moderately robust) function which will work here.

xmax <- function(fit, startx, endx, x='x', within=NA){
    ## find approximate value of variable x where model 
    ## specified by fit takes maximum value, inside interval 
    ## [startx, endx]; precision specified by within

    within <- ifelse(is.na(within), (endx - startx)/100, within)
    testx <- seq(startx, endx, by=within)
    testlist <- list(testx)
    names(testlist)[1] <- x
    testy <- predict(fit, testlist)
    testx[which.max(testy)]
}

Note if your predictor variable were called something other than x, you have to specify it as a string in the x parameter.

So to find the x value where your curve has its peak:

xmax(poly.fit, 50, 80, within=0.1)

Answer 2

Use predict.

 plot( 40:90, predict( poly.fit, list(x=40:90) )