---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Jak widać, model sześcienny jest jeszcze lepszy! Co więcej, dalsze zwiększanie stopnia wie-lomianu dalej owocuje lepszym dopasowaniem (wyników nie przytaczam ze względu na objętość opracowania). Jeszcze dla równania piątego stopnia współczynnik przy x5oraz ANOVA mają istot-ność rzędu p = 0.01. Dopiero w równaniu stopnia szóstego współczynniki stają się nieistotne, zaś test ANOVA między stopniem piątym a szóstym ma wysokie p = 0.736. W praktyce najczęściej najlepszy okazuje się model liniowy lub kwadratowy, rzadko kiedy sześcienny.
Aby narysować wykres zawierający dane oraz krzywą regresji wyższego stopnia, należy się troszkę pomęczyć. Jednak przy okazji istnieje łatwa możliwość dodania „krzywych ufności” dla błędu predykcji. Popatrzmy, jak wygląda dopasowanie liniowe i kwadratowe.
> grid=seq(0,10,.01) # tworzymy sekwencję x > cv=qt(.975,27) # kwantyl dla 95% ufności > p=predict(fit,data.frame(x=grid),se=T)
> matplot(grid,cbind(p$fit,p$fit-cv*p$se,p$fit+cv*p$se),type="l");points(x,y) > p=predict(fit2,data.frame(x=grid),se=T)
> matplot(grid,cbind(p$fit,p$fit-cv*p$se,p$fit+cv*p$se),type="l");points(x,y) Funkcja predict powoduje obliczenie wartości y dla zadanych x na podstawie równania re-gresji. W tym przypadku tworzy siatkę współrzędnych krzywej. W zmiennej p$fit jest wynik tego dopasowania, a w p$se błąd standardowy predykcji. Funkcja matplot powoduje narysowanie kilku wykresów na podstawie współrzędnych w zadanej macierzy, w tym przypadku są to wykre-sy podające y oraz y ± SEt. Na takim wykresie rysujemy punkty odpowiadające danym poprze
points.
7.3. Regresja nieliniowa
W praktyce większość danych doświadczalnych daje się wystarczająco dopasować do wielomianu kwadratowego lub sześciennego. Nie interesuje nas wtedy rzeczywista zależność między zmienny-mi, gdyż wielomian taki jest wystarczającym przybliżeniem tej zależności (co wynika z twierdzenia
Taylora). Czasami jednak zachodzi potrzeba dopasowania danych do konkretnej nieliniowej funk-cji. Umożliwia to polecenie nls. W tym przypadku dokonamy dopasowania naszych danych do równania y = xm+ b.
> fitn=nls(y~x^m+b,start=c(m=0.6,b=1.5)) # musimy zadeklarować
przybliżone wartości startowe > summary(fitn)
Formula: y ~ x^m + b Parameters:
Estimate Std. Error t value Pr(>|t|) m 0.5001452 0.0008658 577.7 <2e-16 *** b 1.0006124 0.0038929 257.0 <2e-16 ***
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.01097 on 28 degrees of freedom Correlation of Parameter Estimates:
m b -0.8575
Na pewno każdy czytelnik zauważy, że estymatory są bardzo „okrągłe”. W tym jest cały se-kret — analizowane przez nas dane specjalnie dobrałem na podstawie funkcji y = x12 + 1. Dlatego też były one tak „paskudne” w dopasowaniu do wielomianów niższego stopnia w poprzednim roz-dziale.
Pakiet R dysponuje całym zestawem funkcji dopasowujących do najczęściej spotykanych modeli nieliniowych. Zaczynają się one literami SS. Dla przykładu dopasujmy nasze dane do modelu Michaelisa-Menten:
> fitmm = nls(y ~ SSmicmen(x,V,K)) > summary(fitmm)
Formula: y ~ SSmicmen(x, V, K) Parameters:
Estimate Std. Error t value Pr(>|t|) V 4.6586 0.1215 38.36 < 2e-16 *** K 1.8435 0.1764 10.45 3.62e-11 ***
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1847 on 28 degrees of freedom Correlation of Parameter Estimates:
V K 0.9206
Do wyboru optymalnego modelu regresji służy gotowa funkcja, wyliczająca współczynnik AIC (Akaike’s Information Criterion). Jest to najczęściej wybierane kryterium optymalnego dopasowa-nia. Wartość tego współczynnika jest zależna nie tylko od sumy kwadratów reszt, ale również od ilości zmiennych w równaniu. Dlatego też zwiększając rząd wielomianu, mimo iż suma kwadratów reszt zawsze będzie maleć, od pewnego momentu współczynnik AIC zacznie rosnąć (i to będzie optymalny stopień wielomianu). Użycie tego współczynnika ma również związek z niemożnością stosowania współczynnika determinacji R2dla modeli nieliniowych. Zobaczmy zatem który model ma najmniejszy AIC:
> AIC(fit,fit2,fit3,fit4,fit5,fit6,fitn,fitmm) df AIC fit 3 -49.94341 fit2 4 -116.03595 fit3 5 -158.82781 fit4 6 -175.68776 fit5 7 -181.84331 fit6 8 -179.99476 fitn 2 -183.69302 fitmm 2 -14.27206
W przypadku regresji wielomianowej potwierdzają się wcześniejsze wyniki — wielomian 5 stop-nia ma najmniejszą wartość, a 6 stopstop-nia nieco większą. Zupełnie najmniejszą wartość ma model nieliniowy y = xm+ b (co się zgadza, w związku ze sposobem wygenerowania tych danych). Model Michaelisa-Menten okazał się jeszcze gorszy, niż zwykła regresja liniowa.
8. Co dalej?
Zapewne większość czytelników jest zaskoczona, że to już koniec. Ale to nieprawda. To dopiero początek. Niniejszy tekst zaznajomił z podstawowymi problemami dotyczącymi stosowania R i jest to wręcz wierzchołek góry lodowej. Dlatego na zakończenie doszedłem do wniosku, że trzeba napisać kilka słów dotyczących zaawansowanej nauki środowiska R.
Autorzy R podkreślają wielokrotnie, że powinno się nazywać ich program nie „pakietem”, co sugerowałoby zamknięty program do konkretnego celu, ale „ językiem”, bądź „środowiskiem”. W pełni na to miano zasługuje. W środowisku R można tworzyć własne polecenia i funkcje, co zdecydowanie upraszcza wykonywanie skomplikowanych, a jednocześnie rutynowych obliczeń. R posiada wszystkie struktury znane z innych języków, jak pętle, instrukcje warunkowe etc. Jest więc całkowicie otwarty, a liczba dostępnych pakietów i rozszerzeń świadczy o możliwości zaadaptowania do wielu zaawansowanych zadań.
Pisząc to opracowanie długo zastanawiałem się nad tym, czy umieszczać w nim podstawy pro-gramowania w R. Po długich namysłach zrezygnowałem z tego. Chciałem pokazać, że podstawowe obliczenia statystyczne można wykonać bez znajomości programowania w R, nawet bez znajomości jakiegokolwiek oprogramowania. Natomiast każdy, kto w życiu programował, nauczy się struktur programistycznych R w ciągu jednego popołudnia z oryginalnej dokumentacji.
Poniżej przedstawiam interesujące dokumenty anglojęzyczne dotyczące R, dostępne bezpłatnie w sieci:
1. „An introduction to R” - to oficjalna dokumentacja środowiska, która jest doskonałym kom-pendium teoretycznym. Omawia wszystkie struktury środowiska R, programowanie, oraz import i eksport danych.
2. „Using R for Data Analysis and Graphics - An introduction” (J.H.Maindonald) — jest doskonałym uzupełnieniem oficjalnej dokumentacji. Opisuje zastosowanie R do różnych ty-powych zadań statystycznych.
3. „simpleR - Using R for Introductory Statistics” (J.Verzani) —bardzo obszerne opracowanie dotyczące podstaw statystyki z użyciem R. Zawiera wyczerpujące przykłady oraz zadania do samodzielnego wykonania.
4. „R for beginners” (E.Paradis) — krótkie opracowanie omawiające podstawowe funkcje aryt-metyczne, graficzne oraz programowanie.
5. „Practical Regression and ANOVA using R” (J.Faraway) — bardzo obszerny tekst o róż-nych metodach regresji i analizy wariancji, z konkretnymi przykładami.
6. Dodatki do dzieła „An R and S-PLUS Companion to Applied Regression” Johna Foxa. Sama praca jest w sieci niedostępna i należy ją kupić. Dodatki są natomiast do pobrania. Dotyczą regresji nieparametrycznej, nieliniowej, regresji danych zmiennych w czasie etc. Mam nadzieję, że mój krótki tekst oraz powyższe dzieła przyczynią się do szybkiego „zaprzy-jaźnienia” z R każdego użytkownika. Życzę tego każdemu, aby szybko docenił możliwości pakietu i stosował go w codziennej pracy.