Nieklasyczne metody statystyki

Regresja nieparametryczna

1 Zad 1

Zestaw danych “Prestige” (z pakietu “car”) zawiera dane nt. prestiżu 102 Kanadyjskich zawodów z 1971 roku. Zawiera on zmienne takie jak:

-education - która określa średnią liczbę lat edukacji wymaganych dla danego zawodu

-income - przedstawiający średni roczny dochód osób pracujących w danym zawodzie

-women - który informuje o procentowej ilości kobiet w każdym zawodzie

-prestige - który odnosi się do społecznie postrzeganego prestiżu zawodu

-census - czyli numer klasyfikacyjny zawodu zgodnie z klasyfikacją spisu powszechnego

-type - który klasyfikuje zawody na kategorie takie jak profesjonalne czy rolnicze

Do zbadania zależności między prestiżem a różnymi zmiennymi wykorzystaną zostane metody regresji nieparametrycznej. Na początku wczytujemy dane.

data("Prestige")
attach(Prestige)

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

1.1 Wykres ukazujący relację między dochodem, a prestiżem.

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. Podsumowując, można zauważyć, że w miarę jak dochód wzrasta, prestiż zawodu nie rośnie proporcjonalnie i pojawia się większe rozproszenie wartości prestiżu.

ggplot(Prestige, aes(x = income, y = prestige)) +
  geom_point(color = "purple", size = 3, alpha = 0.6) +  
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "red") +  
  labs(
    title = "Relacja między dochodem a prestiżem zawodu",
    caption = "Dane źródłowe: Pakiet 'car'"
  ) +
  theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

1.2 Estymacja nieliniowej zależności

Poniżej wykorzystano metodę locpoly() z pakietu KernSmooth do estymacji nieliniowej zależności, gdzie stopień polinomu wynosi zero (średnie jądrowe), a szerokość pasma ustalono na pięć.

fit <- locpoly(prestige, income,
               degree=0, bandwidth=5) %>% as_tibble()
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit, aes(x=x,y=y), col='purple')

Na niższych poziomach prestiżu dochód rośnie powoli, co sugeruje, że dla mniej prestiżowych zawodów wzrost prestiżu nie przekłada się znacząco na dochody. Dla zawodów o średnim poziomie prestiżu (około 25–50 punktów) dochód rośnie bardziej liniowo, co może sugerować umiarkowaną korelację między prestiżem a dochodem w tej grupie zawodów.Dla zawodów o wyższym prestiżu (powyżej 50 punktów) dochód rośnie wykładniczo, co sugeruje, że wśród wysoko prestiżowych zawodów nawet niewielki wzrost prestiżu może wiązać się z bardzo dużym wzrostem dochodów.

1.2.1 Zmiana szerokości pasma

Wartość parametru bandwidth (w tym przypadku 10) kontroluje poziom wygładzenia linii trendu. Mniejsze wartości pasma powodują większą wrażliwość na lokalne wahania danych, co prowadzi do bardziej szczegółowej (czasem zygzakowatej) linii trendu. Z kolei większe wartości bandwidth wygładzają krzywą, co zmniejsza wpływ lokalnych fluktuacji na estymację trendu.

fit <- locpoly(prestige, income,
               degree=0, bandwidth=10) %>% as_tibble()
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit, aes(x=x,y=y), col='purple')

fit <- locpoly(prestige, income,
               degree=0, bandwidth=3) %>% as_tibble()
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit, aes(x=x,y=y), col='purple')

1.2.2 Manipulacja degree wielomianu

Stopień polinomu (degree) wpływa na złożoność modelu - większe wartości pozwalają na bardziej złożone krzywe, które mogą lepiej oddać zawiłości danych.Poniżej znajduje się zmodyfikowany kod, który wykorzystuje degree=2 i szerokość pasma bandwidth=10.

fit <- locpoly(prestige, income,
               degree=2, bandwidth=10) %>% as_tibble()
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  geom_line(data=fit, aes(x=x,y=y), col='purple')

1.3 Interpolacja splotów

Poniżej wykorzytano metodę interpolacji splotów (spline) do wygładzenia zależności między dwiema zmiennymi.Użycie splotów interpolujących pozwala na bardziej elastyczne dopasowanie do danych, szczególnie gdy relacja między zmiennymi jest skomplikowana i nie jest dobrze modelowana przez proste funkcje liniowe.

cv=TRUE: Opcja ta włącza walidację krzyżową. Używając walidacji krzyżowej, funkcja automatycznie wybiera wartość parametru lambda, która minimalizuje błąd predykcji, zapewniając optymalne wygładzenie.

Do zmiennej prestige dodano niewielki szum za pomocą rnorm(length(prestige), 0, 0.01). Jest to zabieg, który pozwala uniknąć problemów związanych z danymi o identycznych wartościach x, które mogłyby wpłynąć na stabilność dopasowania. Wartości x muszą być unikalne lub przynajmniej różnić się minimalnie, aby funkcja smooth.spline() mogła działać poprawnie, ponieważ funkcja ta jest przeznaczona dla zestawów danych o różnorodnych wartościach na osi x.

smr <- smooth.spline(prestige + rnorm(length(prestige), 0, 0.01), income, cv = TRUE)
smr <- data.frame(x=smr$x,y=smr$y)
ggplot(Prestige) +
  geom_point(aes(x=prestige,y=income)) +
  ggtitle("Prestige (Sploty interpolujące, lambda wybrana przez CV)") +
  geom_line(data=smr, aes(x=x, y=y), col='purple')  +
  theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

1.3.1 Sploty naturalne

Naturalne sploty to jedna z form interpolacji splajnowej, która szczególnie dobrze nadaje się do modelowania nieliniowych zależności w danych o naturalnych granicach. W odróżnieniu od splajnowych dopasowań bez ograniczeń, sploty naturalne zakładają, że krzywa dopasowania jest liniowa na krańcach przedziału, co pozwala uniknąć nienaturalnych wygięć na końcach danych. Jest to przydatne w sytuacjach, gdy wartości jednej ze zmiennych mają ograniczenia, jak np. income (dochód), który nie może być ujemny.

fit <- lm(prestige ~ ns(income, df=6), Prestige)
ggplot(Prestige) +
    geom_point(aes(x=income,y=prestige)) +
    ggtitle("Prestige wzgl. dochodu (Naturalne sploty, 6 df)") +
    geom_line(aes(x=income, y=fitted(fit)), col='purple') +
    theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

Na początku linia trendu wzrasta stosunkowo stromo, co sugeruje, że wzrost dochodu w niższym zakresie ma znaczący wpływ na wzrost prestiżu zawodów.W środkowym zakresie dochodów linia trendu wykazuje największe przyspieszenie, co wskazuje na to, że w tym segmencie dochodowym występuje największa zmienność w postrzeganiu prestiżu. W wyższym zakresie dochodów linia trendu staje się bardziej płaska, co może oznaczać, że po pewnym punkcie krytycznym dalsze wzrosty dochodu nie przynoszą proporcjonalnych zysków w postrzeganym prestiżu.

fit <- lm(prestige ~ ns(income, df=12), Prestige)
ggplot(Prestige) +
    geom_point(aes(x=income,y=prestige)) +
    ggtitle("Prestige wzgl. dochodu (Naturalne sploty, 12 df)") +
    geom_line(aes(x=income, y=fitted(fit)), col='purple') +
     theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

1.4 LOESS

Najbardziej znaną implementacją jest loess (lokalnie kwadratowy). Używa się tutaj jądra tri-cube i zmiennej szerokości pasma. Span kontroluje szerokość pasma, a degree jest rzędem wielomianu. Można użyć family=“symmetric” dla solidnego dopasowania.

smr <- loess(prestige ~ income, data=Prestige, span=0.15, degree=2, family="symmetric")
ggplot(Prestige) +
  geom_point(aes(x=income, y=prestige), color="black") +  
  geom_line(aes(x=income, y=fitted(smr)), color='purple') + 
  ggtitle("Prestige względem dochodu (Loess: span=0.15, degree=2, robust fitting)") +
       theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

Niski dochód: Na początku linia LOESS dobrze odwzorowuje lokalne minimum, jednak jest kilka obszarów, gdzie linia znacznie odchyla się od punktów danych, szczególnie przy bardzo niskich dochodach. Średni dochód: W tym zakresie linia dosyć dobrze śledzi główne trendy danych, choć występują punkty, które są daleko od przewidywanej linii trendu. Wysoki dochód: Linia wydaje się być płaska i nie odzwierciedla dobrze dystrybucji danych w tym zakresie, gdzie zmienność danych jest znaczna.

2 Zad 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 wykorzystana została metoda regresji nieparametrycznej.Na początku wczytujemy dane.

data("mcycle")
attach(mcycle)

2.1 Wykres ukazujący relację między czasem, a przyśpieszeniem.

ggplot(mcycle, aes(x = times, y = accel)) +
  geom_point(color = "darkgreen", size = 2, alpha = 0.6) +  
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "red") +  
  labs(
    title = "Relacja między czasem a przyśpieszeniem",
    caption = "Dane źródłowe: Pakiet 'MASS'"
  ) +
  theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

2.2 Estymacja nieliniowej zależności

Poniżej wykorzystano metodę locpoly() z pakietu KernSmooth do estymacji nieliniowej zależności, gdzie stopień polinomu wynosi 2, a szerokość pasma ustalono na pięć.

fit2 <- locpoly(times, accel,
                degree=2, 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("Relacja między czasem a przyspieszeniem")+
  theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

2.3 Interpolacja splotów

times_unique <- times + rnorm(length(times), 0, 1e-6)
smr1 <- smooth.spline(times_unique, accel, cv=TRUE)

## Warning in smooth.spline(times_unique, accel, cv = TRUE): cross-validation with
## non-unique 'x' values seems doubtful

smr1 <- data.frame(x=smr1$x,y=smr1$y)
ggplot(mcycle) +
  geom_point(aes(x=times,y=accel)) +
  ggtitle("Sploty interpolujące, lambda wybrana przez CV") +
  geom_line(data=smr1, aes(x=x, y=y), col='red') +
    theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5)  
  )

Dane wykazują wyraźne oscylacje przyspieszenia w zależności od czasu. Na początku przyspieszenie spada, osiągając minimum około czasu 20, po czym gwałtownie wzrasta do maksimum około czasu 30. Następnie ponownie spada i stabilizuje się na poziomie zbliżonym do zera, wskazując na wygaszenie oscylacji.

Krzywa uzyskana za pomocą metody smooth spline skutecznie odwzorowuje ogólną tendencję w danych, dostosowując się do zmienności przyspieszenia w miarę upływu czasu. Wybór wartości parametru lambda poprzez walidację krzyżową pozwala na optymalne wygładzenie, co oznacza, że krzywa dobrze oddaje ogólny wzorzec bez nadmiernego dopasowania do lokalnych fluktuacji.

Wykorzystanie spline’ów pozwala na wychwycenie i oddanie zmian w przyspieszeniu, uwzględniając zarówno gwałtowne wzrosty i spadki, jak i stabilizację na końcu. Dzięki temu krzywa jest płynna, ale nie traci szczegółów charakterystycznych dla danych, co daje pełniejszy obraz dynamiki przyspieszenia w funkcji czasu.