Na początku wczytujemy dane.
rowery<-read_excel("sklep_rowerowy.xlsx")
SKLEP ROWEROWY
Analizowany zbiór danych zawiera informacje na temat klientów sklepu rowerowego. Z jego pomocą możemy uzyskać informacje na temat wieku, dochodu, płci, liczby posiadanych dzieci czy np. regionu pochodzenia klientów sklepu. Celem projektu jest przeprowadzenie kompleksowej analizy danych, która obejmować będzie kilka kluczowych etapów, a mianowicie: czyszczenie danych, wizualizację, analizę opisową i wnioskowanie statystyczne.
I ETAP - CZYSZCZENIE DANYCH
1. Przegląd danych
Po wyświetleniu początkowych wierszy obserwujemy m.in. spacje w nazwach kolumn naszych zmiennych. Aby się ich pozbyć korzystamy z funkcji clean_names z pakietu ‘janitor’, która służy do przekształcania nazw kolumn w ramce danych w formę bardziej czytelną i zgodą z konwencją.
head(rowery)
## # A tibble: 6 × 13
## ID `Marital Status` Gender Income Children Education Occupation
## <dbl> <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 12496 Married Female 40000 1 Bachelors Skilled Manual
## 2 24107 Married Male 30000 3 Partial College Clerical
## 3 14177 Married Male 80000 5 Partial College Professional
## 4 24381 Single <NA> 70000 0 Bachelors Professional
## 5 25597 Single Male 30000 0 Bachelors Clerical
## 6 13507 Married Female 10000 2 Partial College Manual
## # ℹ 6 more variables: `Home Owner` <chr>, Cars <dbl>, `Commute Distance` <chr>,
## # Region <chr>, Age <dbl>, `Purchased Bike` <chr>
rowery <- clean_names(rowery)
1. Sprawdzenie występowania brakujących obserwacji NA
Poniżej przeprowadzona została analiza występowania wartości brakujących dla poszczególnych zmiennych. Badanie obecności wartości NA jest ważne dla zapewnienia poprawności analizy, utrzymania jakości danych i generowania precyzyjnych wyników. To kluczowy etap w procesie przetwarzania i przygotowywania danych do dalszej analizy. Poprawne zarządzanie wartościami brakującymi pozwala na uniknięcie błędnych interpretacji oraz zapewnia solidne fundamenty do analiz statystycznych czy modelowania danych.
sum(complete.cases(rowery))
## [1] 952
W naszej bazie danych zidentyfikowaliśmy 952 wiersze, w których wszystkie dane są kompletnie wypełnione. Niemniej jednak, istnieje 48 wierszy, w których co najmniej jedna wartość jest brakująca (oznaczona jako NA). Oznacza to, że wiersze te zawierają pewną ilość niekompletnych informacji. Procentowy udział wierszy z brakującymi danymi w stosunku do ogólnej liczby rekordów wynosi 4.8%.
manyNAs(rowery)
## [1] 689
Widzimy, że najwięcej braków danych zawiera wiersz 689.
sum(is.na(rowery))
## [1] 53
Zauważmy, że ogółem mamy 53 braki danych w bazie. Biorąc pod uwagę, że liczba wierszy z brakującymi danymi wynosi 48, można stwierdzić, że niektóre rekordy posiadają więcej niż jedną wartość brakującą.
miss_var_summary(rowery)
## # A tibble: 13 × 3
## variable n_miss pct_miss
## <chr> <int> <dbl>
## 1 gender 11 1.1
## 2 cars 9 0.9
## 3 children 8 0.8
## 4 age 8 0.8
## 5 marital_status 7 0.7
## 6 income 6 0.6
## 7 home_owner 4 0.4
## 8 id 0 0
## 9 education 0 0
## 10 occupation 0 0
## 11 commute_distance 0 0
## 12 region 0 0
## 13 purchased_bike 0 0
Widzimy, że braki wystepują jedynie dla zmiennych gender (11 braków), cars (9 braków), children (8 braków), age (8 braków), marital_status (7 braków), income (6 braków) i home_owner (4 braki). Pozostałe zmienne nie posiadają żadnych wartości brakujących.
rowery %>%
miss_case_table()
## # A tibble: 4 × 3
## n_miss_in_case n_cases pct_cases
## <int> <int> <dbl>
## 1 0 952 95.2
## 2 1 44 4.4
## 3 2 3 0.3
## 4 3 1 0.1
Mamy 44 wiersze, w których brak jest jednej wartości, 3 wiersze z 2 brakami i jeden wiersz z 3 wartościami NA.
which(is.na(rowery$gender))
## [1] 4 155 336 602 689 696 868 909 952 974 998
Dla zmiennej ‘gender’ braki danych występują w wierszach: 4, 155, 336, 602, 689, 696, 868, 909, 952, 974, 998.
which(is.na(rowery$cars))
## [1] 13 197 203 352 449 512 562 616 934
Dla zmiennej ‘cars’ braki danych występują w wierszach: 13, 197, 203, 352, 449, 512, 562, 616, 934.
which(is.na(rowery$children))
## [1] 118 218 387 550 639 689 806 961
Dla zmiennej ‘children’ braki danych występują w wierszach: 118, 218, 387, 550, 639, 689, 806, 961.
which(is.na(rowery$age))
## [1] 10 99 226 372 555 689 771 987
Dla zmiennej ‘age’ braki danych występują w wierszach: 10, 99, 226, 372, 555, 689, 771, 987.
which(is.na(rowery$marital_status))
## [1] 9 28 50 99 151 235 302
Dla zmiennej ‘marital_status’ braki danych występują w wierszach: 9, 28, 50, 99, 151, 235, 302.
which(is.na(rowery$income))
## [1] 10 111 192 302 442 510
Dla zmiennej ‘income’ braki danych występują w wierszach: 10, 111, 192, 302, 442, 510.
which(is.na(rowery$home_owner))
## [1] 7 366 647 944
Dla zmiennej ‘home_owner’ braki danych występują w wierszach: 7, 366, 647, 944.
2. Sprawdzenie typów danych w każdej kolumnie
Sprawdzenie typów danych w każdej kolumnie jest kluczowe dla zachowania poprawności, efektywności oraz precyzji analizy danych. Różne typy danych wymagają różnych operacji, a znajomość ich rodzaju pozwala na zastosowanie odpowiednich metod przetwarzania.
data_class <- data.frame(class = sapply(rowery, class))
data_class
## class
## id numeric
## marital_status character
## gender character
## income numeric
## children numeric
## education character
## occupation character
## home_owner character
## cars numeric
## commute_distance character
## region character
## age numeric
## purchased_bike character
3. Sprawdzenie wiarygodności danych
Poprzez wywołanie podstawowych statystyk możemy sprawdzić, czy nasze dane wydają się być wiarygodne, tzn. czy średnia ilość posiadanych dzieci, czy samochodów wydaje się być prawdopodobna. Ponadto użycie ‘view(dfSummary(rowery))’ dostarcza nam wielu cennych informacji m.in.:
- widzimy, że każdy klient posiada unikalne id,
- otrzymujemy podsumowanie braków obserwacji dla każdej zmiennej,
- dla zmiennych jakościowych widzimy możliwe kategorie wyboru wraz z informacją o częstości jej wystąpienia.
dfSummary(rowery)
## Data Frame Summary
## rowery
## Dimensions: 1000 x 13
## Duplicates: 0
##
## -------------------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ------------------ ------------------------------- ---------------------- --------------------- ---------- ---------
## 1 id Mean (sd) : 19966 (5347.3) 1000 distinct values : . . . 1000 0
## [numeric] min < med < max: : : : : : : : : : : (100.0%) (0.0%)
## 11000 < 19744 < 29447 : : : : : : : : : :
## IQR (CV) : 9180 (0.3) : : : : : : : : : :
## : : : : : : : : : :
##
## 2 marital_status 1. Married 535 (53.9%) IIIIIIIIII 993 7
## [character] 2. Single 458 (46.1%) IIIIIIIII (99.3%) (0.7%)
##
## 3 gender 1. Female 489 (49.4%) IIIIIIIII 989 11
## [character] 2. Male 500 (50.6%) IIIIIIIIII (98.9%) (1.1%)
##
## 4 income Mean (sd) : 56267.6 (31067.8) 16 distinct values : 994 6
## [numeric] min < med < max: : . . (99.4%) (0.6%)
## 10000 < 60000 < 170000 . : : :
## IQR (CV) : 40000 (0.6) : : : :
## : : : : : . .
##
## 5 children Mean (sd) : 1.9 (1.6) 0 : 274 (27.6%) IIIII 992 8
## [numeric] min < med < max: 1 : 169 (17.0%) III (99.2%) (0.8%)
## 0 < 2 < 5 2 : 209 (21.1%) IIII
## IQR (CV) : 3 (0.9) 3 : 133 (13.4%) II
## 4 : 126 (12.7%) II
## 5 : 81 ( 8.2%) I
##
## 6 education 1. Bachelors 306 (30.6%) IIIIII 1000 0
## [character] 2. Graduate Degree 174 (17.4%) III (100.0%) (0.0%)
## 3. High School 179 (17.9%) III
## 4. Partial College 265 (26.5%) IIIII
## 5. Partial High School 76 ( 7.6%) I
##
## 7 occupation 1. Clerical 177 (17.7%) III 1000 0
## [character] 2. Management 173 (17.3%) III (100.0%) (0.0%)
## 3. Manual 119 (11.9%) II
## 4. Professional 276 (27.6%) IIIII
## 5. Skilled Manual 255 (25.5%) IIIII
##
## 8 home_owner 1. No 314 (31.5%) IIIIII 996 4
## [character] 2. Yes 682 (68.5%) IIIIIIIIIIIII (99.6%) (0.4%)
##
## 9 cars Mean (sd) : 1.5 (1.1) 0 : 238 (24.0%) IIII 991 9
## [numeric] min < med < max: 1 : 267 (26.9%) IIIII (99.1%) (0.9%)
## 0 < 1 < 4 2 : 342 (34.5%) IIIIII
## IQR (CV) : 1 (0.8) 3 : 85 ( 8.6%) I
## 4 : 59 ( 6.0%) I
##
## 10 commute_distance 1. 0-1 Miles 366 (36.6%) IIIIIII 1000 0
## [character] 2. 1-2 Miles 169 (16.9%) III (100.0%) (0.0%)
## 3. 10+ Miles 111 (11.1%) II
## 4. 2-5 Miles 162 (16.2%) III
## 5. 5-10 Miles 192 (19.2%) III
##
## 11 region 1. Europe 300 (30.0%) IIIIII 1000 0
## [character] 2. North America 508 (50.8%) IIIIIIIIII (100.0%) (0.0%)
## 3. Pacific 192 (19.2%) III
##
## 12 age Mean (sd) : 44.2 (11.4) 53 distinct values : 992 8
## [numeric] min < med < max: : : : (99.2%) (0.8%)
## 25 < 43 < 89 : : : : :
## IQR (CV) : 17 (0.3) : : : : : .
## : : : : : : :
##
## 13 purchased_bike 1. No 519 (51.9%) IIIIIIIIII 1000 0
## [character] 2. Yes 481 (48.1%) IIIIIIIII (100.0%) (0.0%)
## -------------------------------------------------------------------------------------------------------------------------
#view(dfSummary(rowery))
descr(rowery)
## Non-numerical variable(s) ignored: marital_status, gender, education, occupation, home_owner, commute_distance, region, purchased_bike
## Descriptive Statistics
## rowery
## N: 1000
##
## age cars children id income
## ----------------- -------- -------- ---------- ---------- -----------
## Mean 44.18 1.46 1.91 19965.99 56267.61
## Std.Dev 11.36 1.12 1.63 5347.33 31067.82
## Min 25.00 0.00 0.00 11000.00 10000.00
## Q1 35.00 1.00 0.00 15289.50 30000.00
## Median 43.00 1.00 2.00 19744.00 60000.00
## Q3 52.00 2.00 3.00 24475.50 70000.00
## Max 89.00 4.00 5.00 29447.00 170000.00
## MAD 11.86 1.48 1.48 6848.13 29652.00
## IQR 17.00 1.00 3.00 9180.00 40000.00
## CV 0.26 0.77 0.85 0.27 0.55
## Skewness 0.52 0.42 0.39 0.05 0.75
## SE.Skewness 0.08 0.08 0.08 0.08 0.08
## Kurtosis -0.27 -0.41 -1.02 -1.19 0.50
## N.Valid 992.00 991.00 992.00 1000.00 994.00
## Pct.Valid 99.20 99.10 99.20 100.00 99.40
4. Sprawdzanie spełnienia pewnych reguł dla zbioru danych
Posiadając już podstawową wiedzę na temat naszych danych, chcemy sprawdzić kilka podstawowych reguł:
- czy zmienna wiek na pewno wszędzie przyjmuje wartości dodatnie,
- czy liczba posiadanych dzieci nie jest nigdzie ujemna,
- czy liczba posiadanych samochodów nie jest ujemna,
- czy płeć przyjmuje jedną z dwóch kategorii ‘kobieta’ i ‘mężczyzna’.
Najpierw stwórzmy zbiór powyżej opisanych reguł.
rules <- validator(age > 0, gender %in% c('Female','Male')
, income >= 0, children >= 0, cars >= 0)
rules
## Object of class 'validator' with 5 elements:
## V1: age > 0
## V2: gender %in% c("Female", "Male")
## V3: income >= 0
## V4: children >= 0
## V5: cars >= 0
Sprawdźmy teraz, czy są one spełnione w naszym zbiorze danych.
cf <- confront(rowery, rules, key="id")
summary(cf)
## name items passes fails nNA error warning expression
## 1 V1 1000 992 0 8 FALSE FALSE age > 0
## 2 V2 1000 989 0 11 FALSE FALSE gender %vin% c("Female", "Male")
## 3 V3 1000 994 0 6 FALSE FALSE income - 0 >= -1e-08
## 4 V4 1000 992 0 8 FALSE FALSE children - 0 >= -1e-08
## 5 V5 1000 991 0 9 FALSE FALSE cars - 0 >= -1e-08
Spełnienie reguł możemy również zwizualizować.
barplot(cf, main="rowery")
## Warning: The 'barplot' method for confrontation objects is deprecated. Use
## 'plot' instead
Moglibyśmy to również wykonać w inny sposób, który został przedstawiony poniżej. W razie niespełnienia którejkolwiek z reguł, błędna wartość zostałaby przekształcona na wartość ‘NA’.
RULE <- editset(c("age > 0","gender %in% c('Female','Male')"
, "income >= 0", "children >= 0", "cars >= 0"))
summary(violatedEdits(RULE, rowery))
## Edit violations, 1000 observations, 0 completely missing (0%):
##
## editname freq rel
## dat1 11 1.1%
##
## Edit violations per record:
##
## errors freq rel
## 0 961 96.1%
## 1 37 3.7%
## 2 1 0.1%
## 3 1 0.1%
rowery[localizeErrors(RULE, rowery)$adapt] <- NA
Okazuje się, że wszystkie reguły w naszym zbiorze danych są spełnione. Występujące błędy wynikają jedynie z występowania wartości brakujących, którymi zajmiemy się w następnym kroku.
2. Imputacja danych
Z analiz przeprowadzonych powyżej wynika, że w naszym zbiorze danych zlokalizowano 53 braki danych. W celu skutecznego zarządzania tymi brakami, podjęto decyzję o nieusuwanie żadnego z wierszy, a zamiast tego zastosowano różne metody imputacji wartości NA. Wybór konkretnych metod zależał od natury danych.
1. Imputacja-zmienne ilościowe
Braki w zmiennych dotyczących dochodu, wieku, liczby posiadanych dzieci i samochodów wypełnione zostały za pomocą średniej. Ze względu na specyfikę danych ‘age’, ‘children’ i ‘cars’, wartości te nie powinny mieć żadnych miejsc po przecinku, gdyż nie możemy posiadać np. 1.5 dziecka.
rowery[is.na(rowery$income), "income"] <- mean(rowery$income, na.rm = T)
rowery[is.na(rowery$age), "age"] <- round(mean(rowery$age, na.rm = T), digits=0)
rowery[is.na(rowery$cars), "cars"] <- round(mean(rowery$cars, na.rm = T), digits=0)
rowery[is.na(rowery$children), "children"] <-round(mean(rowery$children, na.rm = T), digits=0)
2. Imputacja-zmienne jakościowe
Po uzupełnieniu braków w zmiennych ilościowych, przystępujemy teraz do korekty brakujących danych w przypadku zmiennych jakościowych. Ponieważ nasze zmienne jakościowe zawierają jedynie dwie kategorie, planujemy zidentyfikować, która z tych kategorii jest dominująca. Następnie, braki danych w tych zmiennych zostaną wypełnione wartościami dominującymi, czyli po prostu dominantą.
Zacznijmy od zmiennej ‘home_owner’. Widzimy, że przyjmuje ona dwie główne kategorie wyboru ‘Yes’ oraz ‘No’. 68,2% klientów sklepu rowerowego posiada dom, więc braki w danych uzupełniamy wartością ‘Yes’.
unique(rowery$home_owner)
## [1] "Yes" "No" NA
rowery %>%
group_by(home_owner) %>%
count(home_owner) %>%
summarize(proporcja=n/1000 * 100)
## # A tibble: 3 × 2
## home_owner proporcja
## <chr> <dbl>
## 1 No 31.4
## 2 Yes 68.2
## 3 <NA> 0.4
rowery[is.na(rowery$home_owner), "home_owner"] <- "Yes"
Dla zmiennej ‘marital_status’ mamy dwie kategorie wyboru ‘Married’ i ‘Single’. 53,5% klientów posiada męża bądź żonę, dlatego braki w danych wypełniamy wartością ‘Married’.
unique(rowery$marital_status)
## [1] "Married" "Single" NA
rowery %>%
group_by(marital_status) %>%
count(marital_status) %>%
summarize(proporcja=n/1000 * 100)
## # A tibble: 3 × 2
## marital_status proporcja
## <chr> <dbl>
## 1 Married 53.5
## 2 Single 45.8
## 3 <NA> 0.7
rowery[is.na(rowery$marital_status), "marital_status"] <- "Married"
Dla zmiennej ‘gender’ mamy dwie kategorie wyboru ‘Male’ i ‘Female’. 50% klientów sklepu to mężczyźni, więc braki w danych uzupełniamy wartością ‘Male’.
unique(rowery$gender)
## [1] "Female" "Male" NA
rowery %>%
group_by(gender) %>%
count(gender) %>%
summarize(proporcja=n/1000 * 100)
## # A tibble: 3 × 2
## gender proporcja
## <chr> <dbl>
## 1 Female 48.9
## 2 Male 50
## 3 <NA> 1.1
rowery[is.na(rowery$gender), "gender"] <- "Male"
Upewnijmy się, że wszystko zadziałało, a w naszym zbiorze danych nie występują już żadne wartości brakujące.
sum(is.na(rowery))
## [1] 0
3. Obserwacje odstające
Sprawdzanie zbioru danych pod kątem wartości odstających jest kluczowe dla utrzymania jakości analizy danych i poprawnego zrozumienia badanego zjawiska. W przypadku identyfikacji wartości odstających, istnieją różne metody ich obsługi, takie jak usuwanie, transformacja, czy stosowanie bardziej zaawansowanych technik modelowania.
Przechodzimy więc do sprawdzenia naszego zbioru danych pod względem występowania wartości odstających.
1. Zlokalizowanie wartości odstających
Możemy przykładowo stworzyć funkcję, która będzie wykrywać obserwacje odstające. Funkcja ta opiera się na metodzie Interquartile Range, gdzie za odstające wartości uznaje się te, które leżą poza pewnym zakresem określonym na podstawie pierwszego i trzeciego kwartyla.
find_outliers <- function(rowery, k = 1.5) {
quantiles <- quantile(rowery, c(0.25, 0.5, 0.75))
diff <- k * (quantiles[3] - quantiles[1])
lb <- quantiles[1] - diff
ub <- quantiles[3] + diff
is_outlier <- rowery[which(rowery <lb | rowery > ub)]
return(is_outlier)
}
find_outliers(rowery$income)
## [1] 160000 170000 170000 150000 160000 150000 160000 150000 170000 150000
Inną metodą jest funkcja ‘boxplot.stats’ za pomocą której możemy wyświetlić wartości odstające. Widzimy, że mamy aż 10 wartości odstających dla zmiennej ‘income’.
out <-boxplot.stats(rowery$income)$out
out
## [1] 160000 170000 170000 150000 160000 150000 160000 150000 170000 150000
Możemy również skorzystać z testów statystycznych. Analizując uzyskaną wartość p-value dla testu Grubbs’a przyjmujemy, że najwyższa wartość dla zmiennej ‘income’ równa 170 000 jest odstająca.
test <- grubbs.test(rowery$income)
test
##
## Grubbs test for one outlier
##
## data: rowery$income
## G = 3.67182, U = 0.98649, p-value = 0.1151
## alternative hypothesis: highest value 170000 is an outlier
Przeprowadzamy test Grubbs’a dla najmniejszej wartości w zbiorze ‘income’.
test1 <- grubbs.test(rowery$income, opposite = TRUE)
test1
##
## Grubbs test for one outlier
##
## data: rowery$income
## G = 1.49374, U = 0.99776, p-value = 1
## alternative hypothesis: lowest value 10000 is an outlier
Zwizualizujmy teraz wartości odstające na wykresie wraz z opisem, które z nich są odstające.
boxplot(rowery$income, col = "blue",
ylab = "income",
main = "Wykres pudełkowy dochodu")
mtext(paste("Outliers: ", paste(out, collapse = ", ")))
Sprawdźmy jeszcze występowanie wartości odstających dla zmiennej ‘age’. Okazuje się, że mamy 4 wartości odstającej dla tej zmiennej.
boxplot.stats(rowery$age)$out
## [1] 78 89 80 78
Wykonajmy jeszcze dodatkowo test Grubbs’a.
test <- grubbs.test(rowery$age)
test
##
## Grubbs test for one outlier
##
## data: rowery$age
## G = 3.96061, U = 0.98428, p-value = 0.03514
## alternative hypothesis: highest value 89 is an outlier
2. Przekształcenie wartości odstających
Zdecydowaliśmy się przekształcić wartości odstające dla zmiennych “age” i “income” za pomocą metody capping. Metoda capping polega na ustaleniu górnej (max) i dolnej (min) granicy wartości dla danej zmiennej, a następnie przypisaniu wszystkim wartościom przekraczającym te granice wartości skrajnych.
qnt <- quantile(rowery$income, probs=c(.25, .75), na.rm = T)
caps <- quantile(rowery$income, probs=c(.05, .95), na.rm = T)
H <- 1.5 * IQR(rowery$income, na.rm = T)
rowery$income[rowery$income < (qnt[1] - H)] <- caps[1]
rowery$income[rowery$income > (qnt[2] + H)] <- caps[2]
Sprawdźmy, czy wszystkie wartości odstające dla ‘income’ zostały przekształcone.
boxplot.stats(rowery$income)$out
## numeric(0)
Tym samym sposobem przekształćmy wartości odstające dla zmiennej ‘age’.
qnt <- quantile(rowery$age, probs=c(.25, .75), na.rm = T)
caps <- quantile(rowery$age, probs=c(.05, .95), na.rm = T)
H <- 1.5 * IQR(rowery$age, na.rm = T)
rowery$age[rowery$age < (qnt[1] - H)] <- caps[1]
rowery$age[rowery$age > (qnt[2] + H)] <- caps[2]
Widzimy, że brak jest już wartości odstających dla zmiennej ‘age’.
boxplot.stats(rowery$age)$out
## numeric(0)
Kończymy pierwszy etap naszego projektu. Dzięki procesowi czyszczenia danych osiągnięto spójność w zbiorze danych. Brakujące wartości zostały uzupełnione, a błędy zostały skorygowane, co pozwali teraz na bardziej precyzyjną analizę. Otrzymaliśmy solidne fundamenty do przeprowadzenia bardziej zaawansowanych etapów projektu, takich jak wizualizacje, analiza opisowa, modelowanie czy wnioskowanie statystyczne.
II ETAP-WIZUALIZACJE
Wizualizacje stanowią potężne narzędzie w każdym projekcie, wspomagają efektywną komunikację, lepsze zrozumienie danych oraz sprzyjają szybszemu i bardziej kreatywnemu podejmowaniu decyzji. Poniżej przedstawiono kilka ciekawych wizualizacji dla analizowanego przez nas zbioru danych.
1.Średni dochód w zależności od wieku i płci
Prezentowany wykres składa się z dwóch paneli, z których lewy reprezentuje dane dotyczące kobiet, natomiast prawy mężczyzn. Na obu panelach zaznaczone są czerwone punkty odzwierciedlające empiryczne dane dotyczące średnich dochodów w poszczególnych grupach wiekowych. Linie trendu w kolorze niebieskim zostały dopasowane do punktów, ukazując ogólny trend wzrostu średnich dochodów w miarę postępującego wieku. Warto zauważyć, że panel mężczyzn charakteryzuje się większym zróżnicowaniem (różnicą) średnich dochodów, jednocześnie prezentując wyższą średnią wartość dochodu. Dodatkowo, dla mężczyzn zauważalny jest silniejszy trend wzrostu średniego dochodu wraz z wiekiem. To sugeruje, że w miarę upływu lat różnice w dochodach między różnymi grupami wiekowymi mężczyzn mogą się nasilać, tworząc bardziej zauważalne tendencje wzrostowe.
dw1<-rowery%>%
group_by(age, gender)%>%
summarize(meaninc=mean(income))
## `summarise()` has grouped output by 'age'. You can override using the `.groups`
## argument.
ggplot(dw1, aes(age, meaninc))+
geom_line()+
geom_point(colour="red")+
facet_wrap(~gender)+labs(x="Wiek", y="Średnie dochody")+geom_smooth(method="lm")+ggtitle("Średni dochód w zależności od wieku oraz płci")
## `geom_smooth()` using formula = 'y ~ x'
2. Występowanie obserwacji o określonym wieku i dochodzie
Na poniższym wykresie przedstawiony został wykres punktowy. Każdy punkt na wykresie reprezentuje jeden przypadek. Wielkość punktu z kolei odpowiada częstości występowania obserwacji o określonym wieku i poziomie dochodów. Analizując wykres, można stwierdzić, że większość obserwacji skupia się w przedziale wiekowym 30-60 lat i niższym przedziale dochodowym.
ggplot(rowery, aes(age, income))+
stat_sum(alpha=0.4)+
scale_size(range=c(1,5))+
ggtitle("Występowanie obserwacji o określonym wieku oraz poziomie dochodów")+labs(x="Wiek", y="Dochód", size="Występowanie")
3. Ilość posidanych samochodów w zależności od regionu
Na poniższym wykresie pokazany jest stosunek liczby posiadanych samochodów w gospodarstwie domowym ze względu na region, z którego pochodzą klienci sklepu rowerowego. Rzuca sie w oczy fakt, że w Ameryce Północnej zdecydowanie najwięcej osób posiada 2 samochody. Region ten posiada również największy odsetek gospodarstw domowych z jednym samochodem. Najwięcej klientów nieposiadających samochodu zamieszkuje Europę. Na Pacyfiku jest tam najmniejszy odsetek gospodarstw bez samochodu. Ten region charakteryzuje się również najmniejszą dysproporcją jeżeli chodzi o liczbę samochodów.
ggplot(rowery, aes(x=cars, fill=region))+
geom_bar(position="dodge")+
labs(x="Ilość posiadanych samochodów", y="Występowanie", fill="Region")
4.Zagęszczenie obserwacji o określonym wieku w zależności od regionu
Poniższy wykres przedstawia rozkład wieku dla trzech różnych regionów: Europy, Ameryki Półnicnej i Pacyfiku. Oś pozioma (x) reprezentuje wiek, a oś pionowa (y) reprezentuje zagęszczenie, czyli estymowaną gęstość prawdopodobieństwa obserwacji w każdym punkcie osi wieku. W Europie największe zagęszczenie obserwacji występuje wśród osób w średnim wieku, z szczytem około 40 lat. Dla Ameryki Północnej rozkład wieku jest szerszy, z szczytem około 50 lat. Region Pacyfiku wydaje się mieć najszerszy rozkład wieku z największym zagęszczeniem wśród osób w wieku około 30 lat i stopniowym spadkiem w kierunku wyższych grup wiekowych. Rozkład wieku w Europie jest bardziej szpiczasty i mniej rozproszony niż w Ameryce Północnej, co może sugerować mniejszą różnorodność wieku w tym regionie. Linie obrysujące każdą z krzywych gęstości dają wyobrażenie o kształcie rozkładu dla każdego regionu, pozwalając dostrzec, gdzie znajdują się największe skupiska obserwacji.
ggplot(rowery, aes(age, fill=region))+
geom_density(alpha=0.25)+
labs(x="Wiek",y="Zagęszczenie", fill="Region")+
ggtitle("Zagęszczenie obserwacji o określonym poziomie wieku w zależności od regionu")
5. Obserwacje o określonym poziomie dochodów w zależności od rodzaju zawodu
Wykres poniżej jednoznacznie ilustruje, że dochody w poszczególnych grupach zawodowych są istotnie zróżnicowane. To zjawisko może być rezultatem różnic w wymaganiach, umiejętnościach, stopniu odpowiedzialności czy też specyfice poszczególnych zawodów. Pracownicy fizyczni najczęściej osiągają najniższe wynagrodzenia, podczas gdy pracownicy biurowi zarabiają nieco więcej. Wykwalifikowani pracownicy fizyczni cieszą się średnimi dochodami. Grupy zawodowe obejmujące pracowników profesjonalnych oraz na stanowiskach kierowniczych prezentują najwyższe wynagrodzenia. To zwykle osoby o zaawansowanych kwalifikacjach, odpowiedzialnościach zarządczych czy specjalistycznych umiejętnościach, co znajduje odzwierciedlenie w ich dochodach.
ggplot(rowery, aes(income, fill=occupation))+
geom_histogram(bins=20, color='black')+
facet_wrap(~occupation, nrow=5)+
labs(x="Dochód", y="Występowanie", fill="Praca")+
theme(legend.position = 'none')+
ggtitle("Wystepowanie obserwacji o określonym poziomie dochodów w zależności od rodzaju zawodu")
6. Wizualizacja poziomu dochodu w zależności od płci
Wykres poniżej pokazuje medianę dochodu, a także jego maksymalny, minimalny i średni poziom w zależności od płci. Mediana i średnie dochody są przedstawione na podobnym poziomie dla obu płci, z lekkim przesunięciem w kierunku wyższych wartości dla mężczyzn. Rozstęp dochodów, czyli różnica między kwartylem górnym a dolnym, również wydaje się być podobny dla obu grup.
ggplot(rowery, aes(gender, income))+
geom_boxplot(aes(fill=gender))+
stat_boxplot(geom="errorbar",position="dodge")+
stat_summary(aes(ymin=after_stat(y),ymax=after_stat(y)),fun=mean)+ggtitle("Mediana dochodów, maksymalny, minimalny oraz średni dochód w zależności od płci")+theme(legend.position='none')
7. Ilość posiadanych samochodów w zależności od posiadania roweru
Wykres poniżej może sugerować, że osoby posidające rower mają tendencję do posiadania mniejszej ilości samochodów niż osoby, które roweru nie posiadają.
ggplot(rowery, aes(cars, fill=purchased_bike))+geom_histogram(bins=5, position="dodge", color='black')+facet_wrap(~purchased_bike)+theme(legend.position='none')+labs(x="Ilość posiadanych samochodów", y="Występowanie")+ggtitle("Ilość posiadanych samochodów w zależności od posiadania roweru")
8. Występowanie obserwacji o poszczególnym rodzaju wykonywanej pracy według wykształcenia
Poniższy rysunek przedstawia występowanie obserwacji o poszczególnym rodzaju wykonywanej pracy według wykształcenia. Najwięcej osób z wykształceniem licencjackim pracuje na stanowisku kierownicznym albo jako profesjonalista. Natomiast najmniej na stanowisku urzędniczym. Jeśli chodzi o osoby z wykształceniem wyższym to także najwięcej osób pracuje na stanowisku kierowniczym (ok. 60 obserwacji). Najmniej osób pracuje na stanowisku pracy fizycznej (ok. 10 obserwacji). Spoglądając na osoby z wykształceniem licealnym to największa ich liczba jest zatrudniona na stanowisku fachowego pracownika fizycznego, jako profesjonalista i pracownik fizyczny (odpowiednio ok. 60, 50 i 40 obserwacji). Natomiast najmniej osób z tym wykształceniem jest zatrudnionych jako menedżer lub pracownik biurowy (ok. 10 i 5 obserwacji). Analizując rodzaj wykonywanej pracy przez osoby z nieskończonymi studiami wyższymi można zauważyć, iż największy udział pracujących to profesjonaliści (ok. 60 obserwacji), a najmniejszy dotyczy menedżerów (ok. 5 obserwacji). Analizując ostatnią kategorię wykształcenia, zaobserwowano tutaj najmniejszą liczbę obserwacji. Najwięcej osób z wykształceniem niepełnym licealnym jest zatrudnionych jako pracownik manualny (ok. 30 obsrwacji), a najmniej na stanowisku menedżera.
ggplot(rowery, aes(occupation, fill=occupation))+
geom_bar(color='black')+
facet_wrap(~education, nrow=5)+
labs(x="Zatrudnienie", y="Występowanie")+
theme(legend.position= 'none')+
ggtitle("Występowanie obserwacji o poszczególnym rodzaju wykonywanej pracy według wykształcenia")
9. Średni dochód według wykształcenia oraz rodzaju pracy
Mapa cieplna jest pomocna w szybkim zrozumieniu, jak różne kombinacje wykształcenia i zawodu wpływają na dochód, gdzie ciemniejsze kolory wskazją na wyższe dochody, a jaśniejsze na niższe. Kategoria ‘Management’ wydaje się oferować najwyższe dochody niezależnie od poziomu wykształcenia. Z kolei prace biurowe i prace fizyczne mają tendencję do przynoszenia niższych dochodów.
ggplot(rowery, aes(occupation, education, fill=income))+
geom_tile(color='black')+
scale_fill_gradient2(low="white", mid="yellow", high="red")+
ggtitle("Średni dochód według wykształcenia oraz rodzaju pracy")+labs(x="Zawód", y="Wykształcenie", fill="Dochód")
10. Analiza wieku według stanu cywilnego
Obserwując wykres można zauważyć, że mediany wieku dla obu grup są różne-dla osób zamężnych/żonatych mediana jest wyższa niż dla osób stanu wolnego. Średni wiek w obu grupach wydaje się być zbliżony do mediany.
ggplot(rowery, aes(marital_status, age))+
geom_boxplot(aes(fill=marital_status))+
stat_boxplot(geom="errorbar", position="dodge")+
stat_summary(aes(ymin=..y..,ymax=..y..), fun=mean)+
labs(x="Stan cywilny", y="Wiek")+
theme(legend.position='none')+
ggtitle("Mediana, maksimum, minimum oraz średnia wieku według stanu cywilnego")
## Warning: The dot-dot notation (`..y..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(y)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
11. Stan cywilny oraz region
Zdecydowanie najwięcej osób zamężnych/żonatych wśród klientów sklepu rowerowego występowało w Ameryce Północnej. W Europie oraz na Pacyfiku występuje niemal identyczna ilość zarówno singli jak i osób po ślubie.
ggplot(rowery, aes(marital_status, fill=marital_status))+
geom_bar(position="dodge")+
facet_wrap(~region)+
theme(legend.position='none')+
labs(x="Stan cywilny", y="Ilość")+
ggtitle("Ilość osób o poszczególnym stanie cywilnym w zależności od regionu")
III ETAP-ANALIZA OPISOWA
Analiza opisowa jest nieodzownym elementem każdego projektu. Dostarcza podstawowych informacji, które stanowią punkt wyjścia do bardziej zaawansowanych analiz statystycznych i pomagają zrozumieć istotę danych.
1. Tabele liczebności
W pierwszym etapie naszej analizy pogrupujemy nasze dane w postaci prostej tabeli częstości. Napotykamy tutaj jednak pewne problemy. Przykładowo zmienna ‘income’ jest zmienną ciągłą. Dokonamy więc diskretyzacji danych-przekształcimy zmienną ciągłą na zmienną dyskretną poprzez podział zakresu wartości na przedziały.
1. Zmienna ciągła ‘income’
Dla zmiennej ‘income’ stwórzmy sześć przedziałów.
etykiety1<-c("0-25 000", "25 000-50 000", "50 000-75 000", "75 000-100 000", "100 000-125 000", "125 000-150 000")
limits1<-cut(rowery$income,seq(0,150000,by=25000),labels=etykiety1)
tabela2<-freq(limits1,type="html")
tabela2
## Frequencies
## limits1
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## --------------------- ------ --------- -------------- --------- --------------
## 0-25 000 147 14.70 14.70 14.70 14.70
## 25 000-50 000 327 32.70 47.40 32.70 47.40
## 50 000-75 000 294 29.40 76.80 29.40 76.80
## 75 000-100 000 157 15.70 92.50 15.70 92.50
## 100 000-125 000 43 4.30 96.80 4.30 96.80
## 125 000-150 000 32 3.20 100.00 3.20 100.00
## <NA> 0 0.00 100.00
## Total 1000 100.00 100.00 100.00 100.00
Możemy to również wykonać za pomocą innego sposobu. Uzyskany wskaźnik TAI jest dosyć wysoki, więc możemy zaakceptować konstrukcję tabeli częstości.
tab1<-classIntervals(rowery$income,n=6,style="fixed",fixedBreaks=seq(0,150000,by=25000))
tab1
## style: fixed
## one of 1,287 possible partitions of this variable into 6 classes
## [0,25000) [25000,50000) [50000,75000) [75000,1e+05) [1e+05,125000)
## 147 287 334 128 72
## [125000,150000]
## 32
jenks.tests(tab1)
## # classes Goodness of fit Tabular accuracy
## 6.0000000 0.9624556 0.7887639
Zwizualizujmy nasz podział.
hist(rowery$income, breaks="FD", col="green", probability = TRUE,
main="INCOME")
2. Zmienna dyskretna ‘age’
Dla zmiennej ‘age’ również utwórzmy sześć przedziałów.
etykiety<-c("20-30 lat","30-40 lat","40-50 lat","50-60 lat","60-70 lat","70-80 lat")
limits<-cut(rowery$age,seq(20,80,by=10),labels=etykiety)
tabela1<-freq(limits,type="html")
tabela1
## Frequencies
## limits
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## --------------- ------ --------- -------------- --------- --------------
## 20-30 lat 108 10.80 10.80 10.80 10.80
## 30-40 lat 313 31.30 42.10 31.30 42.10
## 40-50 lat 305 30.50 72.60 30.50 72.60
## 50-60 lat 174 17.40 90.00 17.40 90.00
## 60-70 lat 93 9.30 99.30 9.30 99.30
## 70-80 lat 7 0.70 100.00 0.70 100.00
## <NA> 0 0.00 100.00
## Total 1000 100.00 100.00 100.00 100.00
Uzyskany wskaźnik TAI jest dosyć wysoki, więc możemy zaakceptować konstrukcję tabeli częstości dla zmiennej ‘age’.
tab2<-classIntervals(rowery$age,n=6,style="fixed",fixedBreaks=seq(20,80,by=10))
tab2
## style: fixed
## one of 2,118,760 possible partitions of this variable into 6 classes
## [20,30) [30,40) [40,50) [50,60) [60,70) [70,80]
## 82 299 322 183 103 11
jenks.tests(tab2)
## # classes Goodness of fit Tabular accuracy
## 6.0000000 0.9412382 0.7494848
Zwizualizujmy nasz podział.
ggplot(rowery, aes(x = age)) +
geom_histogram(binwidth = 2, fill = "blue", color = "red", alpha = 0.8) +
labs(title = "Age Distribution", x = "Age", y = "Frequency")
3. Pozostałe tabele/wykresy liczebności
Liczebność pozostałych zmiennych możemy przedstawić na wiele różnych sposobów, zarówno za pomocą wykresów, jak i tabelek.
Education: W naszym zbiorze danych obserwujemy najwięcej klientów z wykształceniem licencjackim.
ggplot(rowery, aes(x = factor(education), fill = factor(education))) +
geom_bar() +
labs(title = "Education Distribution", x = "Education Level", y = "Count") +
theme_minimal()
Region: Najwięcej klientów pochodziło z Ameryki Północnej.
ggplot(rowery, aes(x = factor(region), fill = factor(region))) +
geom_bar() +
labs(title = "Region Distribution", x = "Region", y = "Count") +
theme_minimal()
Możemy obliczyć ile kobiet i mężczyzn posiadało dany rodzaj wykształcenia.
table(rowery$gender, rowery$education)
##
## Bachelors Graduate Degree High School Partial College
## Female 146 93 80 132
## Male 160 81 99 133
##
## Partial High School
## Female 38
## Male 38
Czasami chcemy sporządzić rozkład liczebności pewnej zmiennej lub grupy zmiennych, ale tylko dla części respondentów, którzy spełniają określone warunki, na przykład respondentów będących w związku małżeńskim lub respondentów, którzy nie mają dzieci.
W poniższym przykładzie widzimy, że najczęściej rowery kupowały osoby na wysokich stanowiskach (profesjonaliści) oraz wykwalifikowani pracownicy fizyczni.
rower_zawod <- rowery %>%
filter(purchased_bike=="Yes") %>%
count(occupation)
rower_zawod
## # A tibble: 5 × 2
## occupation n
## <chr> <int>
## 1 Clerical 88
## 2 Management 73
## 3 Manual 55
## 4 Professional 150
## 5 Skilled Manual 115
Osoby, które posiadały swój własny dom/mieszkanie były zdecydowanie bardziej skłonne do kupienia roweru.
rower_dom <- rowery %>%
filter(purchased_bike=="Yes") %>%
count(home_owner)
rower_dom
## # A tibble: 2 × 2
## home_owner n
## <chr> <int>
## 1 No 155
## 2 Yes 326
Wyświetlając funkcje ‘view(dfSummary(rowery))’ możemy uzyskać liczebność kategorii dla wszystkich zmiennych jakościowych (zarówno liczbową, jak i procentową).
2. Podstawowe statystyki opisowe
Kolejnym etapem będzie przedstawienie podstawowych statystyk opisowych dla zmiennych ilościowych za pomocą zbiorczej tabelki.
descr(rowery)
## Non-numerical variable(s) ignored: marital_status, gender, education, occupation, home_owner, commute_distance, region, purchased_bike
## Descriptive Statistics
## rowery
## N: 1000
##
## age cars children id income
## ----------------- --------- --------- ---------- ---------- -----------
## Mean 44.12 1.45 1.91 19965.99 55877.61
## Std.Dev 11.15 1.12 1.62 5347.33 29892.86
## Min 25.00 0.00 0.00 11000.00 10000.00
## Q1 35.00 1.00 0.00 15289.50 30000.00
## Median 43.00 1.00 2.00 19744.00 60000.00
## Q3 52.00 2.00 3.00 24475.50 70000.00
## Max 74.00 4.00 5.00 29447.00 130000.00
## MAD 11.86 1.48 1.48 6848.13 29652.00
## IQR 17.00 1.00 3.00 9180.00 40000.00
## CV 0.25 0.77 0.85 0.27 0.53
## Skewness 0.44 0.43 0.39 0.05 0.56
## SE.Skewness 0.08 0.08 0.08 0.08 0.08
## Kurtosis -0.55 -0.39 -1.01 -1.19 -0.15
## N.Valid 1000.00 1000.00 1000.00 1000.00 1000.00
## Pct.Valid 100.00 100.00 100.00 100.00 100.00
Możemy również przedstawić poszczególne statystyki opisowe dla zmiennej ‘income’ w zależności od innych zmiennych np. regionu czy wykonywanego zawodu, co zostało przedstawione poniżej.
X <- rowery %>%
group_by(region) %>%
summarize(Min=min(income),
Max=max(income),
średnia=mean(income),
odchylenie = sd(income),
mediana = median(income),
Q1=quantile(income,0.25),
Q3=quantile(income,0.75))
X
## # A tibble: 3 × 8
## region Min Max średnia odchylenie mediana Q1 Q3
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Europe 10000 130000 39996. 32390. 30000 20000 40000
## 2 North America 10000 130000 62552. 23269. 60000 40000 70000
## 3 Pacific 10000 130000 63034. 32063. 70000 30000 80000
Obserwujemy przykładowo, że średni dochód w Europie jest zdecydowanie niższy niż w Ameryce Północnej czy na Pacyfiku.
Y <- rowery %>%
group_by(occupation) %>%
summarize(Min=min(income),
Max=max(income),
średnia=mean(income),
odchylenie = sd(income),
mediana = median(income),
Q1=quantile(income,0.25),
Q3=quantile(income,0.75))
Y
## # A tibble: 5 × 8
## occupation Min Max średnia odchylenie mediana Q1 Q3
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Clerical 10000 56268. 31165. 7924. 30000 30000 40000
## 2 Management 40000 130000 85723. 26658. 80000 70000 110000
## 3 Manual 10000 56268. 16124. 8953. 10000 10000 20000
## 4 Professional 40000 130000 74153. 21878. 70000 60000 80000
## 5 Skilled Manual 10000 90000 51554. 16734. 50000 40000 60000
Na najwyższe średnie zarobki mogą liczyć pracownicy na stanowiskach kierowniczych.
3. Korelacja
Poniżej za pomocą różnych sposobów przedstawiona została korelacja pomiędzy zmiennymi ilościowymi.
Tabelka korelacji: Największy poziom korelacji możemy zaobserwować pomiędzy wiekiem a liczbą posidanych dzieci (0.53).
cor((rowery[,c(4,5,9,12)]), method="pearson")
## income children cars age
## income 1.0000000 0.2647358 0.4313941 0.1747145
## children 0.2647358 1.0000000 0.2753641 0.5328957
## cars 0.4313941 0.2753641 1.0000000 0.1895145
## age 0.1747145 0.5328957 0.1895145 1.0000000
Korelację możemy zwizualizować również za pomocą następujących wykresów.
corrplot(cor(rowery[,c(4,5,9,12)]), method = "number", type = "upper", diag =FALSE)
corr_matrix<-cor(rowery[,c(4,5,9,12)])
corrplot(corr_matrix, method="color")
IV ETAP- WNIOSKOWANIE STATYSTYCZNE
Testowanie hipotez opiera się na założeniu pewnych warunków w populacji, a następnie analizie próby w celu zweryfikowania, czy dane założenie jest prawdziwe. Statystyki testowe i wartości p-value dostarczają nam narzędzi do dokładnego zrozumienia, czy obserwowane różnice między grupami czy parametrami są statystycznie istotne. Jako poziom istotności przyjęto wartość 0.05.
1. Pytania badawcze
- Czy wartość dochodu jest zależna od regionu?
- Czy rodzaj wykonywanej pracy ma wpływ na fakt posiadania własnościowego domu/mieszkania?
- Czy osoby, które nie posiadają samochodu częściej kupują rowery?
- Czy poziom wykształcenia wpływa na przyszłe dochody?
- Czy w małżeństwach rodzi się więcej dzieci?
1. pytanie badawcze
Dla zmiennych ,region’ i ‘income’ zastosowano polecenie ‘ggbetweenstats’, adekwatne dla porównywania zmiennej jakościowej z ilościową. Przyjęto następujące hipotezy badawcze:
H0: Dochód nie różni się istotnie w zależności od zamieszkiwanego regionu. H1: Dochód istotnie różni się w zależności od zamieszkiwanego regionu.
Na podstawie wartości p-value można stwierdzić, że dochód istotnie różni się w zależności od regionu- należy odrzucić hipotezę zerową.
ggbetweenstats(rowery, region, income)+ggtitle("Średni dochód w zależności od regionu")+labs(x="Region", y="Dochód")
2. pytanie badawcze
Dla zmiennych ‘home_owner’ i ‘occupation’ zastosowano funkcję ‘ggbarstats’, odpowiednią w przypadku porównywania ze sobą dwóch zmiennych jakościowych. Przyjęto następujące hipotezy: H0: Odsetek populacji posiadającej własny dom/mieszkanie nie różni się istotnie w zależności od wykonywanej pracy H1: Odsetek populacji posiadającej własny dom/mieszkanie różni się istotnie w zależności od wykonywanej pracy
Na podstawie wartości p-value można stwierdzić, że odsetek populacji posiadającej własny dom bądź mieszkanie istotnie różni się w zależności od rodzaju wykonywanej pracy-należy przyjąć hipotezę alternatywną.
ggbarstats(rowery, home_owner, occupation)+ggtitle("Posiadanie domu w zależności od rodzaju wykonywanej pracy")+labs(x="Zawód", y="Odsetek", fill="Posiadanie domu")
3. pytanie badawcze
Przyjęto następujące hipotezy badawcze: H0: Średnia liczba posidanych samochodów nie jest istotnie niższa wśród osób posidających rower H1: Średnia liczba posiadanych samochodów jest istotnie niższa wśród osób posidających rower
Na podstawie wartości p-value można stwierdzić, że średnia liczba posiadanych samochodów jest istotnie niższa wśród osób posiadających rower- należy odrzucić hipotezę zerową.
ggbetweenstats(rowery, purchased_bike, cars)+ggtitle("Ilość posiadanych samochodów w zależności od posiadania roweru")+labs(x="Posiadanie roweru", y="Ilość posiadanych samochodów")
4. pytanie badawcze
Przyjęto następujące hipotezy badawcze: H0: Średnie dochody nie różnią się istotnie w zależności od poziomu wykształcenia H1: Średnie dochody różnią się istotnie w zależności od poziomu wykształcenia
Na podstawie otrzymanych wartości p-value można przyjąć, że średnie dochody istotnie różnią się w zależności od poziomu wykształcenia- należy przyjąc hipotezę alternatywną.
ggbetweenstats(rowery, education, income)+ggtitle("Średnie dochody w zależności od poziomu wykształcenia")+labs(x="Wykształcenie", y="Dochody")
5. pytanie badawcze
Przyjęto następujące hipotezy badawcze: H0: Ilość posidanych dzieci nie różni się istotnie w zależności od stanu cywilnego H1: Ilość posiadanych dzieci różni się istotnie w zależności od stanu cywilnego
Na podstawie wartości p-value dotyczącej całości obserwacji można przyjąć, że liczba posiadanych dzieci istotnie różni się w zależności od stanu cywilnego ankietowanej jednostki. Jednak przy analizie obserwacji dotyczącej wyłącznie osób posiadających dwójkę bądź trójkę dzieci, relacji pomiędzy liczbą dzieci, a stanem cywilnym nie można uznać za istotnie zróżnicowaną statystycznie z uwagi na wartość p-value przekraczającą przyjęty poziom istotności wynoszący 0.05.
ggpiestats(rowery, marital_status, children)+ggtitle("Ilość posiadanych dzieci w zależności od stanu cywilnego")
2. Model logitowy
Dla zmiennej ‘purchased_bike’ postanowiono stworzyć model logitowy, który ma za zadanie sprawdzić, jaki wpływ mają poszczególne zmienne na fakt, że klient zdecydował się na kupno roweru. Zmienna ‘purchased_bike’ jest jakościowa i przyjmuje dwie kategorie ‘Yes’ oraz ‘No’. Za pomocą model.matrix zmienna ta została przekodowana tak, aby miała postać binarną, czyli przyjmowała tylko dwie wartości 0 dla ‘No’ i 1 dla ‘Yes’. Zmienne o największej wartości p-value były stopniowo usuwane z początkowego modelu zgodnie z metodą a posteriori, tak, aby otrzymać ostateczny model końcowy- model7.
Poniżej przedstawiona została postać wyjściowa modelu.
encoded_data <- model.matrix(~purchased_bike - 1, data = rowery)
model= glm(encoded_data ~ gender + marital_status + region + cars + income + age + occupation + education + home_owner + commute_distance, family = binomial, data=rowery)
summary(model)
##
## Call:
## glm(formula = encoded_data ~ gender + marital_status + region +
## cars + income + age + occupation + education + home_owner +
## commute_distance, family = binomial, data = rowery)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.703e-02 4.452e-01 0.106 0.915874
## genderMale -1.226e-02 1.374e-01 -0.089 0.928917
## marital_statusSingle -6.247e-01 1.531e-01 -4.080 4.50e-05 ***
## regionNorth America 1.529e-01 2.136e-01 0.716 0.474025
## regionPacific -8.040e-01 2.412e-01 -3.333 0.000860 ***
## cars 4.778e-01 8.914e-02 5.360 8.31e-08 ***
## income -1.360e-05 4.189e-06 -3.246 0.001170 **
## age 5.482e-03 7.140e-03 0.768 0.442646
## occupationManagement 2.600e-01 4.152e-01 0.626 0.531141
## occupationManual 5.634e-02 2.850e-01 0.198 0.843298
## occupationProfessional -2.836e-01 3.296e-01 -0.861 0.389481
## occupationSkilled Manual 1.409e-01 2.668e-01 0.528 0.597499
## educationGraduate Degree 3.562e-01 2.198e-01 1.620 0.105157
## educationHigh School -1.371e-01 2.468e-01 -0.556 0.578415
## educationPartial College 2.646e-01 2.133e-01 1.241 0.214739
## educationPartial High School 5.272e-01 3.584e-01 1.471 0.141218
## home_ownerYes -2.592e-01 1.644e-01 -1.577 0.114883
## commute_distance1-2 Miles 1.307e-01 2.107e-01 0.620 0.535288
## commute_distance10+ Miles 1.059e+00 2.909e-01 3.641 0.000272 ***
## commute_distance2-5 Miles -1.406e-02 2.145e-01 -0.066 0.947737
## commute_distance5-10 Miles 6.088e-01 2.399e-01 2.538 0.011150 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1384.9 on 999 degrees of freedom
## Residual deviance: 1246.9 on 979 degrees of freedom
## AIC: 1288.9
##
## Number of Fisher Scoring iterations: 4
Kolejno zgodnie z metodą a posteriori usunięte zostały zmienne: commute_distance, gender, occupation, education, home_owner (w wymienionej kolejności). Poniżej przedstawiona została postać końcowa modelu po eliminacji wszystkich nieistotnych statystycznie zmiennych.
model7= glm(encoded_data ~ marital_status + cars + income + age, family = binomial, data=rowery)
summary(model7)
##
## Call:
## glm(formula = encoded_data ~ marital_status + cars + income +
## age, family = binomial, data = rowery)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.274e-01 3.033e-01 -0.750 0.453550
## marital_statusSingle -4.999e-01 1.367e-01 -3.657 0.000255 ***
## cars 5.112e-01 6.888e-02 7.422 1.16e-13 ***
## income -1.300e-05 2.547e-06 -5.107 3.28e-07 ***
## age 1.182e-02 6.205e-03 1.906 0.056694 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1384.9 on 999 degrees of freedom
## Residual deviance: 1302.4 on 995 degrees of freedom
## AIC: 1312.4
##
## Number of Fisher Scoring iterations: 4
Otrzymane współczynniki modelu nie nadają się do interpretacji, gdyż za ich pomocą możemy określić jedynie w jakim stopniu zmieniał się logit. Liczymy więc ilorazy szans, aby otrzymać wartości do interpretacji.
OR=exp(model7$coefficients)
OR
## (Intercept) marital_statusSingle cars
## 0.7966333 0.6065888 1.6673378
## income age
## 0.9999870 1.0118945
INTERPRETACJE: Jeżeli klient sklepu rowerowego jest singlem (w porównaniu do osób posiadających żonę/męża), to szanse na to, że kupi on rower maleją o około 40%. Z każdym posiadanym przez klienta samochodem, szanse na to, że kupi on rower wzrastają o około 67%. Wraz ze wzrostem wieku o 1 rok, szanse na to, że klient kupi rower wzrastają o około 1,19%. Wraz ze wzrostem dochodu, szanse na to, że klient kupi rower praktycznie się nie zmienniają. Wszelkie opisane zależności następują przy założeniu niezmienności pozostałych czynników- ceteris paribus.
Możemy policzyć również efekty krańcowe, aby dowiedzieć się, jak poszczególne zmienne objaśniające wspływają na prawdopodobieństwo zakupu roweru przez klienta sklepu. Prawdopodobieństwo zsotało wyliczone przy założeniu średniego poziomu pozostałych zmiennych objaśniających.
logitmfx(encoded_data ~ marital_status + cars + income + age, data= rowery)
## Call:
## logitmfx(formula = encoded_data ~ marital_status + cars + income +
## age, data = rowery)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## marital_statusSingle -1.2422e-01 3.3639e-02 -3.6928 0.0002218 ***
## cars 1.2760e-01 1.7195e-02 7.4209 1.163e-13 ***
## income -3.2460e-06 6.3571e-07 -5.1060 3.290e-07 ***
## age 2.9513e-03 1.5486e-03 1.9058 0.0566802 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "marital_statusSingle"
Możemy również sprawdzić stopień dopasowania naszego modelu i stowrzyć wykres obrazujący punkty odcięcia oraz krzywą AUC.
pred=prediction(model7$fitted.values, rowery$purchased_bike)
pref=performance(pred, 'tpr', 'fpr')
plot(pref, colorize= TRUE)
performance(pred, 'auc')@y.values
## [[1]]
## [1] 0.3357208
AUC przyjmuje jednak dosyć niską wartość.
PODSUMOWANIE
Niniejszy projekt przebiegał w kilku kluczowych etapach. Początkowo nastąpiło przygotowanie danych do analizy. Dokonano imputacji 53 braków danych oraz przekształcono za pomocą metody capping wartości odstające dla zmiennych ‘income’ oraz ‘age’. To tak naprawdę późniejsze etapy projektu stanowią źródło wiedzy i przedstawiają pewne wnioski, które mogą okazać się cenne z perspektywy kierownictwa sklepu rowerowego, ale także dla szerokiej gamy interesariuszy, zarówno w branży rowerowej, jak i poza nią z kilku powodów:
Skoro dochód klientów istotnie różni się w zależności od regionu, to strategie cenowe i marketingowe sklepu powinny być dostosowane regionalnie. Sklep może oferować różne gamy produktów w różnych regionach, uwzględniając siłę nabywczą klientów.
Fakt, że osoby posiadające rower mają średnio mniej samochodów, może wskazywać na ekologiczną świadomość lub ograniczenia finansowe. Sklep może więc promować rowery jako ekonomiczną i ekologiczną alternatywę transportu.
Różnice w liczbie posiadanych dzieci w zależności od stanu cywilnego mogą być użyteczne przy tworzeniu ofert skierowanych do rodzin. Na przykład, sklep może oferować promocje na rowery dziecięce lub akcesoria dla rodzin. Analiza tej zależności może pomóc w decydowaniu, gdzie najlepiej otworzyć nowe sklepy lub zwiększyć działania marketingowe, kierując je na obszary o większej liczbie rodzin z dziećmi, które mogą być potencjalnymi klientami.