Lab8 - Regresja liniowa cz. 2

Diagnostyka modelu liniowego

Regresja liniowa opiera się na kilku założeniach, które łatwo zapamiętać jako zasadę LINE:

• Linear trend, liniowa zależność między $Y$ a $X$,
• Independent residuals, niezależne residua,
• Normally distributed residuals, normalny rozkład residuów,
• Equal variance of residuals, homoskedastyczność wariancji (czyli niezależność od $X$ ani $Y$).

Regresja liniowa natomiast nie zakłada, że każdej wartości $Y$ odpowiada tylko jedna wartość $X$, ani na odwrót. Jest wręcz pożądane, żeby dla jednej wartości zmiennej objaśniającej $X$ otrzymać kilka niezależnych pomiarów zmiennej objaśnianej $Y$ (na przykład zmierzyć kilkukrotnie tempo reakcji chemicznej w tej samej temperaturze), ponieważ daje nam to więcej informacji o błędach losowych.
Jest to bardzo częste źródło błędów u studentów (i nie tylko), i spotkałem się wielokrotnie z uśrednianiem wartości $Y$ odpowiadających tym samym wartościom $X$. Nie ma żadnej potrzeby, żeby tak robić, a wręcz może to prowadzić do błędnych wyników.

Oczywiście, tak samo jak w przypadku testów statystycznych, w praktyce nic nie jest idealne i poza nielicznymi sytuacjami żadne z tych założeń tak naprawdę nie jest spełnione. Im bliżej jesteśmy do ich spełnienia, tym bardziej wiarygodne wyniki otrzymamy.

Odstępstwa od poszczególnych założeń zasady LINE mają bardzo różny wpływ na nasze wnioski.
W bardzo zgrubny sposób podsumowuje je poniższa tabela. Przy pierwszym czytaniu będzie prawdopodobnie bardzo niejasna, ale mam nadzieję, że po przerobieniu zadań będzie lepiej.

Założenie	Konsekwencja odstępstwa	Co można z tym zrobić
Liniowość trendu	To podstawowe założenie całej metody, więc kiedy nie jest spełnione, to jest dosyć słabo. Parametry $\beta$ tracą w zasadzie jakiekolwiek znaczenie, tzn. przestają oznaczać średnią zmianę $Y$ gdy zmienimy $X$ o jednostkę.	Na szczęście istnieje dużo metod przekształcenia danych, które często pozwalają uzyskać liniowy trend.
Niezależność błędów	Wszystko zależy od tego, w jaki sposób błędy są zależne. Może się popsuć wszystko albo nic.	Wszystko zależy od struktury błędów. Dla pewnych specyficznych przypadków istnieją specjalne techniki przekształcania danych, na przykład analizowanie różnic pomiędzy kolejnymi obserwacjami.
Normalność błędów	Dopóki wartość średnia błędów wynosi 0, to estymacja $\beta$ jest poprawna. Nie działają testy statystyczne ani przedziały ufności.	Czasami udaje się odpowiednio przekształcić dane, na przykład transformacją Boxa-Coxa, ale dość często po prostu trzeba się z tym pogodzić.
Homoskedastyczność	Obszary na których jest wysoka wariancja błędów mają mocniejszy wpływ na estymację $\beta$, ponieważ dużo silniej wpływają na odległość kwadratową $Y$ od $\hat{Y}$. Możemy przez to mieć większy błąd estymacji. Nie działają testy przedziały ufności.	Wszystko zależy od struktury wariancji. Bardzo często wariancja rośnie wraz ze wzrostem $X$ - wówczas może się udać przekształcić dane tak, aby uzyskać homoskedastyczność.

Jak widać, przekształcanie danych jest bardzo ważnym (i często trudnym) elementem regresji liniowej.
Możemy transformować zarówno zmienną objaśnianą $Y$, jak i zmienne objaśniające $X$.
Wbrew pozorom, w regresji liniowej daje to różne wyniki: modelowanie $\log(Y) \sim X$ to nie to samo, co modelowanie $Y \sim \exp(X)$!

Rozpatrzmy pewną odwracalną funkcję $f$ (na ogół $f(x) = log(x)$ lub $f(x) = x^\alpha$), i załóżmy, że:

• Transformujemy zmienną niezależną, czyli dopasowujemy model $Y = \beta f(X) + \epsilon$. Zmieniamy wówczas trend, ale nie zmieniamy struktury błędów losowych. Stosujemy tę transformację jeśli zależność między $X$ a $Y$ nie jest liniowa, ale pozostałe założenia modelu są spełnione.
• Transformujemy zmienną zależną, czyli dopasowujemy model $f(Y) = \beta X + \epsilon$. Jest to równoważne dopasowaniu $Y = f^{-1}(\beta X + \epsilon)$. W zależności of postaci funkcji $f$, zmieniamy w ten sposób zarówno trend, jak również strukturę błędów losowych: ich rozkład, heteroskedastyczność itp.

Nie ma żadnego gotowego przepisu w jaki sposób przekształcić dane, i często trzeba posługiwać się intuicją lub wiedzą dziedzinową. Zilustrujemy teraz o co w tym wszystkim chodzi na prostym przykładzie.

Zadanie 1. Zadanie przykładowe. Załaduj dane bats.csv, znajdujące się na stronie przedmiotu. Dane zawierają długości ramion nietoperzy oraz ich wiek w momencie pomiaru. Zbadaj zależność pomiędzy wiekiem nietoperza a długością jego ramienia. W jakim wieku ramię nietoperza jest rozwinięte w 90%?

Rozwiązanie. W pierwszym kroku, jak zawsze, wczytujemy dane i wizualizujemy je używając biblioteki ggplot2. Zmienna $Y$ będzie oznaczać długość ramienia w milimetrach, a $X$ wiek nietoperza w dniach.

library(ggplot2)

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

bats <- read.csv('bats.csv', header=T)
ggplot(bats, aes(x=days, y=forearm.length)) + geom_point() + theme_minimal()

Widzimy, że zależność nie jest liniowa. Wygląda za to z grubsza jak $Y = \log(X)$. Zobaczmy zatem, co się stanie, jeśli przestawimy $Y$ w zależności od $\log(X)$:

ggplot(bats, aes(x=log(days), y=forearm.length)) + geom_point() + theme_minimal()

Zależność po transformacji wygląda na bardziej liniową niż wcześniej. Transformacja logarytmiczna ma dodatkowo tę zaletę, że jest łatwo interpretowalna: długość ramienia zwiększy się o 1 mm, gdy wiek wzrośnie około 2.7 raza, ponieważ $log(X\codot e) = log(X) + 1 = Y + 1$.

Inną popularną transformacją jest transformacja Boxa-Coxa, w której dla pewnego parametru $\lambda$ robimy przekształcenie $$Y_i \rightarrow \frac{Y_i^\lambda - 1}{\lambda}$$ lub $$X_i \rightarrow \frac{X_i^\lambda - 1}{\lambda}.$$ Dla $\lambda \rightarrow 0$ dostajemy przekształcenie logarytmiczne, tzn $Y_i \rightarrow \log(Y_i)$. Parametr $\lambda$ jest dobrany tak, aby rozkład residuów był jak najbliższy rozkładowi normalnemu.

Funkcja boxcox z pakietu MASS wykonuje transformację Boxa-Coxa zmiennej zależnej $Y$. Jej składnia jest bardzo podobna do funkcji lm():

library(MASS)
bc <- boxcox(forearm.length~days, data=bats, lambda=seq(0, 4, by=1/10))

Powyższy wykres przedstawia wartość logarytmu funkcji wiarygodności residuów (która ocenia podobieństwo do rozkładu normalnego) w zależności od parametru $\lambda$. Optymalny parametr wybieramy następująco:

best_lambda <- bc$x[which.max(bc$y)]

Następnie przekształcamy naszą zmienną i wizualizujemy wynik:

bats$forearm_boxcox <- (bats$forearm.length**best_lambda - 1)/best_lambda
ggplot(bats, aes(x=days, y = forearm_boxcox)) + geom_point() + theme_minimal()

Wygląda na to, że transformacja Boxa-Coxa na zmiennej objaśnianej nie dała w tym przypadku zbyt dobrych rezultatów. Oczywiście, podobne podejście można stosować również do transformacji logarytmicznej, i badać log-wiarygodność na przykład dla przekształcenia postaci $X \rightarrow \log(X + \lambda)$. Wymaga to jednak większego nakładu pracy. Niejednokrotnie nad jednym modelem liniowym można spędzić kilka tygodni.

Żeby tego uniknąć, na potrzeby tych zajęć zbudujemy model liniowy w oparciu o transformację logarytmiczną, a następnie przeprowadzimy diagnostykę modelu, aby ocenić, na ile odstępujemy od założeń LINE.

bats$log_days <- log(bats$days)
model <- lm(forearm.length ~ log_days, data=bats)
summary(model)

## 
## Call:
## lm(formula = forearm.length ~ log_days, data = bats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1069 -1.8870  0.2385  2.0376  5.0147 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.7853     0.8988   10.89 6.12e-13 ***
## log_days      6.7094     0.3791   17.70  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.745 on 36 degrees of freedom
## Multiple R-squared:  0.8969, Adjusted R-squared:  0.8941 
## F-statistic: 313.2 on 1 and 36 DF,  p-value: < 2.2e-16

Współczynnik $R^2$ jest bardzo wysoki, więc model wyjaśnia dużo wariancji w danych. Na podstawie współczynników możemy uznać, że $Y \approx 9.78 + 6.70\log(X)$. To oznacza, że średnia długość ramienia jednodniowego nietoperza ($X=1$) wynosi około 9.78 mm, następnie rośnie o 6.70 mm gdy jego wiek zwiększa się $e$-krotnie ($e \approx 2.72$). Zwizualizujmy nasz model:

predykcja <- data.frame('days'=bats$days, 'forearm.length' = 9.78 + 6.70*log(bats$days))
ggplot() + geom_point(aes(x=days, y=forearm.length), data=bats) + geom_line(aes(x=days, y=forearm.length), data=predykcja, col='red') + theme_minimal()

Przyjrzymy się teraz wykresom diagnostycznym, które otrzymamy wywołując komendę plot(model). Po wpisaniu tej komendy, wciśnięcie klawisza Enter wyświetli kolejne wykresy. Możemy też skorzystać z argumentu which, który kontroluje, który wykres należy wyświetlić. Więcej informacji możemy uzyskać w dokumentacji funkcji plot.lm, którą uzyskamy wywołując komendę ?plot.lm.

plot(model, which=1)

Powyższy wykres przedstawia zależność $Y$ od $\hat{Y}$, czyli przybliżenia $Y$ funkcją liniową (czyli predykcji $Y$). Służy on do oceny, czy zależność w danych jest liniowa. Na podstawie wykresu można uznać, że zależność jest w przybliżeniu liniowa, chociaż w praktyce prawdopodobnie szukalibyśmy lepszej transformacji.

Ten wykres pozwala również pobieżnie ocenić czy błędy są niezależne. Częstym źródłem zależności pomiędzy błędami jest autokorelacja - czyli korelacja $\varepsilon_{i+1}$ z $\varepsilon_{i}$ lub wcześniejszymi błędami. Na ogół interesuje nas korelacja dodatnia. Sprawia ona, że różnica pomiędzy dwoma sąsiadującymi residuami jest mniejsza niż pomiędzy dwoma wybranymi losowo. Takie zjawisko widać zresztą na powyższym wykresie. W tym przypadku jest ono jednak spowodowane tym, że modelowana przez nas zależność nie jest do końca liniowa. W zadaniach dodatkowych znajduje się konkretny przykład danych w których występuje typowa autokorelacja błędów.

plot(model, which=2)

Powyższy wykres to wykres kwantylowy residuów. Wydaje mi się, ze na tym etapie wiecie już, do czego taki wykres służy. Jak widać, pomimo tego, że trend nie jest do końca liniowy, to rozkład residuów jest bardzo zbliżony do normalnego. Wynika to stąd, że prawdziwe wartości $Y$ raz są większe, a raz mniejsze niż $\hat{Y}$, w związku z czym raz mamy nieco więcej dodatnich, a raz nieco więcej ujemnych residuów, a średnio wychodzi na zero.

plot(model, which=3)

Powyższy wykres przedstawia pierwiastek z wartości bezwzględnych wystandaryzowanych residuów (czyli $\sqrt{| \varepsilon_i - \bar{\varepsilon} |/\text{sd}(\varepsilon))}$) w zależności od predykcji $\hat{Y}$. Służy do oceny homoskedastyczności wariancji, czyli tego, czy zależy ona od wartości zmiennej objaśnianej $Y$. Tam, gdzie błąd losowy ma większą wariancję, pierwiastek ten będzie przyjmował większe wartości.

Czerwona linia na wykresie obrazuje trend. Widać, że wariancja błędu nieznacznie zmniejsza się wraz ze wzrostem $\hat{Y}$ (a więc, w naszym przypadku, wraz ze wzrostem $X$). Im dłuższe ramię nietoperza przewiduje model, tym mniejszy będzie rozrzut wokół wartości średniej. Widać to również na wykresie przedstawiającym zależność długości ramienia od wieku - dla młodych nietoperzy punkty są bardziej rozstrzelone wokół swojej wartości średniej niż dla starych.

Przy interpretacji trendu w tego typu wykresach należy zawsze uważać na obszary dla których mamy mało danych. Trend, czyli czerwona linia na wykresie, jest estymowany lokalnie na podstawie kilku najbliższych punktów. Jeśli takich punktów jest mało, to estymacja trendu jest wysoce niepewna.

plot(model, which=5)

Ostatni wykres służy do wykrywania obserwacji odstających (outlierów). Z obserwacjami odstającymi spotkaliśmy się na trzecich laboratoriach. Są to obserwacje, które pochodzą z innego rozkładu, niż nasze dane. To oznacza, że jakieś nieznane czynniki sprawiają, ze obserwacje te mają nietypowe wartości. W przypadku regresji liniowej należy na nie bardzo uważać, ponieważ mogą bardzo zaburzyć estymację parametrów $\beta$ (wyobraź sobie, że dodajemy do modelu 1-dniowego nietoperza o dwumetrowym ramieniu).

Wykres przedstawia zależność wystandardyzowanych residuów $(\varepsilon_i - \bar{\varepsilon})/\text{sd}(\varepsilon)$ od tak zwanej dźwigni. Dźwignia została omówiona na wykładzie 7. Mierzy to, jak bardzo zmieniły by się wartości estymowanych parametrów, gdybyśmy usunęli daną obserwację. Obserwacje o wysokiej dźwigni i nietypowych wartościach residuów możemy uznać za obserwacje odstające które zaburzają nasz model.

Na wykresie powyżej obserwacje odstające to na ogół pojedyncze, bardzo nietypowo położone punkty. Zazwyczaj są one położone poza obszarem oddzielonym czerwoną przerywaną linią (w tym przypadku jest ona niewidoczna). Możemy uznać, że w danych nie ma obserwacji odstających.

To podsumowuje naszą przygodę z nietoperzami. Utworzyliśmy model, który w dość grubym przybliżeniu wyjaśnia zależność długości ramienia nietoperza od jego wieku. Diagnostyka wykazała, że głównym problemem tego modelu jest niedokładnie oddany trend, oraz lekka heteroskedastyczność, objawiająca się tym, że zmienność długości ramienia obniża się wraz z wiekiem.

Niestety, pomimo tego, że residua mają rozkład normalny, ze względu na silne odstępstwo od liniowości nie możemy zbytnio ufać testom ani przedziałom ufności. Warto to zilustrować, robiąc przedział ufności dla predykcji:

przedzialy <- predict(model, data.frame('log_days'=log(1:50)), interval='prediction')
# Poniewaz funkcja predict zwraca macierz, a ggplot wymaga data.frame, to rzutujemy zmienną na odpowiedni typ:
przedzialy <- data.frame(przedzialy)
przedzialy$days <- 1:50
ggplot() + geom_point(aes(x=days, y=forearm.length), data=bats) + geom_line(aes(x=days, y=forearm.length), data=predykcja, col='red') + geom_ribbon(aes(x=days, ymin=lwr, ymax=upr), data=przedzialy, alpha=0.2, fill='orange') +  theme_minimal() 

Jak widać, obserwacje raz się koncentrują w górnej części przedziału ufności (dla nietoperzy noworodków), a raz w dolnej (dla nietoperzy młodocianych). W tym przypadku na szczęście tracimy jedynie na dokładności nasze predykcji, ponieważ przedziały są szersze niż by mogły być. Gorsza sytuacja występuje, gdy zależność jest na tyle nieliniowa, że w pewnych obszarach dużo obserwacji wychodzi poza przedział ufności. Wówczas stosowanie takiego modelu mogłoby być bardzo ryzykowne - na przykład gdyby modelował on ryzyko inwestycji na giełdzie.

Więcej na temat transformacji Boxa-Coxa i przykład jej zastosowania na innych danych można obejrzeć tutaj.

Materiał ponadprogramowy dla osób zainteresowanych modelowaniem matematycznym: Transformację logarytmiczną wybraliśmy, ponieważ wizualnie pasuje do zależności na wykresie i daje dość dobry wynik.
Takie podejście można określić jako modelowanie oparte na danych, lub data driven development.

Z drugiej strony, możemy zastanowić się, jak taka zależność powinna wyglądać w rzeczywistości. Ramię nietoperza do pewnego wieku powinno rosnąć dość szybko, następnie wzrost powinien stopniowo zwalniać aż w końcu się zatrzymać w momencie gdy nietoperz osiąga dojrzałość. Tego typu własności ma funkcja logistyczna. Możemy zatem utworzyć prosty model matematyczny wzrostu ramienia nietoperza, który zakłada, że $Y = Y_{max}\exp(aX + b)/(1 + \exp(aX + b))$, gdzie $Y_{max}$ to średnia długość ramienia dorosłego nietoperza, parametr $a$ wpływa na tempo wzrostu ramienia, a parametr $b$ wpływa na długość ramienia nietoperza noworodka (dla którego $X=0$). Takie podejście można określić jako modelowanie wychodzące od pierwszych przesłanek (modelling based on first principles) lub modelowanie ab initio.

Niestety, dokładne wyznaczenie stałych w takim modelu wymaga zastosowania metod optymalizacji numerycznej które wykraczają poza ramy tego wykładu. Z tego powodu nie omówimy ich dokładnie. Przedstawię jedynie sposób dopasowania modelu do danych w R i wyniki, jakie można za jego pomocą otrzymać, w celu porównania ich z regresją liniową:

library(stats)
logistic_model <- nls(forearm.length~SSlogis(days, Asym, xmid, scal), bats)
# Wartosci parametrow:
# Ymax = Asym, xmid = -b/a, scal = a
# Stad, Ymax = 34.434, a = 6.044, b = -27.089
predykcja  <- data.frame('forearm.length' = predict(logistic_model, 
                                                    data.frame('days'=1:50)
                                                    )
                         )
predykcja$days <- 1:50
ggplot() + geom_point(aes(x=days, y=forearm.length), data=bats) + geom_line(aes(x=days, y=forearm.length), data=predykcja, col='red') + theme_minimal()

Jak widać, model wyprowadzony z pierwszych przesłanek bardzo dobrze pasuje do danych. Konstrukcją takich modeli zajmuje się dziedzina zwana modelowaniem matematycznym. Wadą takiego modelu jest natomiast to, że wszystkie narzędzia statystyczne taką jak testy i przedziały ufności trzeba wyprowadzać na nowo. Ceną za nieco lepsze dopasowanie do danych może być zatem kilka dodatkowych lat pracy.

Diagnostyka sportowców

Zadanie 2. Załaduj zbiór danych dotyczący sportowców, który analizowaliśmy na poprzednich zajęciach. Utwórz model liniowy wyjaśniający zmienność wagi (Wt) w zależności od wysokości (Ht) i przeprowadź diagnostykę modelu. Czy założenia pozwalające na sensowną estymację parametrów $\beta$ są spełnione? Co z założeniami, które pozwalają na ocenę tego, czy na podstawie danych możemy stwierdzić, że $\beta \neq 0$? Patrząc na odstępstwa od założeń LINE, jakich efektów możemy się spodziewać, jeśli chodzi o przedziały ufności dla predykcji? Zweryfikuj swoje przypuszczenia, tworząc wykres zależności wagi od wysokości, na którym zaznaczysz prostą regresji i przedziały ufności.

Zadanie 3. Ćwiczenie ggplot2. Korzystając z wyników poprzedniego zadania, utwórz wykres punktowy zmiennej Ht w zależności od Wt, w którym: - Zamiast punktów, wykreślisz numery obserwacji, - Zaznaczysz linię reprezentującą predykcję modelu, - Zaznaczysz przedział ufności dla predykcji, - Dla każdej obserwacji poprowadzisz linię od prawdziwej wartości zmiennej Ht do wartości predykowanej.

Wykorzystaj fakt, że funkcja lm zwraca wartości predykowane w polu o nazwie fitted.values. Do stworzenia wykresu wykorzystaj geometrie geom_segment, geom_text (lub geom_label lub geom_shadowtext z pakietu shadowtext) oraz geom_ribbon. Zwróć uwagę na to, że kolejność dodawania warstw ma znaczenie: każda kolejna jest nakładana na poprzednie.

Wykresy z zaznaczonymi numerami obserwacji są bardzo przydatne gdy chcemy zidentyfikować obserwacje odstające.

Zadanie 4. Utwórz model wyjaśniający zmienność procentowej zawartości tkanki tłuszczowej w ciele (zmienna X.Bfat) za pomocą wszystkich pozostałych zmiennych, wliczając w to zmienne jakościowe (wykorzystaj formułę X.Bfat ~ .). Następnie: - Przeprowadź diagnostykę modelu. Jakie są główne odstępstwa od zasady LINE? Czy w danych znajdują się obserwacje odstające? - Sprawdź wyniki funkcji summary. Które zmienne mają znaczenie dla przewidywania procentu tkanki tłuszczowej?

Zwróć szczególną uwagę na zmienną Sport i spróbuj zinterpretować wynik. Pamiętaj, że parametr $\beta$ odpowiadający danemu poziomowi zmiennej jakościowej oznacza zmianę zmiennej zależnej, która ma miejsce wtedy, gdy zmienimy poziom zmiennej objaśnianej z poziomu bazowego na ten odpowiadający parametrowi $\beta$.

Zadanie 5. Przypomnij sobie informacje dotyczące błędu treningowego z Laboratorium 6 i odpowiadającego mu wykładu. Następnie:

• Utwórz model wyaśniający procent tkanki tłuszczowej, ale tym razem wykorzystaj jedynie co drugą obserwację ze zbioru danych. Możesz wykorzystać argument subset funkcji lm.
  
• Wykorzystaj funkcję predict do przewidzenia wartości zmiennej X.Bfat dla pozostałych obserwacji. Tutaj niestety nie ma argumentu subset.
  
• Oblicz błąd średniokwadratowy treningowy oraz testowy (jeśli nie pamiętasz, co to takiego, to przypomnij sobie Wykład 5).

Po wykonaniu powyższych punktów, wytrenuj nowy model, objaśniający zmienną X.Bfat za pomocą zmiennych Sport, Sex, SSF, LBM, Wt. Nowy model wytrenuj na tym samym zbiorze treningowym co poprzedni. Następnie:

• Porównaj wartości $R^2$ w pełnym i w mniejszym modelu. Czy mniejszy model wyjaśnia dużo mniej zmienności zmiennej objaśnianej niż pełny model?
• Ponownie wykorzystaj funkcję predict do otrzymania predykcji dla zbioru testowego.
• Oblicz błąd średniokwadratowy treningowy oraz testowy.

Który błąd wzrósł, a który zmalał? Dlaczego? Jak nazywa się to zjawisko?

Wskazówka 1. Żeby łatwo wybrać potrzebne podzbiory danych, warto stworzyć wektor indeksów zbioru treningowego to_train <- seq(2, nrow(ais), by=2). Wówczas zbiór treningowy wybieramy pisząc ais[to_train, ], a zbiór testowy wybieramy pisząc ais[-to_train, ].
Wskazówka 2. Żeby otrzymać predykcję w mniejszym modelu, należy wybrać z danych obserwacje należące do zbioru testowego, ale nie musimy wybierać kolumn odpowiadających zmiennym objaśniającym. Funkcja predict sama wybierze te kolumny, których potrzebuje. Wystarczy zatem napisać small_model_fit <- predict(small_model, ais[-to_train, ]).

Zadania dodatkowe.

Zadanie 7. Dyfuzja - przykład autokorelacji. W tym zadaniu rozpatrzymy cząstkę zawieszoną w płynie, która porusza się ruchem Browna z dryfem. Założymy mianowicie, że cząstka z jednej strony opada ruchem jednostajnie przyspieszonym, a z drugiej, że wysokość na jakiej się znajduje podlega chaotycznym zaburzeniom. Wysymulujemy ruch cząstki i zastanowimy się, w jaki sposób wykorzystać regresję liniową do estymacji jej przyspieszenia. Pomimo tego, że wstęp do tego zadania jest bardzo długi, to jest ono tak naprawdę bardzo proste.

Standardowy ruch Browna $B_t$ modeluje się za pomocą procesu Wienera, czyli w uproszczeniu takiego procesu, że $B(t+h) - B(t) \sim \mathcal{N}(0, h)$ (tutaj $h$ oznacza wariancję), a ponadto $B(t+h) - B(t)$ nie zależy od żadnej wcześniejszej zmiany $B(t) - B(t-h')$. Szczegółów można się dowiedzieć na kursie Wstęp do Analizy Stochastycznej lub w artykule na Wikipedii.
My założymy dodatkowo, że cząsteczka opada, tak że średnia wysokość, na jakiej znajduje się w chwili $t$, to $-at^2/2$ dla dodatniej stałej $a$.
Połączymy opadanie z ruchem Browna, zakładając, że położenie cząstki $H(t)$ spełnia zależności $H(t+h) - H(t) \sim \mathcal{N}(-ah^2/2 -ath, h)$. Ponadto założymy, że w chwili 0 mamy na pewno $H(0) = 0$. Zauważ, że wynika stąd, że $H(t) \sim \mathcal{N}(-at^2/2, t)$.

Przyjmiemy stały krok czasowy, czyli $h = \text{const.}$, tak że nasze czasy $t$ mają postać $hi$ dla $i = 0, 1, \dots, n$.
Oznaczmy dla uproszczenia położenie cząstki $i$-tym kroku jako $H_i = H(hi)$. Mamy wówczas H_{i+1} - H_i (-ah^2(i + 1/2), h)$.

Zmienne $H_1, \dots, H_n$ symulujemy w trzech prostych krokach:

• Obliczamy średnie położenia $m_i = -at^2/2 = -a(hi)^2/2$ dla $i = 1, 2, \dots, n$.
• Generujemy ciąg niezależnych zmiennych losowych $X_1, X_2, \dots, X_n \sim \mathcal{N}(0, h)$. Trzeba pamiętać, że $h$ jest tutaj wariancją, podczas gdy w funkcji rnorm() podajemy odchylenie standardowe.
• Obliczamy sumę kumulatywną $Y_i = X_1 + X_2 + \dots + X_i$. Korzystamy w tym celu z funkcji cumsum().

Zmienna $H_i$ jest równa $m_i + Y_i$ dla $i = 1, 2, \dots, n$. Zmienna $H_0$ jest równa 0.

Poprawność tej procedury wynika z faktu, że $H_{i} - H_{i-1} \sim \mathcal{N}(-ah^2(i + 1/2), h) = -ah^2(i + 1/2) + \mathcal{N}(0, h) = -ah^2(i+1/2) + X_i$. Wystarczy teraz przekształcić $H_i = (H_i - H_{i-1}) + (H_{i-1} - H_{i-2}) + \dots + (H_1 - H_0) + H_0 = -a(hi)^2/2 + X_1 + \dots + X_i$.

Wybierz dowolną wartość przyspieszenia $a$ i wysymuluj zmienne $H_t$ dla $t = 0, 0.1, 0.2, 0.3, \dots, 5$.
Zapisz zmienne $t$ oraz $H_t$ w postaci data frame. Dołącz również kolumnę zawierającą kwadrat zmiennej $t$.
Przedstaw zależność $H_t$ od $t$ na wykresie liniowym.
Wynik powinien wyglądać tak jak poniżej (dla $a = 5$).

Spróbuj teraz wyestymować przyspieszenie $a$ z otrzymanych w ten sposób danych. W tym celu zbuduj model liniowy, który wyjaśni zmienność $H_t$ za pomocą $t^2$. Uwaga: ważne, żeby tę ostatnią wartość mieć jako osobną kolumnę w danych - formuły typu H ~ t^2 mają swoje własne, bardzo specyficzne znaczenie.
Zwróć uwagę, że wyestymowany parametr $\beta$ przybliża wartość $-a/2$, wobec czego $\hat{a} = -2\hat{\beta}$. Możesz wykorzystać wiedzę o tym, że $H_0 = 0$, dopasowując model bez wyrazu wolnego. Formuła modelu bez wyrazu wolnego to H ~ 0 + t lub H ~ t - 1.

Żeby sprawdzić swoje wyniki, przedstaw otrzymaną zależność na wykresie.

Stwórz teraz przedział ufności dla parametru $\beta$. Powtórz kilkukrotnie analizę i zobacz jak z symulacji na symulację zmieniają się wartości parametru $\beta$ oraz przedziały ufności. Czy coś się nie zgadza?

Czemu wyniki są złe? Przeprowadź diagnostykę modelu. Spójrz na wykres Residuals vs Fitted i podziwiaj autokorelację błędów.

Jak temu zaradzić? W tym przypadku możemy usunąć autokorelację, robiąc model, który wyjaśni zmienność różnic $H_{i+1} - H_{i}$. Zgodnie z naszymi założeniami, takie różnice są niezależnymi zmiennymi losowymi, co gwarantuje spełnienie zasady LINE. Jest to częsta technika w analizie szeregów czasowych, czyli takich danych, które dotyczą procesu losowego obserwowanego przez pewien okres czasu.