##### Rozdział 6: metody regresji -------------------

#### Część 1: regresja liniowa -------------------

## Regresja ----
## Przykład: dane dotyczące startów wahadłowca ----
launch <- read.csv("challenger.csv")

# ręczne oszacowanie współczynnika beta
b <- cov(launch$temperature, launch$distress_ct) / var(launch$temperature)
b

# ręczne oszacowanie współczynnika alfa
a <- mean(launch$distress_ct) - b * mean(launch$temperature)
a

# obliczamy korelację danych
r <- cov(launch$temperature, launch$distress_ct) /
       (sd(launch$temperature) * sd(launch$distress_ct))
r
cor(launch$temperature, launch$distress_ct)

# tworzymy prostą funkcję regresji wielokrotnej
reg <- function(y, x) {
  x <- as.matrix(x)
  x <- cbind(Intercept = 1, x)
  b <- solve(t(x) %*% x) %*% t(x) %*% y
  colnames(b) <- "estimate"
  print(b)
}

# badamy dane dotyczące startów
str(launch)

# testujemy model regresji z prostą regresją liniową
reg(y = launch$distress_ct, x = launch[2])

# używamy modelu regresji z regresją wieloraką
reg(y = launch$distress_ct, x = launch[2:4])

# możesz potwierdzić, że niestandardowa funkcja regresji wielorakiej działa
# poprawnie, przez porównanie jej z wbudowaną funkcją lm języka R
model <- lm(distress_ct ~ temperature + field_check_pressure + flight_num, data = launch)
model

## Przykład: przewidywanie kosztów likwidacji szkód ----
## Etap 2: badanie i przygotowywanie danych ----
insurance <- read.csv("autoinsurance.csv", stringsAsFactors = TRUE)
str(insurance)

# podsumowujemy zmienną wydatków
summary(insurance$expenses)

# histogram wydatków
hist(insurance$expenses)

# tabele z cechami kategorycznymi
table(insurance$geo_area)
table(insurance$vehicle_type)

# badamy związki między cechami: macierz korelacji
cor(insurance[c("age", "est_value", "miles_driven", "expenses")])

# wizualizujemy związki między cechami: macierz wykresów punktowych
pairs(insurance[c("age", "est_value", "miles_driven",
                  "expenses")], pch = ".")

# użyteczniejsza macierz wykresów punktowych
library(psych)
pairs.panels(insurance[c("age", "est_value", "miles_driven",
                         "expenses")], pch = ".")

## Etap 3: trenowanie modelu na danych ----
ins_model <- lm(expenses ~ age + geo_area + vehicle_type +
                  est_value + miles_driven +
                  college_grad_ind + speeding_ticket_ind +
                  hard_braking_ind + late_driving_ind +
                  clean_driving_ind,
                data = insurance)

ins_model <- lm(expenses ~ ., data = insurance) # równoważne powyższemu poleceniu

# wypisujemy szacowane współczynniki beta
options(scipen = 999) # wyłączamy notację naukową
ins_model

## Etap 4: ewaluacja modelu ----
# wypisujemy bardziej szczegółowe informacje o szacowanych współczynnikach beta
summary(ins_model)

## Etap 5: poprawianie działania modelu ----

# dodajemy wyraz "age" wyższego rzędu
insurance$age2 <- insurance$age^2

# tworzymy ostateczny model
ins_model2 <- lm(expenses ~ . + hard_braking_ind:late_driving_ind,
                 data = insurance)

summary(ins_model2)

# dokonujemy prognoz za pomocą modelu regresji
insurance$pred <- predict(ins_model2, insurance)
cor(insurance$pred, insurance$expenses)

plot(insurance$pred, insurance$expenses)
abline(a = 0, b = 1, col = "red", lwd = 3, lty = 2)

predict(ins_model2,
        data.frame(age = 30, age2 = 30^2, geo_area = "rural", 
                   vehicle_type = "truck", est_value = 25000,
                   miles_driven = 14000, college_grad_ind = 0,
                   speeding_ticket_ind = 0, hard_braking_ind = 0,
                   late_driving_ind = 0, clean_driving_ind = 1))

predict(ins_model2,
        data.frame(age = 30, age2 = 30^2, geo_area = "rural", 
                   vehicle_type = "truck", est_value = 25000,
                   miles_driven = 14000, college_grad_ind = 0,
                   speeding_ticket_ind = 0, hard_braking_ind = 0,
                   late_driving_ind = 0, clean_driving_ind = 0))

predict(ins_model2,
        data.frame(age = 30, age2 = 30^2, geo_area = "rural", 
                   vehicle_type = "truck", est_value = 25000,
                   miles_driven = 24000, college_grad_ind = 0,
                   speeding_ticket_ind = 0, hard_braking_ind = 0,
                   late_driving_ind = 0, clean_driving_ind = 0))

## Przewidywanie odpływu ubezpieczonych z wykorzystaniem regresji logistycznej ----

churn_data <- read.csv("insurance_churn.csv")

prop.table(table(churn_data$churn)) # wypisujemy procentowy odpływ

# tworzymy model regresji logistycznej
churn_model <- glm(churn ~ . -member_id, data = churn_data,
                   family = binomial(link = "logit"))

# badamy estymatory modelu regresji logistycznej
summary(churn_model)

# wczytujemy testowy zbiór danych
churn_test <- read.csv("insurance_churn_test.csv")

# dokonujemy prognoz na zbiorze testowym
churn_test$churn_prob <- predict(churn_model, churn_test,
                                 type = "response")

# badamy przewidziane wartości
summary(churn_test$churn_prob)

# wypisujemy ubezpieczonych o największym prawdopodobieństwie odpływu
churn_order <- order(churn_test$churn_prob, decreasing = TRUE)
head(churn_test[churn_order, c("member_id", "churn_prob")], n = 5)

#### Część 2: drzewa regresji i drzewa modeli -------------------

## Drzewa regresji i drzewa modeli ----
## Przykład: obliczanie SDR ----
# przygotowujemy dane
tee <- c(1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7)
at1 <- c(1, 1, 1, 2, 2, 3, 4, 5, 5)
at2 <- c(6, 6, 7, 7, 7, 7)
bt1 <- c(1, 1, 1, 2, 2, 3, 4)
bt2 <- c(5, 5, 6, 6, 7, 7, 7, 7)

# obliczamy SDR
sdr_a <- sd(tee) - (length(at1) / length(tee) * sd(at1) + length(at2) / length(tee) * sd(at2))
sdr_b <- sd(tee) - (length(bt1) / length(tee) * sd(bt1) + length(bt2) / length(tee) * sd(bt2))

# porównujemy SDR dla każdego podziału
sdr_a
sdr_b

## Przykład: ocenianie jakości win ----
## Etap 2: badanie i przygotowywanie danych ----
wine <- read.csv("whitewines.csv")

# badanie danych o winach
str(wine)

# rozkład jakości win
hist(wine$quality)

# zbiorcze statystyki danych
summary(wine)

wine_train <- wine[1:3750, ]
wine_test <- wine[3751:4898, ]

## Etap 3: trenowanie modelu na danych ----
# budujemy drzewo regresji za pomocą rpart
library(rpart)
m.rpart <- rpart(quality ~ ., data = wine_train)

# wypisujemy podstawowe informacje o drzewie
m.rpart

# wypisujemy bardziej szczegółowe informacje o drzewie
summary(m.rpart)

# używamy pakietu rpart.plot do utworzenia wizualizacji
library(rpart.plot)

# podstawowy diagram drzewa decyzyjnego
rpart.plot(m.rpart, digits = 3)

# kilka poprawek diagramu
rpart.plot(m.rpart, digits = 4, fallen.leaves = TRUE, type = 3, extra = 101)

## Etap 4: ewaluacja modelu ----

# generujemy prognozy na testowym zbiorze danych
p.rpart <- predict(m.rpart, wine_test)

# porównujemy rozkład wartości przewidzianych i rzeczywistych
summary(p.rpart)
summary(wine_test$quality)

# porównujemy korelację
cor(p.rpart, wine_test$quality)

# funkcja do obliczania średniego błędu bezwzględnego
MAE <- function(actual, predicted) {
  mean(abs(actual - predicted))  
}

# średni błąd bezwzględny między wartościami przewidzianymi i rzeczywistymi
MAE(p.rpart, wine_test$quality)

# średni błąd bezwzględny między wartościami rzeczywistymi a wartością średnią
mean(wine_train$quality) # wynik = 5.87
MAE(5.87, wine_test$quality)

## Etap 5: poprawianie działania modelu ----
# trenujemy drzewo modeli Cubist
library(Cubist)
m.cubist <- cubist(x = wine_train[-12], y = wine_train$quality)

# wypisujemy podstawowe informacje o drzewie modeli
m.cubist

# wyświetlamy samo drzewo
summary(m.cubist)

# generujemy prognozy modelu
p.cubist <- predict(m.cubist, wine_test)

# statystyki zbiorcze prognoz
summary(p.cubist)

# korelacja między wartościami przewidzianymi a rzeczywistymi
cor(p.cubist, wine_test$quality)

# średni błąd bezwzględny między wartościami przewidzianymi a rzeczywistymi
# (używa zdefiniowanej wyżej niestandardowej funkcji)
MAE(wine_test$quality, p.cubist)
