Ajustando o modelo polinomial aos dados em R

83

Li as respostas a esta pergunta e são bastante úteis, mas preciso de ajuda principalmente na R.

Tenho um exemplo de conjunto de dados em R da seguinte forma:

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)

Quero ajustar um modelo a esses dados para que y = f(x). Eu quero que seja um modelo polinomial de 3ª ordem.

Como posso fazer isso em R?

Além disso, o R pode me ajudar a encontrar o modelo de melhor ajuste?

r curve-fitting data-analysis polynomial-math Mehper C. Palavuzlar
fonte

Respostas:

98

Para obter um polinômio de terceira ordem em x (x ^ 3), você pode fazer

lm(y ~ x + I(x^2) + I(x^3))

ou

lm(y ~ poly(x, 3, raw=TRUE))

Você poderia ajustar um polinômio de 10ª ordem e obter um ajuste quase perfeito, mas deveria?

EDIT: poly (x, 3) é provavelmente uma escolha melhor (veja @hadley abaixo).

Greg
fonte
6
está certo ao perguntar "você deveria". Os dados da amostra possuem apenas 8 pontos. Os graus de liberdade são muito baixos aqui. Os dados da vida real podem ter muito mais, é claro.
JD Long
1
Obrigado pela sua resposta. Que tal fazer com que R encontre o modelo de melhor ajuste? Existem funções para isso?
Mehper C. Palavuzlar
5
Depende da sua definição de "melhor modelo". O modelo que fornece o maior R ^ 2 (que um polinômio de 10ª ordem daria) não é necessariamente o "melhor" modelo. Os termos em seu modelo precisam ser escolhidos de maneira razoável. Você pode obter um ajuste quase perfeito com vários parâmetros, mas o modelo não terá poder preditivo e será inútil para outra coisa senão desenhar a linha de melhor ajuste através dos pontos.
Greg,
10
Por que você está usando raw = T? É melhor usar variáveis ​​não correlacionadas.
hadley
2
Fiz isso para obter os mesmos resultados que lm(y ~ x + I(x^2) + I(x^3)). Talvez não seja o ideal, apenas fornecendo dois meios para o mesmo fim.
Greg
45

Qual modelo é o "modelo de melhor ajuste" depende do que você entende por "melhor". R tem ferramentas para ajudar, mas você precisa fornecer a definição de "melhor" para escolher entre elas. Considere os seguintes dados e código de exemplo:

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')

Qual desses modelos é o melhor? argumentos poderiam ser feitos para qualquer um deles (mas eu, pelo menos, não gostaria de usar o roxo para interpolação).

Greg Snow
fonte
15

Com relação à pergunta 'R pode me ajudar a encontrar o modelo de melhor ajuste', provavelmente há uma função para fazer isso, supondo que você possa declarar o conjunto de modelos para testar, mas esta seria uma boa primeira abordagem para o conjunto de n-1 polinômios de grau:

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)

Notas

  • A validade desta abordagem dependerá de seus objetivos, as premissas de optimize() e AIC()e se AIC é o critério que você deseja usar,

  • polyfit()pode não ter um único mínimo. verifique isso com algo como:

    for (i in 2:length(x)-1) print(polyfit(i))

  • Usei a as.integer()função porque não está claro para mim como interpretaria um polinômio não inteiro.

  • para testar um conjunto arbitrário de equações matemáticas, considere o programa 'Eureqa' revisado por Andrew Gelman aqui

Atualizar

Veja também a stepAICfunção (no pacote MASS) para automatizar a seleção do modelo.

David LeBauer
fonte
Como posso fazer a interface do Eurequa com o R?
adam.888
@ adam.888 ótima pergunta - não sei a resposta, mas você poderia publicá-la separadamente. Esse último ponto foi um pouco uma digressão.
David LeBauer de
Nota: AIC é o Critério de Informação Akaike , que recompensa um ajuste estreito e penaliza um grande número de parâmetros de um modelo, de uma forma que se mostrou ótima em vários sentidos. en.wikipedia.org/wiki/Akaike_information_criterion
Evgeni Sergeev
5

A maneira mais fácil de encontrar o melhor ajuste em R é codificar o modelo como:

lm.1 <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + ...)

Depois de usar regressão AIC de redução

lm.s <- step(lm.1)
Matthew Fidler
fonte
5
Usar I(x^2), etc. não fornece polinômios ortogonais apropriados para ajuste.
Brian Diggs,