ANALIZA DANYCH-PROJEKT

Maciej Fleks, Małgorzata Dudanowicz, Krzysztof Kowalski

2024-01-25

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.