##### Rozdział 13: trudne dane --------------------

## Regresja krokowa ----

# wczytujemy dane i przeprowadzamy proste przygotowania
library(tidyverse)
titanic_train <- read_csv("titanic_train.csv") |>
  mutate(
    Age_MVI = if_else(is.na(Age), 1, 0),
    Age = if_else(is.na(Age), mean(Age, na.rm = TRUE), Age),
    Cabin = if_else(is.na(Cabin), "X", Cabin),
    Embarked = factor(if_else(is.na(Embarked), "X", Embarked)),
    Sex = factor(Sex)
  )

# określamy najprostszy model regresji logistycznej
simple_model <- glm(Survived ~ 1, family = binomial, data = titanic_train)
  
# określamy pełny model regresji logistycznej
full_model <- glm(Survived ~ Age + Age_MVI + Embarked + Sex + Pclass + SibSp + Fare,
                  family = binomial, data = titanic_train)

# regresja krokowa z selekcją w przód
sw_forward <- stats::step(simple_model, scope = formula(full_model),
                          direction = "forward")

# uzyskujemy formułę końcowego modelu
formula(sw_forward)

# współczynniki regresji końcowego modelu
sw_forward$coefficients

# eliminacja wsteczna
sw_backward <- stats::step(full_model, direction = "backward")

## Selekcja cech za pomocą pakietu Boruta ----

set.seed(12345) # gwarantuje powtarzalność wyników
# tworzymy cechę o losowych wartościach, aby zademonstrować bezużyteczną cechę
titanic_train$rand_vals <- runif(n = 891, min = 1, max = 100)

# wykonujemy algorytm Boruta na zbiorze danych o pasażerach Titanica
# (może to trwać dość długo w przypadku większych zbiorów danych)
library(Boruta)
titanic_boruta <- Boruta(Survived ~ PassengerId + Age + 
                           Sex + Pclass + SibSp + rand_vals,
                         data = titanic_train, doTrace = 1)
# sprawdzamy wynik
titanic_boruta

# kreślimy ważność cech
plot(titanic_boruta)

## Analiza składowych głównych (PCA) ----

library(tidyverse) # wczytujemy zestaw pakietów tidyverse

sns_data <- read_csv("snsdata.csv") # wczytujemy dane o nastoletnich użytkownikach mediów społecznościowych

# wybieramy tylko 36 kolumn, od kolumny o nazwie 'basketball' do kolumny o nazwie 'drugs'
# każda kolumna wskazuje, ile razy użyto danego wyrazu na profilu społecznościowym
sns_terms <- sns_data |> select(basketball:drugs)

# biblioteka irlba zawiera efektywniejszą funkcję PCA, niż prcomp() wbudowana w R
library(irlba)

# przeprowadzamy PCA - zauważ, że środkujemy i przeskalowujemy dane
set.seed(2023) # gwarantuje zgodność wyników z książką
sns_pca <- sns_terms |> 
  prcomp_irlba(n = 10, center = TRUE, scale = TRUE) # znajdujemy 10 pierwszych składowych głównych w danych SNS

# tworzymy wykres piargowy analizy PCA danych SNS
screeplot(sns_pca, npcs = 10, type = "lines",
          main = "Wykres piargowy składowych głównych w danych SNS")

# używamy funkcji summary, aby zobaczyć składowe i proporcję wyjaśnianej zmienności
summary(sns_pca)

# badamy obiekt PCA -- najbardziej interesują nas komponenty $x i $rotation 
str(sns_pca)

# komponent $x to przekształcona wersja pierwotnych danych
str(sns_pca$x)

# $x to nasze pierwotne dane przekształcone tak, aby miały nowe "cechy" -- składowe główne
nrow(sns_pca$x) # powinno być 30 000 wierszy
head(sns_pca$x) # powinno być 10 kolumn

# tworzymy "długą" wersję zbioru danych PCA do celów wizualizacji
sns_pca_long <- tibble(SNS_Term = colnames(sns_terms), as_tibble(sns_pca$rotation)) |> # dodajemy etykiety wierszy
  pivot_longer(PC1:PC10, names_to = "PC", values_to = "Contribution") # przechodzimy od szerokiego do długiego zbioru danych

# używamy ggplot, aby zwizualizować wyrazy ważne dla PC3
library(ggplot2)

sns_pca_long |>
  filter(PC == "PC3") |>
  top_n(15, abs(Contribution)) |>
  mutate(SNS_Term = reorder(SNS_Term, Contribution)) |>
  ggplot(aes(SNS_Term, Contribution, fill = SNS_Term)) +
  geom_col(show.legend = FALSE, alpha = 0.8) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1,
                                   vjust = 0.5), axis.ticks.x = element_blank()) +
  labs(x = "Słowo z mediów społecznościowych",
       y = "Wzgledna ważność dla składowej głównej",
       title = "15 słów o największym wkładzie w PC3")

# tworzymy funkcję do wizualizacji pozostałych czterech składowych
plot_pca <- function(component) {
  sns_pca_long |>
    filter(PC == component) |>
    top_n(15, abs(Contribution)) |>
    mutate(SNS_Term = reorder(SNS_Term, Contribution)) |>
    ggplot(aes(SNS_Term, Contribution, fill = SNS_Term)) +
    geom_col(show.legend = FALSE, alpha = 0.8) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
          axis.ticks.x = element_blank()) + 
    labs(x = "Słowo z mediów społecznościowych",
         y = "Wzgledna ważność dla składowej głównej",
         title = paste("15 słów o największym wkładzie w", component))
}

# używamy funkcji
plot_pca("PC1")
plot_pca("PC2")
plot_pca("PC4")
plot_pca("PC5")

# możemy użyć składowych głównych do przewidywania liczby znajomych
sns_data_pca <- cbind(sns_data[1:4], sns_pca$x) # łączymy składowe główne z pierwotnymi danymi

# tworzymy model regresji liniowej przewidujący liczbę znajomych na podstawie składowych głównych
m <- lm(friends ~ PC1 + PC2 + PC3 + PC4 + PC5, data = sns_data_pca)

m # wyświetlamy współczynniki modelu

## Zmiana odwzorowania rozrzedzonych danych kategorycznych ----

library(tidyverse)

# wczytujemy zbiór danych o pasażerach Titanica i tworzymy cechę Title (z rozdziału 12.)
titanic_train <- read_csv("titanic_train.csv") |>
  mutate(Title = str_extract(Name, ", [A-z]+\\.")) |>
  mutate(Title = str_replace_all(Title, "[, \\.]", ""))

# cecha Title ma dużą liczbę kategorii
table(titanic_train$Title, useNA = "ifany")

# grupujemy kategorie o podobnym znaczeniu
titanic_train <- titanic_train |>
  mutate(TitleGroup = fct_collapse(Title, 
    Mr = "Mr",
    Mrs = "Mrs",
    Master = "Master",
    Miss = c("Miss", "Mlle", "Mme", "Ms"),
    Noble = c("Don", "Sir", "Jonkheer", "Lady"),
    Military = c("Capt", "Col", "Major"),
    Doctor = "Dr",
    Clergy = "Rev",
    other_level = "Other")
  ) |>
  mutate(TitleGroup = fct_na_value_to_level(TitleGroup,
                                            level = "Unknown"))

# badamy przekodowane dane
table(titanic_train$TitleGroup)

# przyglądamy się liczebnościom i proporcjom wszystkich poziomów, posortowanych od największego do najmniejszego
look at the counts and proportions of all levels, sorted largest to smallest
fct_count(titanic_train$Title, sort = TRUE, prop = TRUE)

# zlepiamy wszystko poza trzema najliczniejszymi poziomami
table(fct_lump_n(titanic_train$Title, n = 3))

# zlepiamy wszystko o proporcjach mniejszych niż 1%
table(fct_lump_prop(titanic_train$Title, prop = 0.01))

# zlepiamy wszystko, co ma mniej niż 5 obserwacji
table(fct_lump_min(titanic_train$Title, min = 5))

## Dzielenie rozrzedzonych danych liczbowych na przedziały ----

# badamy opłaty za podróż Titanikiem
head(titanic_train$Fare)
summary(titanic_train$Fare)

# tworzymy zmienną binarną dla pierwszej/drugiej klasy
titanic_train <- titanic_train |> mutate(
  fare_firstclass = if_else(Fare >= 31, 1, 0, missing = 0)
)

# tabelaryzujemy wartości binarne
table(titanic_train$fare_firstclass)

# tworzymy trójpoziomową cechę za pomocą case_when()
titanic_train <- titanic_train |>
  mutate(
    fare_class = case_when(
      Fare >= 31 ~ "1st Class",
      Fare >= 15 ~ "2nd Class",
      TRUE ~ "3rd Class"
    )
  )

# badamy wynik
table(titanic_train$fare_class)

# funkcja cut() może zrobić to samo, co powyższa instrukcja case_when()
table(cut(titanic_train$Fare, breaks = c(-Inf, 15, 31, Inf),
          right = FALSE))

# używamy cut() w połączeniu z seq(), aby wygenerować interwały równej szerokości
table(cut(titanic_train$Fare, right = FALSE,
          breaks = seq(from = 0, to = 550, by = 50)))

# używamy cut() w połączeniu z quantiles() i seq(), aby utworzyć przedziały o równej liczbie przykładów
table(cut(titanic_train$Fare, right = FALSE,
          breaks = quantile(titanic_train$Fare,
                            probs = seq(0, 1, 0.20))))

# używamy funkcji ntile() z tidyverse do utworzenia pięciu przedziałów
table(ntile(titanic_train$Fare, n = 5))

# przekształcamy grupy ntile() w czynnik
titanic_train <- titanic_train |>
  mutate(fare_level = factor(ntile(Fare, n = 11)))

table(titanic_train$fare_level)

## Imputacja brakujących wartości ----

library(readr)
titanic_train <- read_csv("titanic_train.csv")

# zastępujemy brakujące dane kategoryczne dowolnym łańcuchem tekstu
titanic_train <- titanic_train |>
  mutate(
    Cabin = if_else(is.na(Cabin), "X", Cabin),
    Embarked = if_else(is.na(Embarked), "Unknown", Embarked)
  )

# zastępujemy brakujące wartości wieku wartością średnią i tworzymy wskaźnik braku wartości
impute mean value and create missing value indicator for age
titanic_train <- titanic_train |>
  mutate(
    Age_MVI = if_else(is.na(Age), 1, 0),
    Age = if_else(is.na(Age), mean(Age, na.rm = TRUE), Age)
  )

## Proste strategie przywracania równowagi danych ----

# wczytujemy i przygotowujemy zbiór danych o nastoletnich użytkownikach mediów społecznościowych
library(tidyverse)

snsdata <- read_csv("snsdata.csv") |>
  mutate(
    gender = fct_recode(gender, Female = "F", Male = "M"),
    gender = fct_na_value_to_level(gender, level = "Unknown"),
    age = ifelse(age < 13 | age > 20, NA, age) # zastępujemy odstające wartości wieku
  ) |>
  group_by(gradyear) |>
  mutate(age_imp = if_else(is.na(age), median(age, na.rm = TRUE), age)) |>
  ungroup() |>
  select(gender, friends, gradyear, age_imp, basketball:drugs)

# badamy początkową nierównowagę klas
fct_count(snsdata$gender, prop = TRUE)

# podpróbkowujemy klasy większościowe
library(caret)
sns_undersample <- downSample(x = snsdata[2:40], y = snsdata$gender, yname = "gender")
fct_count(sns_undersample$gender, prop = TRUE)

# nadpróbkowujemy klasy mniejszościowe
library(caret)
sns_oversample <- upSample(x = snsdata[2:40], y = snsdata$gender, yname = "gender")
fct_count(sns_oversample$gender, prop = TRUE)

## Generowanie syntetycznego zrównoważonego zbioru danych z wykorzystaniem algorytmu SMOTE ----

# tworzymy zbiór danych zrównoważony pod względem płci za pomocą SMOTE
library(themis)
sns_balanced <- snsdata |> smote("gender") # prosta składnia (bez normalizacji)

# sprawdzamy, czy zbiór danych jest teraz zrównoważony pod względem płci
table(sns_balanced$gender)

# aby uzyskać lepsze wyniki SMOTE, definiujemy funkcję normalizacji (wprowadzoną w rozdziale 3.)
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

# funkcja unnormalize() (wprowadzona w rozdziale 7.) przywraca pierwotną skalę danych 
unnormalize <- function(norm_vals, col_name) {
  old_vals <- snsdata[col_name]
  unnormalized_vals <- norm_vals * (max(old_vals) - min(old_vals)) + min(old_vals)
  
  # zaokrąglamy wszystkie kolumny z wyjątkiem age_imp do wartości całkowitych
  rounded_vals <- if(col_name != "age_imp") { round(unnormalized_vals) }
                  else {unnormalized_vals}
  
  return (rounded_vals)
}

# bardziej zaawansowany proces smote() ze znormalizowanymi danymi
snsdata_balanced <- snsdata |>
  mutate(across(where(is.numeric), normalize)) |> # normalizujemy dane liczbowe
  smote("gender") |>
  mutate(across(where(is.numeric), ~unnormalize(.x, cur_column()))) # odwracamy normalizację danych

# sprawdzamy, czy równoważenie zbioru danych zadziałało poprawnie
table(snsdata$gender)
table(snsdata_balanced$gender)
