## Statystyka praktyczna w data science (R)
## Rozdział 3. Eksperymenty statystyczne i testowanie istotności
# > (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck

# Importuje wymagane pakiety R.

# Pakiet lmPerm jest dostępny w repozytorium CRAN
# conda install r-devtools
# devtools::install_github('mtorchiano/lmPerm')

library(ggplot2)
library(dplyr)
library(lmPerm)
library(pwr)

# Definiuje ścieżki do zestawów danych. Jeżeli przechowujesz dane w innym katalogu, wprowadź poniżej stosowne zmiany.

PSDS_PATH <- file.path(dirname(dirname(getwd())))


session_times <- read.csv(file.path(PSDS_PATH, 'dane', 'web_page_data.csv'))
session_times[,2] <- session_times[,2] * 100
four_sessions  <- read.csv(file.path(PSDS_PATH, 'dane', 'four_sessions.csv'))
click_rate <-  read.csv(file.path(PSDS_PATH, 'dane', 'click_rates.csv'))
imanishi <-  read.csv(file.path(PSDS_PATH, 'dane', 'imanishi_data.csv'))

## Test randomizacyjny
### Przykład: licznik odwiedzin strony

graph <- ggplot(session_times, aes(x=Page, y=Time)) +
  geom_boxplot() +
  labs(y='Czas (w sekundach') +
  theme_bw()
graph

mean_a <- mean(session_times[session_times['Page'] == 'Page A', 'Time'])
mean_b <- mean(session_times[session_times['Page'] == 'Page B', 'Time'])
mean_b - mean_a

## Przykład testu permutacyjnego
perm_fun <- function(x, nA, nB)
{
  n <- nA + nB
  idx_b <- sample(1:n, nB)
  idx_a <- setdiff(1:n, idx_b)
  mean_diff <- mean(x[idx_b]) - mean(x[idx_a])
  return(mean_diff)
}

perm_diffs <- rep(0, 1000)
for (i in 1:1000) {
  perm_diffs[i] = perm_fun(session_times[, 'Time'], 21, 15)
}
par(mar=c(4,4,1,0)+.1)
hist(perm_diffs, xlab='Różnice w czasach trwania sesji (w sekundach)', main='')
abline(v=mean_b - mean_a, lty=2, lwd=1.5)
text('  Zaobserwowana\n różnica', x=mean_b - mean_a,  y=par()$usr[4]-20, adj=0)

mean(perm_diffs > (mean_b - mean_a))

## Istotność statystyczna i p-wartości

obs_pct_diff <- 100 * (200 / 23739 - 182 / 22588)
conversion <- c(rep(0, 45945), rep(1, 382))
perm_diffs <- rep(0, 1000)
for (i in 1:1000) {
  perm_diffs[i] = 100 * perm_fun(conversion, 23739, 22588)
}

hist(perm_diffs, xlab='Poziom zwrotu (w procentach)', main='')
abline(v=obs_pct_diff, lty=2, lwd=1.5)
text('   Zaobserwowana\n   różnica', x=obs_pct_diff,  y=par()$usr[4]-20, adj=0)

# Zaobserwowana różnica
obs_pct_diff

### p-wartość

mean(perm_diffs > obs_pct_diff)

prop.test(x=c(200,182), n=c(23739,22588), alternative='greater')

## Test t

t.test(Time ~ Page, data=session_times, alternative='less')

## ANOVA

graph <- ggplot(four_sessions, aes(x=Page, y=Time)) +
  geom_boxplot() +
  labs(y='Czas (w sekundach)') +
  theme_bw()
graph

library(lmPerm)
summary(aovp(Time ~ Page, data=four_sessions))

### Statystyka F

summary(aov(Time ~ Page, data=four_sessions))

## Test chi-kwadrat
### Test chi-kwadrat: podejście randomizacyjne

clicks <- matrix(click_rate$Rate, nrow=3, ncol=2, byrow=TRUE)
dimnames(clicks) <- list(unique(click_rate$Headline), unique(click_rate$Click))

chisq.test(clicks, simulate.p.value=TRUE)

chisq.test(clicks, simulate.p.value=FALSE)

### Rozkład chi kwadrat (wykres)

x <- seq(1, 30, length=100)
chi <- data.frame(df = factor(rep(c(1, 2, 5, 20), rep(100, 4))),
                  x = rep(x, 4),
                  p = c(dchisq(x, 1), dchisq(x, 2), dchisq(x, 5), dchisq(x, 20)))

graph <- ggplot(chi, aes(x=x, y=p)) +
  geom_line(aes(linetype=df)) +
  geom_text(aes(x=25, y=0.06, label='st. swob.=20'), size=2.5) +
  geom_text(aes(x=8, y=0.1, label='st. swob.=5'), size=2.5) +
  geom_text(aes(x=4, y=0.2, label='st. swob.=2'), size=2.5) +
  geom_text(aes(x=1.5, y=0.05, label='st. swob.=1'), size=2.5) +
  theme_bw() + theme(legend.position = "none") +
  labs(x='Wartość statystyki chi-kwadrat', y='Prawdopodobieństwo')
graph

### Dokładny test Fishera

fisher.test(clicks)

#### Oszustwo naukowe

## Rysunek 3.8
imanishi$Digit <- factor(imanishi$Digit)

graph <- ggplot(imanishi, aes(x=Digit, y=Frequency)) +
  geom_bar(stat='identity') +
  theme_bw()
graph

## Moc i rozmiar próby
### Rozmiar próby

effect_size = ES.h(p1=0.0121, p2=0.011)
pwr.2p.test(h=effect_size, sig.level=0.05, power=0.8, alternative='greater')

effect_size = ES.h(p1=0.0165, p2=0.011)
pwr.2p.test(h=effect_size, sig.level=0.05, power=0.8, alternative='greater')
