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

Dla zawodów z zarobkami poniżej 10 tys. dolarów prestiż rośnie wyraźnie wraz z dochodem. Natomiast w przedziale od 10 do 25 tysięcy ta zależność jest już dużo słabsza i mniej oczywista. Efekt zaczyna wygasać.

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 = 'red') +
  geom_line(data = fit2, aes(x = x, y = y), col = 'blue') +
  ggtitle("Kernel smoothing")

Czerwona linia (bandwidth = 1) nadmiernie się dopasowuje do szumów. Niebieska linia bandwidth = 5 wygładza dane i lepiej ukazuje ogólny trend.

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

loess1 <- loess(income ~ prestige, span=0.25, degree=2,
  family="gaussian", Prestige)
loess2 <- 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(loess1)), col='pink')+
  geom_line(aes(x=prestige, y=fitted(loess2)), col='grey')+
  ggtitle("Loess 0.25 vs 0.75")

Szara linia (spann = 0,75) wydaje się bardziej redukować szumy oraz wygładzać trend. Linia różowa (span = 0,25) posiada wiele nieregularnych odchyleń dopsaowując się do lokalnych wahań. ### 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.75)
## `geom_smooth()` using formula = 'y ~ x'

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'

Więcej punktów znajduje się w przedziale ufności przy span 0,25 ze względu na większą reakcję na wahania lini trendu. ## 3. Sploty interpolujące

interpolacja <- smooth.spline(Prestige$income, Prestige$prestige, cv=TRUE)
## Warning in smooth.spline(Prestige$income, Prestige$prestige, cv = TRUE):
## krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
interpolacja <- data.frame(x=interpolacja$x,y=interpolacja$y)
ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=interpolacja, aes(x=x, y=y), col='yellow')

4. Sploty naturalne

snatural1 <- lm(income ~ ns(prestige, df=6), Prestige)
snatural2 <- 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(snatural1)), col='darkblue')+
    geom_line(aes(x=prestige, y=fitted(snatural2)), col='pink')

W tym przypadku niebieska linia z tych samych powodów co wcześniej szara jest lepsza. Różowa zbyt dopasowana.

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)
## 
## Dołączanie pakietu: 'MASS'
## Następujący obiekt został zakryty z 'package:dplyr':
## 
##     select
data("mcycle")

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

Wyraźny spadek prędkości w ok. 20 sekundzie mogący świadczyć o momencie zderzenia.

1. Kernel smoothing

fit3 <- 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.
fit4 <- locpoly(mcycle$times, mcycle$accel,         
                degree=0, bandwidth=5) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times, y=accel)) +
  geom_line(data=fit3, aes(x=x, y=y), col='orange')+ 
  geom_line(data=fit4, aes(x=x, y=y), col='green')+ 
  ggtitle("Kernel smoothing - motocykl")

Pomarańczowa linia, przedstawiająca szerokość pasma 4 zbyt dopasowuje się do pojedynczych wahań. Zielona zdecydowanie lepiej ukazuje ogólny trend, ale nie reprezentuje ogromego spadku prędkości w 20 sekundzie.

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

Loes1moto <- loess(accel ~ times, span=0.2, degree=2,   family="gaussian", mcycle)
Loes2moto <- 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(Loes1moto)), col='red')+ 
  geom_line(aes(x=times, y=fitted(Loes2moto)), col='blue')+ 
  ggtitle("Loess MOTO 0.25 vs 0.75")

Podobna sytaucja jak przy poprzednim wykresie. Niebieska linia (loess 0,75) przedstawia o wiele bardziej wygładzone dane jednak nie odwzrorowuje odchylenia przy spadku prędkości. Jednak w tym przypadku dodatkowo nie odpowiednio wizualizuje trend przy końcowej fazie między 45 - 60 sek. Tym razem czerwona linia (loess 0,25) wydaje się być lepiej dopasowana.

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ą.

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

3. Sploty interpolujące

interpolacjamoto <- smooth.spline(mcycle$times, mcycle$accel, cv=TRUE)
## Warning in smooth.spline(mcycle$times, mcycle$accel, cv = TRUE): krzyżowa
## walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
interpolacjamoto<- data.frame(x=interpolacjamoto$x, y=interpolacjamoto$y)

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

4. Sploty naturalne

snaturalmoto<- lm(accel ~ ns(times, df=6), mcycle)
snaturalmoto2 <- 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(snaturalmoto)), col='brown')+ 
    geom_line(aes(x=times, y=fitted(snaturalmoto2)), col='darkgreen')

Zielona linia df=12 znajduje się bliżej punktów w lepszy sposób dopasowując się do danych.

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 pokazuje, jak zmieniają się zależności między zmiennymi i poziom niepewności w danych. Szerokie przedziały ufności na krańcach sugerują, że w tych miejscach trzeba uważać przy interpretacji wyników.