Dopasowanie modelu wielomianowego do danych w R.

86

Przeczytałem odpowiedzi na to pytanie i są one bardzo pomocne, ale potrzebuję pomocy szczególnie w R.

Mam przykładowy zestaw danych w R w następujący sposób:

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)

Chcę dopasować model do tych danych, żeby y = f(x). Chcę, żeby był to model wielomianowy trzeciego rzędu.

Jak mogę to zrobić w R?

Ponadto czy R może mi pomóc znaleźć najlepiej dopasowany model?

Mehper C. Palavuzlar
źródło

Odpowiedzi:

102

Aby uzyskać wielomian trzeciego rzędu w x (x ^ 3), możesz to zrobić

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

lub

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

Możesz dopasować wielomian dziesiątego rzędu i uzyskać prawie idealne dopasowanie, ale czy powinieneś?

EDYCJA: poly (x, 3) jest prawdopodobnie lepszym wyborem (patrz @hadley poniżej).

Greg
źródło
6
jest na miejscu w pytaniu „czy powinieneś”. Przykładowe dane mają tylko 8 punktów. Stopnie swobody są tutaj dość niskie. Rzeczywiste dane mogą mieć oczywiście znacznie więcej.
JD Long,
1
Dzięki za odpowiedź. A co powiesz na to, żeby R znalazł najlepiej dopasowany model? Czy są do tego jakieś funkcje?
Mehper C. Palavuzlar
5
To zależy od Twojej definicji „najlepszego modelu”. Model, który daje największe R ^ 2 (który dałby wielomian dziesiątego rzędu) niekoniecznie jest „najlepszym” modelem. Warunki w Twoim modelu muszą być rozsądnie wybrane. Możesz uzyskać prawie idealne dopasowanie z wieloma parametrami, ale model nie będzie miał mocy predykcyjnej i będzie bezużyteczny do niczego innego niż rysowanie linii najlepszego dopasowania przez punkty.
Greg,
11
Dlaczego używasz raw = T? Lepiej jest używać nieskorelowanych zmiennych.
hadley,
2
Zrobiłem to, aby uzyskać takie same wyniki, jak lm(y ~ x + I(x^2) + I(x^3)). Być może nie optymalne, po prostu dając dwa środki do tego samego celu.
Greg,
46

To, który model jest „najlepiej dopasowanym modelem”, zależy od tego, co rozumiesz przez „najlepszy”. R ma narzędzia, które mogą Ci pomóc, ale aby dokonać wyboru, musisz podać definicję „najlepszego”. Rozważmy następujące przykładowe dane i kod:

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

Który z tych modeli jest najlepszy? można by argumentować dla każdego z nich (ale ja nie chciałbym używać fioletowego do interpolacji).

Greg Snow
źródło
15

Jeśli chodzi o pytanie `` czy R może mi pomóc znaleźć najlepiej dopasowany model '', prawdopodobnie istnieje funkcja, która to zrobi, zakładając, że możesz określić zestaw modeli do przetestowania, ale byłoby to dobre pierwsze podejście dla zbioru n-1 wielomiany stopnia:

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

Uwagi

  • Słuszność tego podejścia zależy od celów, założenia optimize()i AIC()jeśli AIC jest kryterium, które chcesz użyć,

  • polyfit()może nie mieć jednego minimum. sprawdź to za pomocą czegoś takiego:

    for (i in 2:length(x)-1) print(polyfit(i))
    
  • Użyłem tej as.integer()funkcji, ponieważ nie jest dla mnie jasne, jak zinterpretować wielomian niecałkowity.

  • aby przetestować dowolny zestaw równań matematycznych, rozważ program „Eureqa” omówiony przez Andrew Gelmana tutaj

Aktualizacja

Zobacz także stepAICfunkcję (w pakiecie MASS) do automatyzacji wyboru modelu.

David LeBauer
źródło
Jak mogę połączyć Eurequa z R?
adam.888,
@ adam.888 świetne pytanie - nie znam odpowiedzi, ale możesz je opublikować osobno. Ten ostatni punkt był trochę dygresją.
David LeBauer,
Uwaga: AIC to kryterium informacyjne Akaike'a , które nagradza za ścisłe dopasowanie i penalizuje większą liczbę parametrów modelu, w sposób, który okazał się optymalny pod różnymi względami. en.wikipedia.org/wiki/Akaike_information_criterion
Evgeni Sergeev
4

Najłatwiejszym sposobem znalezienia najlepszego dopasowania w R jest zakodowanie modelu jako:

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

Po zastosowaniu regresji obniżania AIC

lm.s <- step(lm.1)
Matthew Fidler
źródło
6
Użycie I(x^2)itp. Nie daje odpowiednio ortogonalnych wielomianów do dopasowania.
Brian Diggs