Raport 3

Regresja nieparametryczna

1 Zadanie 1 - dane Prestige

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")
attach(Prestige)

## Następujący obiekt został zakryty z package:datasets:
## 
##     women

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

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

Metody regresji nieparametrycznej.

1.1 Local polynomial regression (locpoly) - metoda lokalnego wygładzania

Główne parametry, takie jak degree i bandwidth, mają kluczowy wpływ na sposób, w jaki wygładzana jest zależność.

Parametr degree określa stopień wielomianu użytego do lokalnej aproksymacji. Degree=0 oznacza użycie lokalnych średnich (co odpowiada regresji najbliższych sąsiadów).

Degree=1 oznacza wykorzystanie lokalnej regresji liniowej

Degree=2 oznacza wykorzystanie lokalnej regresji kwadratowej.

Parametr bandwidth decyduje o rozpiętości jądra, czyli definiuje rozmiar obszaru w jakim obserwacje wpływają na funkcję gęstości.Im większa wartość bandwidth, tym bardziej wygładzona jest funkcja. Zarówno zbyt niska wartość bandwidth może być nieodpowiednia (overfitting- model zbyt szczegółowy), jak i również zbyt wysoka wartość tego parametru (underfitting- zbyt duże wygładzenie, zbyt duże uogólnienie).

1.1.1 Porównanie różnych bandwidth

fit1 <- locpoly(prestige, income,
        degree=0, bandwidth=1) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue')+
  ggtitle("Kernel smoothing - szerokość pasma równa 1")

Zwiększam szerokość pasma do 5

fit2 <- locpoly(prestige, income,
        degree=0, bandwidth=5) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit2, aes(x=x,y=y), col='red')+
  ggtitle("Kernel smoothing - szerokość pasma równa 5")

Porównanie:

Niebieska linia: szerokość pasma 1 Czerwona linia: szerokość pasma 5

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue')+
  geom_line(data=fit2, aes(x=x,y=y), col='red')+
  ggtitle("Kernel smoothing - porówananie dwóch szerokości pasm")

Aby dopasować szerokośc pasma - cross validation

dopasowana <- dpill(prestige, income)
fit3 <- locpoly(prestige, income,
        degree=0, bandwidth=dopasowana) %>% as.tibble

paste("Dopasowana szerokość pasma wynosi:", round(dopasowana, 4))

## [1] "Dopasowana szerokość pasma wynosi: 1.6936"

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income))+
  geom_line(data=fit3, aes(x=x,y=y), col='green',size=1)

1.1.2 Porównanie różnych degree (stopień wielomianu)

fit1 <- locpoly(prestige, income,
        degree=0, bandwidth=dopasowana) %>% as.tibble
fit2 <- locpoly(prestige, income,
        degree=1, bandwidth=dopasowana) %>% as.tibble
fit3 <- locpoly(prestige, income,
        degree=2, bandwidth=dopasowana) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue',size=0.8)+
  geom_line(data=fit2, aes(x=x,y=y), col='red',size=0.8)+
  geom_line(data=fit3, aes(x=x,y=y), col='green',size=0.8)+
  ggtitle("Kernel smoothing - porówananie trzech stopni wielomianu")

Można zauważyć, że wybór degree wpływa na funkcję gęstości na wykresie. W przypadku degree=2 występuje overfitting. Lepszym wyborem byłoby degree równe 0 lub 1.

1.2 Loess - lokalnie kwadratowy

smr <- loess (income~prestige,span=0.75, degree=2,
               family="gaussian")

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income))+
  geom_line(aes(x=prestige,y=fitted(smr)), col='skyblue',size=1)

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income))+
  geom_smooth(aes(x=prestige,y=income),method='loess',span=0.22)

## `geom_smooth()` using formula = 'y ~ x'

1.3 Sploty interpolujące - smooth.spline

1.3.1 Zmienianie parametru ‘spar’

fit1 <-smooth.spline(prestige, income, spar=0.2) 
smr1 <- data.frame(x=fit1$x,y=fit1$y)

fit2 <-smooth.spline(prestige, income, spar=0.8) 
smr2 <- data.frame(x=fit2$x,y=fit2$y)

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=smr1, aes(x=x,y=y), col='blue',size=0.8)+
  geom_line(data=smr2, aes(x=x,y=y), col='red',size=0.8)+
  ggtitle("Sploty interpolujące - porównanie różnych wartości parametru 'spar'")

Spar=0.2 powoduje zbyt mocne dopasowanie do danych. Spar=0.8 prezentuje zbyt duże uogólnienie. Wniosek: Zmniejszanie ‘spar’ powoduje większe dopasowanie, a zwiększanie - większe uogólnienie.

Aby automatycznie wybrać szerokośc pasma można zastosować cross validation

1.3.2 Zastosowanie cross validation (CV=TRUE)

smr2 <- smooth.spline(
  prestige,
  income,
  cv=TRUE)
smr2 <- data.frame(x=smr2$x,y=smr2$y)
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=smr2, aes(x=x, y=y), col='blue')

1.4 Porównanie funkcji locpoly i smooth.spline

Niebieska linia - locpoly Czerwona linia - smooth.spline

fit <- locpoly(prestige, income,
        degree=0, bandwidth=dopasowana) %>% as.tibble

smr <- smooth.spline(
  prestige,
  income,
  cv=TRUE)
smr <- data.frame(x=smr$x,y=smr$y)

ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit, aes(x=x,y=y), col='blue')+
  geom_line(data=smr, aes(x=x,y=y), col='red')+
  ggtitle("Porównanie funkcji locpoly i smooth.spline")

1.5 Naturalne sploty

library(splines)

fit <-lm(income~ns(prestige,df=6),Prestige)
fit2 <-lm(income~ns(prestige,df=12),Prestige)

ggplot(Prestige)+
  geom_point(aes(x=prestige,y=income))+
  ggtitle("Naturalne sploty, 6 df(niebieski) i 12 df(czerwony)")+
  geom_line(aes(x=prestige,y=fitted(fit)),col='blue')+
  geom_line(aes(x=prestige,y=fitted(fit2)),col='red')

1.6 Sploty oraz geom_smooth()

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

fit <- lm(income ~ ns(prestige, df = 6), Prestige)
pred <- data.frame(prestige = Prestige$prestige,
                   predict(fit, interval = "confidence"))
ggplot(Prestige, aes(x = prestige, y = income)) +
    geom_point() +
    geom_line(data = pred, aes(y = fit), color = "blue") +
    geom_ribbon(data = pred, aes(ymin = lwr, ymax = upr), alpha = 0.2)

2 Zadanie 2 - dane mcycle

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")
attach(mcycle)

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

Zależność wygląda na nieliniową.

Przyspieszenie jest stabilne od 0-15 ms, spada od ok. 15-20 ms, rośnie od 20-30 ms, spada od 30-40 ms, a następnie zaczyna się stabilizować.

2.1 Local polynomial regression (locpoly) - metoda lokalnego wygładzania

2.1.1 Porównanie różnych bandwidth

fit1 <- locpoly(times, accel,
        degree=0, bandwidth=1) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue')+
  ggtitle("Kernel smoothing - szerokość pasma równa 1")

Zwiększam szerokość pasma do 5

fit2 <- locpoly(times, accel,
        degree=0, bandwidth=5) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit2, aes(x=x,y=y), col='red')+
  ggtitle("Kernel smoothing - szerokość pasma równa 5")

Porównanie:

Niebieska linia: szerokość pasma 1 Czerwona linia: szerokość pasma 5

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue')+
  geom_line(data=fit2, aes(x=x,y=y), col='red')+
  ggtitle("Kernel smoothing - porówananie dwóch szerokości pasm")

Aby dopasować szerokośc pasma - cross validation

dopasowana <- dpill(times, accel)
fit3 <- locpoly(times, accel,
        degree=0, bandwidth=dopasowana) %>% as.tibble

paste("Dopasowana szerokość pasma wynosi:", round(dopasowana, 4))

## [1] "Dopasowana szerokość pasma wynosi: 1.4453"

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel))+
  geom_line(data=fit3, aes(x=x,y=y), col='green',size=1)

2.1.2 Porównanie różnych degree (stopień wielomianu)

fit1 <- locpoly(times, accel,
        degree=0, bandwidth=dopasowana) %>% as.tibble
fit2 <- locpoly(times, accel,,
        degree=1, bandwidth=dopasowana) %>% as.tibble
fit3 <- locpoly(times, accel,
        degree=2, bandwidth=dopasowana) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit1, aes(x=x,y=y), col='blue',size=0.8)+
  geom_line(data=fit2, aes(x=x,y=y), col='red',size=0.8)+
  geom_line(data=fit3, aes(x=x,y=y), col='green',size=0.8)+
  ggtitle("Kernel smoothing - porówananie trzech stopni wielomianu")

2.2 Loess - lokalnie kwadratowy

smr <- loess (accel~times,span=0.75, degree=2,
               family="gaussian")

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel))+
  geom_line(aes(x=times,y=fitted(smr)), col='skyblue',size=1)

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel))+
  geom_smooth(aes(x=times,y=accel),method='loess',span=0.22)

## `geom_smooth()` using formula = 'y ~ x'

2.3 Sploty interpolujące - smooth.spline

2.3.1 Zmienianie parametru ‘spar’

fit1 <-smooth.spline(times, accel, spar=0.2) 
smr1 <- data.frame(x=fit1$x,y=fit1$y)

fit2 <-smooth.spline(times, accel,, spar=0.8) 
smr2 <- data.frame(x=fit2$x,y=fit2$y)

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=smr1, aes(x=x,y=y), col='blue',size=0.8)+
  geom_line(data=smr2, aes(x=x,y=y), col='red',size=0.8)+
  ggtitle("Sploty interpolujące - porównanie różnych wartości parametru 'spar'")

Spar=0.2 powoduje zbyt mocne dopasowanie do danych. Spar=0.8 prezentuje zbyt duże uogólnienie. Wniosek: Zmniejszanie ‘spar’ powoduje większe dopasowanie, a zwiększanie - większe uogólnienie.

Aby automatycznie wybrać szerokośc pasma można zastosować cross validation

2.3.2 Zastosowanie cross validation (CV=TRUE)

smr2 <- smooth.spline(
  times,
  accel,
  cv=TRUE)
smr2 <- data.frame(x=smr2$x,y=smr2$y)
ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=smr2, aes(x=x, y=y), col='blue')

2.4 Porównanie funkcji locpoly i smooth.spline

fit <- locpoly(times, accel,
        degree=0, bandwidth=dopasowana) %>% as.tibble

smr <- smooth.spline(
  times,
  accel,
  cv=TRUE)
smr <- data.frame(x=smr$x,y=smr$y)

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit, aes(x=x,y=y), col='blue')+
  geom_line(data=smr, aes(x=x,y=y), col='red')+
  ggtitle("Porównanie funkcji locpoly i smooth.spline")

2.5 Naturalne sploty

library(splines)

fit <-lm(accel~ns(times,df=6),mcycle)
fit2 <-lm(accel~ns(times,df=12),mcycle)

ggplot(mcycle)+
  geom_point(aes(x=times,y=accel))+
  ggtitle("Naturalne sploty, 6 df(niebieski) i 12 df(czerwony)")+
  geom_line(aes(x=times,y=fitted(fit)),col='blue')+
  geom_line(aes(x=times,y=fitted(fit2)),col='red')

2.6 Sploty oraz geom_smooth()

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

fit <- lm(accel ~ ns(times, df = 6), mcycle)
pred <- data.frame(times = mcycle$times,
                   predict(fit, interval = "confidence"))
ggplot(mcycle, aes(x = times, y = accel)) +
    geom_point() +
    geom_line(data = pred, aes(y = fit), color = "blue") +
    geom_ribbon(data = pred, aes(ymin = lwr, ymax = upr), alpha = 0.2)