I've read the answers to this question and they are quite helpful, but I need help particularly in R.
I have an example data set in R as follows:
x <- c(32,64,96,118,126,144,152.5,158)
y <- c(99.5,104.8,108.5,100,86,64,35.3,15)
I want to fit a model to these data so that y = f(x). I want it to be a 3rd order polynomial model.
How can I do that in R?
Additionally, can R help me to find the best fitting model?
This question is related to
r
curve-fitting
data-analysis
polynomial-math
Regarding the question 'can R help me find the best fitting model', there is probably a function to do this, assuming you can state the set of models to test, but this would be a good first approach for the set of n-1 degree polynomials:
polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)
Notes
The validity of this approach will depend on your objectives, the assumptions of optimize() and AIC() and if AIC is the criterion that you want to use,
polyfit() may not have a single minimum. check this with something like:
for (i in 2:length(x)-1) print(polyfit(i))
I used the as.integer() function because it is not clear to me how I would interpret a non-integer polynomial.
for testing an arbitrary set of mathematical equations, consider the 'Eureqa' program reviewed by Andrew Gelman here
Update
Also see the stepAIC function (in the MASS package) to automate model selection.
The easiest way to find the best fit in R is to code the model as:
lm.1 <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + ...)
After using step down AIC regression
lm.s <- step(lm.1)
Which model is the "best fitting model" depends on what you mean by "best". R has tools to help, but you need to provide the definition for "best" to choose between them. Consider the following example data and code:
x <- 1:10
y <- x + c(-0.5,0.5)
plot(x,y, xlim=c(0,11), ylim=c(-1,12))
fit1 <- lm( y~offset(x) -1 )
fit2 <- lm( y~x )
fit3 <- lm( y~poly(x,3) )
fit4 <- lm( y~poly(x,9) )
library(splines)
fit5 <- lm( y~ns(x, 3) )
fit6 <- lm( y~ns(x, 9) )
fit7 <- lm( y ~ x + cos(x*pi) )
xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1, data.frame(x=xx)), col='blue')
lines(xx, predict(fit2, data.frame(x=xx)), col='green')
lines(xx, predict(fit3, data.frame(x=xx)), col='red')
lines(xx, predict(fit4, data.frame(x=xx)), col='purple')
lines(xx, predict(fit5, data.frame(x=xx)), col='orange')
lines(xx, predict(fit6, data.frame(x=xx)), col='grey')
lines(xx, predict(fit7, data.frame(x=xx)), col='black')
Which of those models is the best? arguments could be made for any of them (but I for one would not want to use the purple one for interpolation).
Source: Stackoverflow.com