Zadanie 1: Prestiż vs dochód

Zaczynam od zaciągnięcia danych i przedstawienia zależności między prestiżem a dochodem.

data_prestige <- carData::Prestige
library(ggplot2)
## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3
ggplot(data=data_prestige, aes(x=income, y=prestige)) + 
         geom_point() + 
         ggtitle("Presige vs income")

Przedstawiona relacja wskazuje na związek nieliniowy, bardziej przypomina wykres logarytmu. Należy zauważyć, że na pewnym odcinku, to jest do 10000 dolarów dochodu, dane charakteryzuje silna pozytywna liniowa zależność z prestiżem. Po przekroczeniu tej granicy, liniowa prawidłowość nie ma miejsca.

#install.packages('KernSmooth')
library(KernSmooth)
## Warning: pakiet 'KernSmooth' został zbudowany w wersji R 4.3.3
## KernSmooth 2.23 załadowane
## Prawa autorskie M. P. Wand 1997-2009
x <- data_prestige$income
y <- data_prestige$prestige

kernel_smoothing_linear <- locpoly(x, y, degree=1, bandwidth = 2000) 
kernel_smoothing_squared <- locpoly(x, y, degree=2, bandwidth = 2000)

data_linear <- data.frame(x=kernel_smoothing_linear$x, y=kernel_smoothing_linear$y)
data_squared <- data.frame(x=kernel_smoothing_squared$x, y=kernel_smoothing_squared$y)

colors_locpoly <- c('Linear' = 'green', 'Squared'='orange')

ggplot(data_prestige) + 
  geom_point(aes(x=income, y=prestige)) + 
  ggtitle("Presige vs income local polynomial") + 
  geom_line(data=data_linear, aes(x=x, y=y, color='Linear')) + 
  geom_line(data=data_squared, aes(x=x, y=y, color='Squared')) + 
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_locpoly)

Na początku użyta została regresja nieparametryczna w postaci localnego wielomianu, a konkretnie dwóch, odpowiednio pierwszego i drugiego stopnia. Z powyższego wykresu można wywnioskować, że wielomian kwadratowy lepiej radzi sobie z dopasowaniem do danych, oczywiście w zdroworozsądkowych ramach szerokości pasma.

smr_025 <- loess(y ~ x, span=0.25,  degree=2, family='gaussian', data_prestige)
smr_075 <- loess(y ~ x, span=0.75,  degree=2, family='gaussian', data_prestige)
smr_15 <- loess(y ~ x, span=1.5,  degree=2, family='gaussian', data_prestige)
colors_loess <- c('span_025'='blue', 'span_075'='orange', 'span_15'='red')

ggplot(data_prestige) +
  geom_point(aes(x=income,y=prestige)) +
  ggtitle("Presige vs income (Loess)") +
  geom_line(aes(x=x, y=fitted(smr_025), color='span_025')) + 
  geom_line(aes(x=x, y=fitted(smr_075), color='span_075')) + 
  geom_line(aes(x=x, y=fitted(smr_15), color='span_15')) + 
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_loess)

Przedstawiona powyżej została implementacja funkcji loess (lokalnie kwadratowa), czyli odpowiednika wcześniejszego wielomianu drugiego stopnia. Co jednak ciekawe, przedstawione modelowanie loess opiera się na zmianie argumentu span, który pozwala użytkownikowi kontrolować szerokość pasma, określoną w tej funkcji jako procent pokrytych danych. Zgodnie z teorią, niższa wartość argumentu span powoduje lepsze dopasowanie, zaś wyższa pozwala uogólnić, czyli wygładzić.

Powyższy wykres można również osadzić w ramach tzw. przedziału ufności. W tym celu wybrany został tylko jeden wielomian, aby pokazać ideę.

ggplot(data_prestige) + 
  geom_point(aes(x=income,y=prestige)) + 
  geom_smooth(aes(x=income,y=prestige, color='span_025'), method='loess', span=0.25) +
  ggtitle("Presige vs income (Loess) with confidence interval") + 
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_loess)
## `geom_smooth()` using formula = 'y ~ x'

smr_025 <- smooth.spline(
  data_prestige$income, 
  data_prestige$prestige, 
  spar=0.25,
  cv=TRUE)
## Warning in smooth.spline(data_prestige$income, data_prestige$prestige, spar =
## 0.25, : krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_075 <- smooth.spline(
  data_prestige$income, 
  data_prestige$prestige, 
  spar=0.75,
  cv=TRUE)
## Warning in smooth.spline(data_prestige$income, data_prestige$prestige, spar =
## 0.75, : krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_15 <- smooth.spline(
  data_prestige$income, 
  data_prestige$prestige, 
  spar=1.5,
  cv=TRUE)
## Warning in smooth.spline(data_prestige$income, data_prestige$prestige, spar =
## 1.5, : krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_025_df <- data.frame(x=smr_025$x, y=smr_025$y)
smr_075_df <- data.frame(x=smr_075$x, y=smr_075$y)
smr_15_df <- data.frame(x=smr_15$x, y=smr_15$y)
ggplot(data_prestige) + 
  geom_point(aes(x=income, y=prestige)) + 
  ggtitle("Income vs Prestige (Sploty interpolujące, lambda wybrana przez CV)") + 
  geom_line(data=smr_025_df, aes(x=x, y=y, color='span_025')) +
  geom_line(data=smr_075_df, aes(x=x, y=y, color='span_075')) +
  geom_line(data=smr_15_df, aes(x=x, y=y, color='span_15')) +
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_loess)

Przedstawione powyżej sploty interpolujące, w którym lambda została wybrana metodą walidacji krzyżowej, pokazały bardzo ciekawe wyniki. Okazało się, ze im wyższa wartość szerokości pasma, tym interpolacja przebiega liniowo i prawie w ogóle nie łączy obserwacji.

library(splines)
colors_splots <- c('ns_3'='blue', 'ns_10'='red')

natural_splots_3 <- lm(prestige ~ ns(income, df=3), data_prestige)
natural_splots_10 <- lm(prestige ~ ns(income, df=10), data_prestige)

ggplot(data_prestige) +
  geom_point(aes(x=income,y=prestige)) +
  ggtitle("Income vs Prestige (Natural splots)") +
  geom_line(aes(x=x, y=fitted(natural_splots_3), color='ns_3')) + 
  geom_line(aes(x=x, y=fitted(natural_splots_10), color='ns_10')) +
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_splots)

ggplot(data_prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_smooth(aes(x=income,y=prestige, color='ns_3'), method='gam',
              formula = y ~ s(x,k=3)) + 
  geom_smooth(aes(x=income,y=prestige, color='ns_10'), method='gam',
              formula = y ~ s(x,k=10)) +
  ggtitle('Prestige vs Income (Natural splots with confidence interval)') + 
  labs(x='Income',
       y='Prestige',
       color='Legend') + 
  scale_color_manual(values=colors_splots)

Zadanie 2: Wypadek motocyklowy

Zaczynam od zaciągnięcia danych i przedstawienia zależności między czasem a przyspieszeniem.

data_moto <- MASS::mcycle
ggplot(data=data_moto, aes(x=times, y=accel)) + 
         geom_point() + 
         ggtitle("Acceleration of head vs times") + 
         labs(x='Times (ms)', 
              y='Acceleration of head (g)')

Przedstawiona relacja wskazuje na wyraźny związek nieliniowy, gdzie można wyróżnić odcinki, które przypominają proste fuynkcje. Warto wskazać, że do 15 ms, tak naprawdę zmiana czasu nie niesie zmian w przyspieszeniu. Po około 15 ms przyspieszenie rośnie wykładniczo w dół, osiągając dołek w okolicach 20 ms, a następnie silnie rośnie do ponad 50 g w czasie 20-30 ms, następnie spada w rejony 0 do 40 ms i od tego momentu wartości się stabilizują.

x_moto <- data_moto$times
y_moto <- data_moto$accel

kernel_smoothing_linear_moto <- locpoly(x_moto, y_moto, degree=1, bandwidth = 5) 
kernel_smoothing_squared_moto <- locpoly(x_moto, y_moto, degree=2, bandwidth = 5) 

data_linear_moto <- data.frame(x=kernel_smoothing_linear_moto$x, y=kernel_smoothing_linear_moto$y)
data_squared_moto <- data.frame(x=kernel_smoothing_squared_moto$x, y=kernel_smoothing_squared_moto$y)

ggplot(data_moto) + 
  geom_point(aes(x=times, y=accel)) + 
  ggtitle("Acceleration vs times local polynomial") + 
  geom_line(data=data_linear_moto, aes(x=x, y=y, color='Linear')) + 
  geom_line(data=data_squared_moto, aes(x=x, y=y, color='Squared')) + 
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_locpoly)

Na początku użyta została regresja nieparametryczna w postaci localnego wielomianu, a konkretnie dwóch, odpowiednio pierwszego i drugiego stopnia. Z powyższego wykresu można wywnioskować, że wielomian kwadratowy lepiej radzi sobie z dopasowaniem do danych, oczywiście w zdroworozsądkowych ramach szerokości pasma :).

smr_033_moto <- loess(y_moto ~ x_moto, span=0.33,  degree=2, family='gaussian', data_moto)
smr_066_moto <- loess(y_moto ~ x_moto, span=0.66,  degree=2, family='gaussian', data_moto)
smr_1_moto  <- loess(y_moto ~ x_moto, span=1,  degree=2, family='gaussian', data_moto)
colors_loess_moto <- c('span_033'='green', 'span_066'='purple', 'span_1'='brown')

ggplot(data_moto) +
  geom_point(aes(x=times, y=accel)) +
  ggtitle("Acceleration vs times local polynomial") +
  geom_line(aes(x=x_moto, y=fitted(smr_033_moto), color='span_033')) + 
  geom_line(aes(x=x_moto, y=fitted(smr_066_moto), color='span_066')) + 
  geom_line(aes(x=x_moto, y=fitted(smr_1_moto), color='span_1')) + 
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_loess_moto)

Przedstawiona powyżej została implementacja funkcji loess (lokalnie kwadratowa), czyli odpowiednika wcześniejszego wielomianu drugiego stopnia. Co jednak ciekawe, przedstawione modelowanie loess opiera się na zmianie argumentu span, który pozwala użytkownikowi kontrolować szerokość pasma, określoną w tej funkcji jako procent pokrytych danych. Zgodnie z teorią, niższa wartość argumentu span powoduje lepsze dopasowanie, zaś wyższa pozwala uogólnić, czyli wygładzić.

Powyższy wykres można również osadzić w ramach tzw. przedziału ufności. W tym celu wybrany został tylko jeden wielomian, aby pokazać ideę.

ggplot(data_moto) + 
  geom_point(aes(x=times, y=accel)) + 
  geom_smooth(aes(x=times, y=accel, color='span_033'), method='loess', span=0.33) +
  ggtitle("Acceleration vs times (Loess) with confidence interval") + 
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_loess_moto)
## `geom_smooth()` using formula = 'y ~ x'

smr_033_moto <- smooth.spline(
  data_moto$times, 
  data_moto$accel, 
  spar=0.33,
  cv=TRUE)
## Warning in smooth.spline(data_moto$times, data_moto$accel, spar = 0.33, :
## krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_066_moto <- smooth.spline(
  data_moto$times, 
  data_moto$accel, 
  spar=0.66,
  cv=TRUE)
## Warning in smooth.spline(data_moto$times, data_moto$accel, spar = 0.66, :
## krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_1_moto <- smooth.spline(
  data_moto$times, 
  data_moto$accel, 
  spar=1,
  cv=TRUE)
## Warning in smooth.spline(data_moto$times, data_moto$accel, spar = 1, cv =
## TRUE): krzyżowa walidacja z nieunikalnymi wartościami 'x' wydaje się wątpliwa
smr_033_df_moto <- data.frame(x=smr_033_moto$x, y=smr_033_moto$y)
smr_066_df_moto <- data.frame(x=smr_066_moto$x, y=smr_066_moto$y)
smr_1_df_moto   <- data.frame(x=smr_1_moto$x, y=smr_1_moto$y)
ggplot(data_moto) + 
  geom_point(aes(x=times, y=accel)) + 
  ggtitle("Acceleration vs times (Sploty interpolujące, lambda wybrana przez CV)") + 
  geom_line(data=smr_033_df_moto, aes(x=x, y=y, color='span_033')) +
  geom_line(data=smr_066_df_moto, aes(x=x, y=y, color='span_066')) +
  geom_line(data=smr_1_df_moto, aes(x=x, y=y, color='span_1')) +
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_loess_moto)

Przedstawione powyżej sploty interpolujące, w którym lambda została wybrana metodą walidacji krzyżowej, pokazały bardzo ciekawe wyniki. Okazało się, ze im wyższa wartość szerokości pasma, tym interpolacja przebiega jako wielomian drugiego stopnia i prawie w ogóle nie łączy obserwacji.

#library(splines)
colors_splots_moto <- c('ns_2'='green', 'ns_8'='black')

natural_splots_2 <- lm(accel ~ ns(times, df=2), data_moto)
natural_splots_8 <- lm(accel ~ ns(times, df=8), data_moto)

ggplot(data_moto) +
  geom_point(aes(x=times,y=accel)) +
  ggtitle("Acceleration vs times (Natural splots)") +
  geom_line(aes(x=x_moto, y=fitted(natural_splots_2), color='ns_2')) + 
  geom_line(aes(x=x_moto, y=fitted(natural_splots_8), color='ns_8')) +
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_splots_moto)

ggplot(data_moto) +
  geom_point(aes(x=times,y=accel)) +
  geom_smooth(aes(x=times,y=accel, color='ns_2'), method='gam',
              formula = y ~ s(x,k=2)) + 
  geom_smooth(aes(x=times,y=accel, color='ns_8'), method='gam',
              formula = y ~ s(x,k=8)) +
  ggtitle('Acceleration vs times (Natural splots with confidence interval)') + 
  labs(x='Times (ms)',
       y='Acceleration (g)',
       color='Legend') + 
  scale_color_manual(values=colors_splots_moto)
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): wymiar podstawy, k, zwiększył się do minimalnego możliwego