0

I'm using R and I want to translate the results from lm() to an equation.

My model is:

 Residuals:
  Min        1Q    Median        3Q       Max 
-0.048110 -0.023948 -0.000376  0.024511  0.044190 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          3.17691    0.00909  349.50  < 2e-16 ***
poly(QPB2_REF1, 2)1  0.64947    0.03015   21.54 2.66e-14 ***
poly(QPB2_REF1, 2)2  0.10824    0.03015    3.59  0.00209 ** 
B2DBSA_REF1DONSON   -0.20959    0.01286  -16.30 3.17e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.03015 on 18 degrees of freedom
Multiple R-squared: 0.9763, Adjusted R-squared: 0.9724 
F-statistic: 247.6 on 3 and 18 DF,  p-value: 8.098e-15 

Do you have any idea?

I tried to have something like

 f <- function(x) {3.17691 + 0.64947*x +0.10824*x^2 -0.20959*1 + 0.03015^2}

but when I tried to set a x, the f(x) value is incorrect.

zinon
  • 4,427
  • 14
  • 70
  • 112
  • 1
    what do you mean by <>? – Mamoun Benghezal Jul 15 '15 at 15:17
  • Note that `poly()` returns orthogonal polynomials by default which may not be what you expect. But you really should provide a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5963610) with sample input data that we can run to make the problem more clear. As already answered, you probably should be using `predict()` – MrFlick Jul 15 '15 at 15:28
  • 1
    Try `poly(QPB2_REF1, 2, raw = TRUE)`. – Roland Jul 15 '15 at 15:29

3 Answers3

1

If you want to calculate y-hat based on the model, you can just use predict!

Example:

set.seed(123)
my_dat <- data.frame(x=1:10, e=rnorm(10))
my_dat$y <- with(my_dat, x*2 + e)
my_lm <- lm(y~x, data=my_dat)

summary(my_lm)

Result:

Call:
lm(formula = y ~ x, data = my_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1348 -0.5624 -0.1393  0.3854  1.6814 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.5255     0.6673   0.787    0.454    
x             1.9180     0.1075  17.835    1e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9768 on 8 degrees of freedom
Multiple R-squared:  0.9755,    Adjusted R-squared:  0.9724 
F-statistic: 318.1 on 1 and 8 DF,  p-value: 1e-07

Now, instead of making a function like 0.5255 + x * 1.9180 manually, I just call predict for my_lm:

predict(my_lm, data.frame(x=11:20))

Same result as this (not counting minor errors from rounding the slope/intercept estimates):

0.5255 + (11:20) * 1.9180
arvi1000
  • 9,393
  • 2
  • 42
  • 52
1

Your output indicates that the model includes use of the poly function which be default orthogonalizes the polynomials (includes centering the x's and other things). In your formula there is no orthogonalization done and that is the likely difference. You can refit the model using raw=TRUE in the call to poly to get the raw coefficients that can be multiplied by $x$ and $x^2$.

You may also be interested in the Function function in the rms package which automates creating functions from fitted models.

Edit

Here is an example:

library(rms)

xx <- 1:25
yy <- 5 - 1.5*xx + 0.1*xx^2 + rnorm(25)
plot(xx,yy)

fit <- ols(yy ~ pol(xx,2))

mypred <- Function(fit)
curve(mypred, add=TRUE)

mypred( c(1,25, 3, 3.5))

You need to use the rms functions for fitting (ols and pol for this example instead of lm and poly).

Greg Snow
  • 48,497
  • 6
  • 83
  • 110
  • Thank you for your answer. Can you please give me an example of how to use `Function` function from rms package? – zinon Jul 16 '15 at 08:40
  • @zinon, I added an example above. Your can see more examples and the help by doing `?Function.rms`. – Greg Snow Jul 16 '15 at 16:35
1

If you are looking for actually visualizing or writing out a complex equation (e.g. something that has restricted cubic spline transformations), I recommend using the rms package, fitting your model, and using the latex function to see it in latex

my_lm <- ols(y~x, data=my_dat)
latex(my_lm)

Note you will need to render the latex code so as to see your equation. There are websites and, if you are using a Mac, Mac Tex software, that will render it for you.

MT-FreeHK
  • 2,462
  • 1
  • 13
  • 29