Zadanie 1

Zestaw danych “Prestige” (z pakietu “car”) zawiera dane nt. prestiżu n=102 Kanadyjskich zawodów z 1971 roku, a także średni dochód w danym zawodzie. Do zbadania zależności między prestiżem a dochodem wykorzystaj metody regresji nieparametrycznej.

Najpierw załaduj dane i zwizualizuj relację pomiędzy dochodem (X) a prestiżem (Y).

data("Prestige")
ggplot(Prestige)+
  geom_point(aes(x=income, y=prestige))

Związek wygląda na nieliniowy. Dla zawodów, które zarabiają mniej niż $10K, istnieje silna (pozytywna) liniowa zależność pomiędzy dochodem a prestiżem. Jednak w przypadku zawodów, które zarabiają od 10 do 25 tysięcy dolarów, związek ten ma znacznie inne (osłabione) nachylenie.

1. Kernel smoothing

fit <- locpoly(x = Prestige$prestige, y = Prestige$income, degree = 0, bandwidth = 1, gridsize = 10000)
fit <- tibble(x = fit$x, y = fit$y)

fit2 <- locpoly(x = Prestige$prestige, y = Prestige$income, degree = 0, bandwidth = 5)
fit2 <- tibble(x = fit2$x, y = fit2$y)

ggplot(data = Prestige, aes(x = prestige, y = income)) +
  geom_point() +
  geom_line(data = fit, aes(x = x, y = y), col = 'orange') +
  geom_line(data = fit2, aes(x = x, y = y), col = 'green') +
  ggtitle("Kernel smoothing - szerokość pasma 1 vs 5")

Pomarańczowa linia, przedstawiająca szerokość pasma 4 wydaje się być zbyt szczegółowa, tzw. overfitting. Zielona linia bandwidth=4 wydaje się być bardiej odpowiednia.

PS Zamieniliśmy miejscami na wykresie x=prestige, y=income, ze względnu na niespodziewane, nieporządane łączenie punktów na wykresie.

2. Loess - wielomian lokalnie kwadratowy (dagree=2)

fit3 <- loess(income ~ prestige, span=0.25, degree=2,
  family="gaussian", Prestige)
fit4 <- loess(income ~ prestige, span=0.75, degree=2,
  family="gaussian", Prestige)

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(aes(x=prestige, y=fitted(fit3)), col='pink')+
  geom_line(aes(x=prestige, y=fitted(fit4)), col='darkred')+
  ggtitle("Loess 0.25 vs 0.75")

Na powyższym wykresie bordowa linia, przedstawiająca szerokość pasma na poziomie 0.75 wydaje się być zbyt wygładzona. Zwłaszcza przy poziomie “prestige” powyżej 65, gdzie “nie łapie” obserwacji daleko od siebie odsuniętych. Za tą linią, dochody wachają się aż od ok 18tys. dolarów, do zaledwie 5tys. dolarów.

Loess z przedziałem ufności

ggplot(Prestige) +
    geom_point(aes(x=prestige,y=income)) +
    geom_smooth(aes(x=prestige,y=income), method='loess',
      span=0.25)
## `geom_smooth()` using formula = 'y ~ x'

3. Sploty interpolujące

fit5 <- smooth.spline(Prestige$income, Prestige$prestige, cv=TRUE)
## Warning in smooth.spline(Prestige$income, Prestige$prestige, cv = TRUE):
## cross-validation with non-unique 'x' values seems doubtful
fit5 <- data.frame(x=fit5$x,y=fit5$y)
ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=fit5, aes(x=x, y=y), col='yellow')

4. Sploty naturalne

fit6 <- lm(income ~ ns(prestige, df=6), Prestige)
fit7 <- lm(income ~ ns(prestige, df=12), Prestige)

ggplot(Prestige) +
    geom_point(aes(x=prestige,y=income)) +
    ggtitle("Naturalne sploty, 6 vs 12 df") +
    geom_line(aes(x=prestige, y=fitted(fit6)), col='darkblue')+
    geom_line(aes(x=prestige, y=fitted(fit7)), col='pink')

Różowa linia (df = 12) lepiej odwzorowuje szczegóły w danych, ale może wprowadzać większe ryzyko nadmiernego dopasowania do szumu w danych.

Splotu naturalne z przedziałem ufności

ggplot(Prestige) +
    geom_point(aes(x=income,y=prestige)) +
    geom_smooth(aes(x=income,y=prestige), method='gam',
      formula = y ~ s(x,k=12))

Zadania 2 - wypadek motocyklowy

Zbiór danych “mcycle” (z pakietu MASS) zawiera n=133 pary punktów czasowych (w ms) i obserwowanych przyspieszeń głowy (w g), które zostały zarejestrowane w symulowanym wypadku motocyklowym.

Do zbadania zależności między czasem a przyspieszeniem wykorzystaj metody regresji nieparametrycznej.

Najpierw wczytaj dane i zwizualizuj zależność między czasem (X) a przyspieszeniem (Y).

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
data("mcycle")

ggplot(mcycle)+
  geom_point(aes(x=times,y=accel))

1. Kernel smoothing

fit <- locpoly(mcycle$times, mcycle$accel,
               degree=0, bandwidth=1, gridsize= 10000) %>% as.tibble
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
fit2 <- locpoly(mcycle$times, mcycle$accel,         
                degree=0, bandwidth=5) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times, y=accel)) +
  geom_line(data=fit, aes(x=x, y=y), col='orange')+ 
  geom_line(data=fit2, aes(x=x, y=y), col='green')+ 
  ggtitle("Kernel smoothing - szerokość pasma 1 vs 4")

Pomarańczowa linia, przedstawiająca szerokość pasma 4 wydaje się być dobrze dopasowana w charakterystycznych dla danych dołku i picku. Z drugiej strony, po przekroczeniu times=30 wyniki przestają równo się układać i linia regresji jest zbyt szczegółowa, zbyt wrażliwa na pojedyncze obserwacje.

2. Loess - wielomian lokalnie kwadratowy (degree=2)

fit3 <- loess(accel ~ times, span=0.2, degree=2,   family="gaussian", mcycle)
fit4 <- loess(accel ~ times, span=0.75, degree=2,   family="gaussian", mcycle)

ggplot(mcycle) +
  geom_point(aes(x=times, y=accel)) +
  geom_line(aes(x=times, y=fitted(fit3)), col='pink')+ 
  geom_line(aes(x=times, y=fitted(fit4)), col='darkred')+ 
  ggtitle("Loess 0.25 vs 0.75")

Na powyższym wykresie bordowa linia, przedstawiająca szerokość pasma na poziomie 0.75 wydaje się być zbyt wygładzona. Totalnie nie odwzorowuje kształtu układania się zmiennych. Natomiast span=0.75 wydaje się być zbyt wrażliwe na pojedyncze obserwacje.

Loess z przedziałem ufności

ggplot(mcycle) +
    geom_point(aes(x=times, y=accel)) +
    geom_smooth(aes(x=times, y=accel), method='loess', span=0.25)
## `geom_smooth()` using formula = 'y ~ x'

Uwzględniając przedział ufności przy span=0.25 widać, że duża część obserwacji mieści się w nim. Do tej pory wydaje się to być najlepiej dopasowaną metodą.

3. Sploty interpolujące

fit5 <- smooth.spline(mcycle$times, mcycle$accel, cv=TRUE)
## Warning in smooth.spline(mcycle$times, mcycle$accel, cv = TRUE):
## cross-validation with non-unique 'x' values seems doubtful
fit5 <- data.frame(x=fit5$x, y=fit5$y)

ggplot(mcycle) +
  geom_point(aes(x=times, y=accel)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=fit5, aes(x=x, y=y), col='yellow')

4. Sploty naturalne

fit6 <- lm(accel ~ ns(times, df=6), mcycle)
fit7 <- lm(accel ~ ns(times, df=12), mcycle)

ggplot(mcycle) +
    geom_point(aes(x=times, y=accel)) +
    ggtitle("Naturalne sploty, 6 vs 12 df") +
    geom_line(aes(x=times, y=fitted(fit6)), col='darkblue')+ 
    geom_line(aes(x=times, y=fitted(fit7)), col='pink')

Różowa linia (df = 12) lepiej odwzorowuje szczegóły w danych, wyłapując skoki, spadki, ale też wypłaszczenia w układających się na wykresie obserwacjach.

Sploty naturalne z przedziałem ufności

ggplot(mcycle) +
    geom_point(aes(x=times, y=accel)) +
    geom_smooth(aes(x=times, y=accel), method='gam', formula = y ~ s(x, k=12))

Wykres dobrze pokazuje zmiany w zależnościach między zmiennymi oraz poziom niepewności w danych. Natomiast obszary o szerokim przedziale ufności (np. na krańcach wykresu) sugerują potrzebę ostrożności przy interpretacji wyników w tych zakresach.