Raport 3

Regresja nieparametryczna

library(car)
library(KernSmooth)
library(tidyverse)
library(ggplot2)
library(splines)
library(MASS)

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

Prestige$income <- as.numeric(Prestige$income)
Prestige$prestige <- as.numeric(Prestige$prestige)

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

Patrząc na wykres powiązania pomiędzy income a prestige widzim że kropki układają w sposób sugerujący że związek tych czynników jest nieliniowy. Dla zawodów, które zarabiają mniej niż $10K,widzimy , że wartości zależności pomiędzy dochodem a prestiżem układają się nieliniowo. W przypadku zawodów, które zarabiają od 10 do 25 tysięcy dolarów, związek ten ma znacznie inne nachylenie.

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

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. Parametr bandwidth określa szerokość pasma, które jest używane do wygładzania.

Degree=0 ->lokalna średnich (regresja najbliższych sąsiadów)

Degree=1 -> lokalnej regresja liniowa

Degree=2 ->lokalna regresja kwadratowa

1.1.1 Porównanie różnych bandwidth

fit1a <- locpoly(income,prestige, 
        degree=0, bandwidth=800,gridsize = 100000) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_line(data=fit1a, aes(x=x,y=y), col='blue')+
  ggtitle("Szerokość pasma równa 800")+
  labs(subtitle = "Kernel Smothing")

Zwiększam szerokość pasma do 5

fit1b <- locpoly(income,prestige,
        degree=0, bandwidth=1000,gridsize = 260000) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_line(data=fit1b, aes(x=x,y=y), col='red')+
  ggtitle("Szerokość pasma równa 1000")+
  labs(subtitle = "Kernel Smothing")

Porównanie:

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

ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_line(data=fit1a, aes(x=x,y=y), col='blue')+
  geom_line(data=fit1b, aes(x=x,y=y), col='red')+
  ggtitle("Porówananie dwóch szerokości pasm")+
  labs(subtitle = "Kernel Smothing")

Aby dopasować szerokośc pasma - cross validation

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

ggplot(Prestige) +
  geom_point(aes(x=income, y=prestige))+
  geom_line(data=fit1d, aes(x=x,y=y), col='lightgreen',size=1)+
ggtitle("Dopasowana szerokość pasma wynosi:", round(dopasowana, 4))

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

fit2a <- locpoly(income,prestige,
        degree=0, bandwidth=dopasowana) %>% as.tibble
fit2b <- locpoly(income,prestige,
        degree=1, bandwidth=dopasowana) %>% as.tibble
fit2c <- locpoly(income,prestige,
        degree=2, bandwidth=dopasowana) %>% as.tibble

ggplot(Prestige) +
  geom_point(aes(income,prestige)) +
  geom_line(data=fit2a, aes(x=x,y=y), col='darkblue',size=0.8)+
  geom_line(data=fit2b, aes(x=x,y=y), col='hotpink',size=0.8)+
  geom_line(data=fit2c, aes(x=x,y=y), col='darkgoldenrod2',size=0.8)+
  ggtitle("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=0 występuje underfitting. Lepszym wyborem byłoby degree równe 1 lub 2.

1.2 Lokalnie kwadratowy-Loess

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

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

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

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

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

1.3 Sploty interpolujące - smooth.spline

1.3.1 Zmienianie parametru ‘spar’

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

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

ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_line(data=smr3a, aes(x=x,y=y), col='lawngreen',size=0.8)+
  geom_line(data=smr3b, aes(x=x,y=y), col='plum3',size=0.8)+
  ggtitle("Porównanie różnych wartości parametru 'spar'")+
  labs(subtitle = "Sploty interpolujące")

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

1.3.2 Zastosowanie cross validation (CV=TRUE)

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


ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV:", round(lambda_value, 4)) +
  geom_line(data=smr3d, aes(x=x, y=y), col='darkgoldenrod2')

1.4 Porównanie funkcji locpoly i smooth.spline

Fioletowa linia - locpoly Różowa linia - smooth.spline

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

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

ggplot(Prestige) +
  geom_point(aes(x=income,y=prestige)) +
  geom_line(data=FIT, aes(x=x,y=y), col='blueviolet')+
  geom_line(data=SMR, aes(x=x,y=y), col='hotpink')+
ggtitle("Porównanie funkcji locpoly i smooth.spline") + 
scale_color_manual(values = c("locpoly" = "blueviolet", "smooth.spline" = "hotpink"), 
                   labels = c("locpoly", "smooth.spline")) + labs(color = "Function")

1.5 Naturalne sploty

Naturalne sploty, 6 df oznaczone kolorem niebieskim
10 df oznaczone kolorem fioletowym 12 df oznaczone kolorem zielonym

Fit1 <-lm(prestige~ns(income,df=6),Prestige)
Fit2 <-lm(prestige~ns(income,df=10),Prestige)
Fit3 <-lm(prestige~ns(income,df=12),Prestige)

ggplot(Prestige)+
  geom_point(aes(x=income,y=prestige))+
  ggtitle("Naturalne sploty")+
  geom_line(aes(x=income,y=fitted(Fit1)),col='darkblue')+
  geom_line(aes(x=income,y=fitted(Fit2)),col='darkmagenta')+
  geom_line(aes(x=income,y=fitted(Fit3)),col='darkgreen')

1.6 Sploty oraz geom_smooth()

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

Predykcjia z przedziałem ufności 0,1

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

2 Zadanie 2

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

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

Przyspieszenie jest

-stabilne od 0-15 ms

-spada gwałtownie od ok. 15-20 ms

-rośnie szybko od 20-30 ms

-spada ponownie od 30-40 ms,

-zaczyna się stabilizować

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

2.1.1 Porównanie różnych bandwidth

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

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit_1-1, aes(x=x,y=y), col='hotpink')+
  ggtitle("Szerokość pasma równa 1")+
  labs(subtitle = "Kernel Smothing")

Zwiększam szerokość pasma do 5

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

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit_2, aes(x=x,y=y), col='darkblue')+
  ggtitle("Szerokość pasma równa 5")+
  labs(subtitle = "Kernel Smothing")

Porównanie:

Różowa linia: szerokość pasma 1

Niebieska linia: szerokość pasma 5

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit_1, aes(x=x,y=y), col='hotpink')+
  geom_line(data=fit_2, aes(x=x,y=y), col='darkblue')+
  ggtitle("Porówananie dwóch szerokości pasm")+
  labs(subtitle = "Kernel smoothing")

2.1.2 Cross validation

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

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

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

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

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

fit_1a<- locpoly(times, accel,
        degree=0, bandwidth=dopasowana_2) %>% as.tibble
fit_2a <- locpoly(times, accel,,
        degree=1, bandwidth=dopasowana_2) %>% as.tibble
fit_3a <- locpoly(times, accel,
        degree=2, bandwidth=dopasowana_2) %>% as.tibble

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit_1a, aes(x=x,y=y), col='darkblue',size=0.8)+
  geom_line(data=fit_2a, aes(x=x,y=y), col='hotpink',size=0.8)+
  geom_line(data=fit_3a, aes(x=x,y=y), col='darkmagenta',size=0.8)+
  ggtitle("Porówananie trzech stopni wielomianu")+
  labs(subtitle = "Kernel Smothing")

2.2 Loess - lokalnie kwadratowy

smr_1 <- 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_1)), col='darkgoldenrod2',size=1)

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

2.3 Sploty interpolujące - smooth.spline

2.3.1 Zmienianie parametru ‘spar’

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

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

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=smr_4, aes(x=x,y=y), col='blueviolet',size=0.8)+
  geom_line(data=smr_5, aes(x=x,y=y), col='darkgoldenrod2',size=0.8)+
  ggtitle("Porównanie różnych wartości parametru 'spar'")+
  labs(subtitle = "Sploty interpolujące")

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

2.3.2 Zastosowanie cross validation (CV=TRUE)

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

2.4 Porównanie funkcji locpoly i smooth.spline

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

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

ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  geom_line(data=fit_6, aes(x=x,y=y), col='skyblue')+
  geom_line(data=smr_6, aes(x=x,y=y), col='hotpink')+
  ggtitle("Porównanie funkcji locpoly i smooth.spline")

2.5 Naturalne sploty

fit_7 <-lm(accel~ns(times,df=6),mcycle)
fit_7b <-lm(accel~ns(times,df=12),mcycle)
fit_7c <-lm(accel~ns(times,df=12),mcycle)

ggplot(mcycle)+
  geom_point(aes(x=times,y=accel))+
  ggtitle("Naturalne sploty")+
  labs(subtitle = "6 kolor różowy, 10 kolor fioletowy,12 kolor granatowy")+
  geom_line(aes(x=times,y=fitted(fit_7)),col='hotpink')+
  geom_line(aes(x=times,y=fitted(fit_7b)),col='darkmagenta')+
  geom_line(aes(x=times,y=fitted(fit_7c)),col='darkblue')

2.6 Sploty oraz geom_smooth()

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

fit_8 <- lm(accel ~ ns(times, df = 6), data = mcycle)
pred_8 <- data.frame(times = mcycle$times,
                     predict(fit_8, interval = "confidence"))
colnames(pred_8) <- c("times", "fit_8", "lwr", "upr")

ggplot(mcycle, aes(x = times, y = accel)) +
  geom_point() +
  geom_line(data = pred_8, aes(x = times, y = fit_8), color = "blueviolet") +
  geom_ribbon(data = pred_8, aes(x = times, ymin = lwr, ymax = upr), alpha = 0.1)