## Statystyka praktyczna w data science (R)
## Rozdział 1. Badania eksploracyjne
# > (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck

# Importuje wymagane pakiety.

library(dplyr)
library(tidyr)
library(ggplot2)
library(vioplot)
library(corrplot)
library(gmodels)
library(matrixStats)

# Importuje zestawy danych używane w rozdziale 1.

PSDS_PATH <- file.path(dirname(dirname(getwd())))
 
state <- read.csv(file.path(PSDS_PATH, 'dane', 'state.csv'))
dfw <- read.csv(file.path(PSDS_PATH, 'dane', 'dfw_airline.csv'))
sp500_px <- read.csv(file.path(PSDS_PATH, 'dane', 'sp500_data.csv.gz'), row.names=1)
sp500_sym <- read.csv(file.path(PSDS_PATH, 'dane', 'sp500_sectors.csv'), stringsAsFactors = FALSE)
kc_tax <- read.csv(file.path(PSDS_PATH, 'dane', 'kc_tax.csv.gz'))
lc_loans <- read.csv(file.path(PSDS_PATH, 'dane', 'lc_loans.csv'))
airline_stats <- read.csv(file.path(PSDS_PATH, 'dane', 'airline_stats.csv'), stringsAsFactors = FALSE)
airline_stats$airline <- ordered(airline_stats$airline, 
                                 levels=c('Alaska', 'American', 'Jet Blue', 'Delta', 'United', 'Southwest'))

## Miary położenia
### Przykład: miara położenia dla wielkości populacji i wskaźnika morderstw

# Tabela 1.2
state_asc <- state
state_asc[['Population']] <- formatC(state_asc[['Population']], format='d', digits=0, big.mark=',')
state_asc[1:8,]

mean(state[['Population']])
mean(state[['Population']], trim=0.1)
median(state[['Population']])

weighted.mean(state[['Murder.Rate']], w=state[['Population']])
library('matrixStats')
weightedMedian(state[['Murder.Rate']], w=state[['Population']])

## Miary rozproszenia

sd(state[['Population']])
IQR(state[['Population']])
mad(state[['Population']])

### Percentyle i boxploty

quantile(state[['Murder.Rate']], p=c(.05, .25, .5, .75, .95))

boxplot(state[['Population']]/1000000, ylab='Populacja (w milionach)')

### Tablica częstości i histogramy

breaks <- seq(from=min(state[['Population']]),
              to=max(state[['Population']]), length=11)
pop_freq <- cut(state[['Population']], breaks=breaks,
                right=TRUE, include.lowest=TRUE)
state['PopFreq'] <- pop_freq
table(pop_freq)

options(scipen=5)
hist(state[['Population']], breaks=breaks)

### Szacowanie gęstości
# Gęstość stanowi rozwiązanie alternatywne dla histogramów; dzięki niej widać wyraźniej rozkład punktów danych.

hist(state[['Murder.Rate']], freq=FALSE )
lines(density(state[['Murder.Rate']]), lwd=3, col='blue')

## Badanie danych binarnych i skategoryzowanych

barplot(as.matrix(dfw) / 6, cex.axis=0.8, cex.names=0.7,
        xlab='Powód opóźnienia', ylab='Liczba')

## Korelacja

telecom <- sp500_px[, sp500_sym[sp500_sym$sector == 'telecommunications_services', 'symbol']]
telecom <- telecom[row.names(telecom) > '2012-07-01',]
telecom_cor <- cor(telecom)

# Następnie skoncentrujmy się na zwrotach dla głównych funduszy (sector == 'etf').

etfs <- sp500_px[row.names(sp500_px) > '2012-07-01',
                 sp500_sym[sp500_sym$sector == 'etf', 'symbol']]
corrplot(cor(etfs), method='ellipse')

### Wykres punktowy

# plot(telecom$T, telecom$VZ, xlab='T', ylab='VZ', cex=.8)
plot(telecom$T, telecom$VZ, xlab='ATT (T)', ylab='Verizon (VZ)')
abline(h=0, v=0, col='grey')
dim(telecom)

## Badanie dwóch lub więcej zmiennych
# Wczytujemy zestaw danych kc_tax i filtrujemy go na podstawie zróżnicowania kryteriów:

kc_tax0 <- subset(kc_tax, TaxAssessedValue < 750000 &
                  SqFtTotLiving > 100 &
                  SqFtTotLiving < 3500)
nrow(kc_tax0)

### Wykres hexagonal binning i kontury
#### Przedstawianie danych numerycznych względem danych numerycznych

# W przypadku dużej liczby punktów danych wykres punktowy przestaje być czytelny. Przydatniejsze stają się metody wizualizujące gęstość. Jednym z lepszych rozwiązań jest metoda `stat_hexbin`.

graph <- ggplot(kc_tax0, (aes(x=SqFtTotLiving, y=TaxAssessedValue))) +
  stat_binhex(color='white') +
  theme_bw() +
  scale_fill_gradient(low='white', high='blue') +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
  labs(x='Powierzchnia domu', y='Wartość podatku')
graph

# Wizualizacja dwywymiarowego rozwinięcia wykresu gęstości.

graph <- ggplot(kc_tax0, aes(SqFtTotLiving, TaxAssessedValue)) +
  theme_bw() +
  geom_point(color='blue', alpha=0.1) +
  geom_density2d(color='white') +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
  labs(x='Powierzchnia domu', y='Wartość podatku')
graph

### Dwie zmienne skategoryzowane
# Wczytuje zestaw danych `lc_loans`:

x_tab <- CrossTable(lc_loans$grade, lc_loans$status,
                    prop.c=FALSE, prop.chisq=FALSE, prop.t=FALSE)

### Dane kategoryzowane i numeryczne
# Boxploty kolumnowe mogą być grupowane na podstawie innej kolumny.

boxplot(pct_carrier_delay ~ airline, data=airline_stats, ylim=c(0, 50),
        cex.axis=.6, ylab='Procent opóźnionych lotów')

# Odmiana boxplotów o nazwie _wykres skrzypcowy_. 

graph <- ggplot(data=airline_stats, aes(airline, pct_carrier_delay)) +
  ylim(0, 50) +
  geom_violin(draw_quantiles = c(.25,.5,.75), linetype=2) +
  geom_violin(fill=NA, size=1.1) +
  llabs(x='', y='Procent opóźnionych lotów') +
  theme_bw()
graph

### Wizualizacja wielu zmiennych

graph <- ggplot(subset(kc_tax0, ZipCode %in% c(98188, 98105, 98108, 98126)),
                aes(x=SqFtTotLiving, y=TaxAssessedValue)) +
  stat_binhex(colour='white') +
  theme_bw() +
  scale_fill_gradient(low='gray95', high='black') +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
  labs(x='Powierzchnia domu', y='Wartość podatku') +
  facet_wrap('ZipCode')
graph
