# Comparison of Samples
## Preamble

```{r}
# Preamble
## Install Pacman
load.pac <- function() {
  
  if(require("pacman")){
    library(pacman)
  }else{
    install.packages("pacman")
    library(pacman)
  }
  
  pacman::p_load(xts, sp, gstat, ggplot2, rmarkdown, reshape2, ggmap,
                 parallel, dplyr, plotly, tidyverse, reticulate, UsingR, Rmpfr,
                 swirl, corrplot, gridExtra, mise, latex2exp, tree, rpart, lattice)
  
}

load.pac()
mise()
load(file = "./TAD.rdata")
```

Load the Data

```{r}
(birthwt <- as_tibble(read.csv(file = "./0datasets/birthwt.csv", header = TRUE, sep = ",")))

 # Always use doubles / numeric unless the data is legitamately
   # restricted to the integers, it will break functions ( and so it should)
birthwt$bwt <- as.numeric(birthwt$bwt)

birthwt$smoke <- c(FALSE, TRUE)[birthwt$smoke]

summary(birthwt)
str(birthwt)
dim(birthwt)
nrow(birthwt)
ncol(birthwt)
```


## Smoking and Birth Weight
### Summarise Data
#### Table

```{r}
table(birthwt$smoke)
table(birthwt$smoke) %>% barplot(main = "Birth Weight and Smoking")
```

#### Barplot

```{r}

desc_stats <- function(x) {
  mean(x)
  median(x)
  var(x)
  sd(x)
}

(desc_stats <- data.frame(
mean = apply(birthwt, 2, mean),
median = apply(birthwt, 2, median),
var = apply(birthwt, 2, var),
sd = apply(birthwt, 2, sd)
))

```

##### Range

```{r}
range(birthwt$bwt)
range(birthwt$bwt) %>% diff()
max(birthwt$bwt) - min(birthwt$bwt)
```

###### Quantile

The quantile function returns $x$ -axis values corresponding to a what proportion of the data is specified, so for example, for a standard normal distribution $\mathcal{N}~\left( 0, 1 \right)$, 2.5% of the observations lie below 2 and another 2.5% lie above 2.

```{r}
quantile(rnorm(1000), 0.025)

quantile(birthwt$bwt, 0.25)
quantile(birthwt$bwt, 0.75)
```

###### Inter-Quartile Range

This can be calculated thusly:

```{r}
IQR(birthwt$bwt)
```


For normally distributed data we would expect:

$$
\textsf{IQR} = 1.349 \times \sigma
$$

Remember tha the normal distribution is modelled using calculus:

$$
\begin{aligned}
f(x) &= - \sqrt{\frac{k}{2\pi}}  \cdot   e^{k\cdot  \frac{\left( x- \mu \right)^2}{2}} \\
    f\left( x \right)&= \sqrt{\frac{1}{2\pi}}\cdot  \sum^{\infty}_{n= 0}   \left[\frac{ \left( - \frac{1}{2}z^2 \right)^n}{n!} \right]  \\
\int f\left( x \right) \mathrm{d}x &= \frac{1}{\sqrt{2\pi} }\int \sum^{\infty}_{n= 0}   \left[ \frac{\left( - \frac{1}{2}z^2 \right)^n}{n!} \right]  \mathrm{d}z \\
&= \frac{1}{\sqrt{2\pi} }\cdot  \sum^{\infty}_{n= 0}   \left[ \int \frac{\left( - 1 \right)^{- 1}z^{2n}}{2^n\cdot  n!} \mathrm{d}z  \right] \\
&= \frac{1}{\sqrt{2\pi} }\cdot  \sum^{\infty}_{n= 0}   \left[ \frac{\left( - 1 \right)^n \cdot  z^{2n+ 1}}{2^n\left( 2n+  1 \right)n!} \right] 
\end{aligned}
$$


#### Histograms

A histogram would offer a better understanding of the data:

```{r}
x <- birthwt$bwt
hist(birthwt$bwt, main = "Histogram of Birth Weight Given Smoking", col = "lightblue", border = "indianred", freq = FALSE)
# curve(dnorm(x, mean(x), sd(x)), add = TRUE)
```

Adding a Density urve is extremely difficult in base plot, it's so much easier to use ggplot2:

```{r}
birthwt_pretty <- birthwt
birthwt_pretty$smoke <- ifelse(birthwt$smoke, "Smoking", "non\nSmoking")

hist <- ggplot(birthwt_pretty, aes(x = bwt, fill = smoke, col = "black", y = ..density..)) +
         theme_classic() + 
         labs(x = "Birth Weight", y = "Density")

plots <- list()

# Dodge
plots[[1]] <- hist + geom_histogram(position = "dodge2", col = "blue", binwidth = 300)

# Overlay
plots[[2]] <-  hist + geom_histogram(binwidth = 300, col = "black")

# Single Histogram
plots[[3]] <- hist + geom_histogram(binwidth = 300, col = "black", aes(group = 1), fill = "lightblue") 

# Facet Grid (i.e. Split Chart)
 plots[[4]] <- hist + geom_histogram(binwidth = 300, col = "black") +
         facet_grid(. ~ smoke) +
         guides(fill = FALSE)

  
  # Colour it
  # Make a Facet Grid
  # Add a Density Curve

```

```{r}
layout <- matrix(c(1, 1, 2, 3, 4, 4), byrow = TRUE, nrow = 3)
# arrangeGrob(grobs = plots, layou_tmatrix = layout)
grid.arrange(grobs = plots, layout_matrix = layout)


```


#### Splitting the Data Up

The `aggregate()` function can split up data by a specified variable given a `formula`:

```{r}
bwt_form <- bwt~smoke

# Formula as a variable
aggregate(bwt_form, birthwt, mean)

# Otherwise
aggregate(bwt~smoke, birthwt, sd)
aggregate(bwt~smoke, birthwt, median)
```
It would also be possible to do this manually but it's actually quite tricky:
> Tibble's will always remain as tibble's when a column is extracted and will not yield to `as.vector()`, it is necessary to instead use the `dplyr::pull()` function or first conver the tibble to a matrix.


```{r}

matrix(
  c(
"Mean_Smokers"  = as.matrix(birthwt)[birthwt$smoke==1,"bwt"] %>% mean(),
"Mean_NonSmoke" = as.matrix(birthwt)[birthwt$smoke==0,"bwt"] %>% mean()
))

```


#### Split Charts

##### Base

Split Charts are essentially a `facet grid` of two historgrams corresponding to a variable, the can be generated using the `lattice` library which is usually bundled with base **_R_**: 

```{r}
library(lattice)
histogram(~bwt|smoke, data = birthwt)
```



##### GGPlot2 

This was shown previously:

```{r}
# Facet Grid

ggplot(birthwt_pretty, aes(x = bwt, fill = smoke, col = "black", y = ..density..)) +
         geom_histogram(binwidth = 300, col = "black") +
         facet_grid(. ~ smoke) +
         guides(fill = FALSE) +
         theme_classic() + 
         labs(x = "Birth Weight", y = "Density")
```

#### Box Plots

##### Base 
This is varily easy if the data is such that one variable is a factor/logical, the syntax uses a `formula`, so in this case `Y ~ X` would by *birthweights depends on smoking* [^1]:

[^1]: As opposed to smoking depending on birthweight.
###### Vertical
```{r}
birthwt_pretty <- birthwt
birthwt_pretty$smoke <- ifelse(birthwt$smoke, "Smoking", "Non-Smoking")
boxplot(bwt ~ smoke, birthwt_pretty,
        col = "lightblue", main = "Comparison of Birthweight Given Smoking",
        xlab = "", ylab = "Birth Weight", cex = 2, pch = 3)
```

###### Horizontal

```{r}
birthwt_pretty <- birthwt
birthwt_pretty$smoke <- ifelse(birthwt$smoke, "Smoking", "Non-Smoking")
boxplot(bwt ~ smoke, birthwt_pretty,
        col = "lightblue", main = "Comparison of Birthweight Given Smoking",
        xlab = "", ylab = "Birth Weight", cex = 2, pch = 3, horizontal = TRUE)
```

##### GGPlot2

```{r}
ggplot(birthwt, aes(x = smoke, y = bwt, fill = smoke)) +
  geom_boxplot(show.legend = FALSE) +
  scale_x_discrete(labels = c("Non\nSmoking", "Smoking")) + 
  theme_classic() + 
  labs(x = "", y = "Birth Weight", title = "Comparison of Birthweight Given Smoking Habits")
  
  
  
```

### Difference in Means

```{r}
(birthwt_agg <- aggregate(bwt ~ smoke, data = birthwt, FUN = mean))

# Diff is the difference like the delta or the derivative in calculus or real anlysis so X2-X1
 # this represents the change from those that smoke to those that dont

# yes - no
diff(birthwt_agg$bwt)
```

In order to simulate a test statistic consider the difference in the mean values, under the null hypothesis this would be zero.

#### Simulation

Under the Null Hypothesis the Birth Weight is identical regardless of smoking
status, so the simulation would simply consider random permutations of smoking
across the birth weights.

The number of times the birthwt is lower in the observation that that expected
by chance is the rate of false positives, the expected value of which is the
p-value:

> recall that the false positive rate is the rate at which false positives occur
> from all negative responses, i.e. the rate of false positives assuming that
> the null hypothesis is true.

```{r}
smoke.sim <- birthwt
smoke.sim$smoke <- sample(birthwt$smoke, replace = FALSE)

(smoke.sim_agg <- aggregate(bwt ~ smoke, smoke.sim, mean))


# yes - no
diff(smoke.sim_agg$bwt)
```

A False Positive occurs when the random value is even more extreme than what we observed, this measures the frequency of false positives for this observation, i.e. the p-value corresponding to this sample.

#### Replicate the Simulation

Replicate uses a `apply` and automatically creates a vector (statically not
dynamicall), so it's usually better to use than a for loop.

```{r}
## Create the Simulation

sim_diff <- birthwt
# smoke.sim <- birthwt
simulation <- function() {
  replicate(n = 10^3, expr = {
    comparison()
  }
 )
}

comparison <- function(){
smoke.sim$smoke <- sample(birthwt$smoke, replace = FALSE)
smoke.sim_agg   <- aggregate(bwt ~ smoke, smoke.sim, mean)

diff(smoke.sim_agg$bwt) 
}

# using Abs is a good healthy sanity check
obs     <- aggregate(bwt ~ smoke, birthwt, FUN = mean) %>%
  pull(bwt) %>%
  diff() %>%
  abs()*-1

simdata <- simulation()
simdata.tb <- tibble::enframe(simdata, value = "Weight")
```


##### Create a Histogram
```{r}
# Plot a Histogram
  ## Base
   ### Be mindful to set xlim
simdata %>% hist(main = "Expected Difference in Birthweights\n Assuming No Population Difference",
                 xlim = c(-abs(obs), abs(obs)))
abline(v = obs)

## ggplot
ggplot(data = simdata.tb, aes(x = Weight, y = ..density..)) +
    geom_histogram(fill = "#Eea2ad", col = "#1e90ff" ) +
    theme_bw() +
    labs(y = "Frequency", title = "Simulated Samples Under null hypothesis") +
    # Use aes() to make legend
    geom_vline(lwd = 0.8, lty = 2, aes(xintercept = obs, col = "Observation")) +
    scale_color_manual(name = "Observation", values = c(Observation= "#D02090", mean = "red"), label = "Observed Value")
```


##### P-value

The rate at which our observation exceeds the simulation (which is performed under the assumption that there is no difference) is the false positive rate, this is the expected value of the p-value.

The p-value is the probability of a type 1 error, the probability of rejecting the null hypothesis, under the assumption that it is true.

```{r}
(pval <- mean(sim_diff < obs))
```

Hence the probability of asserting that:

  * smoking decreases birthweight

Assuming that smoking caused no difference, is too small to detect and hence we conlude that smoking decreases birthweight.


### Simulating Data (All Together)

```{r}
## Create the Simulation
smoke.sim <- birthwt
sim_diff <- replicate(n = 10 ^ 3, expr = {
  # Shuffle the Smoking Status
  smoke.sim$smoke <- sample(birthwt$smoke, replace = FALSE)
  # Aggregate the average Value
  smoke.sim_agg   <- aggregate(bwt ~ smoke, smoke.sim, mean)
  
  # Determine the Difference
  diff(smoke.sim_agg$bwt) 
})

# Select for a data frame; Pull out a vector
obs <- aggregate(bwt ~ smoke, birthwt, mean) %>%
  pull(bwt) %>%
  diff()

# How often is a value equally extreme observed?
# Question on Two Sided vs One Sided

pval <- mean(sim_diff < -obs)

```

## A note on Hypotheses Formation

The null hypothesis $H_0$ is always that the two populations, as indicated by the mean value of the sample statistic, is zero, this is independent of the alternative hypothesis.

The alternative hypothesis depends entirely on the direction of the teststatistic, if the false positives are tested for above below or both above/below determines the whether or not the alternative hypothesis is upper/lower tail or two-sided.


## Challenge Data


## East and West

### What is the average difference?

```{r, warning=FALSE}
load()
# attach(salesEW) # (would warn me (but don't use the variable assignments it makes,
                   # that's asking for trouble, but do use the autocomplete)
names(salesEW)
xbar_obs <- aggregate(sales ~ office, data = salesEW, FUN = mean)[,2] %>% diff() %>% abs() 
```

### What is the simulated average difference

```{r}
  sales_shuffle <- salesEW
xbar_sim <- replicate(10^3, {
  sales_shuffle$office <- sample(salesEW$office)
  aggregate(sales ~ office, sales_shuffle, mean)[,2] %>% diff() # NEVER Abs Here it doesn't make sense!%>% abs
})

```

### What is the P-Value

```{r}
#How often is the simulated difference really negative?
p_12 <- mean(xbar_sim < -xbar_obs) + mean(xbar_sim > xbar_obs)
#How often is the simulated difference really Positive?
p_12 <-mean(xbar_sim > xbar_obs)

# Putting them together what's the total FPR for any difference?

(pval <- mean(xbar_sim < -xbar_obs) + mean(xbar_sim > xbar_obs))
```

Hence the probability of incorrectly asserting there is a difference in sales between the two locations (assumint there isn't one) is about 5% and so the the null hypothesis is rejected and it is concluded that there is indeed a difference in sales between the two offices.

### Plot the Data

```{r}
ggplot(salesEW, aes(x = office, y = sales)) +
  geom_boxplot()
```



## Spiders

```{r}
load(file = "~/Notes/DataSci/ThinkingAboutData/TAD.rdata")
spider <- read.csv(file = "./0datasets/Spider.csv")
str(spider)
attach(spider)
```

```{r}
g <- diff(aggregate(Anxiety ~ Group, spider, mean)$Anxiety)
g0 <- replicate(10^3, {
  groupShuf <- sample(spider$Group)
  diff(aggregate(Anxiety~groupShuf, spider, mean)$Anxiety)
})

(pVal <- mean(g0 > g))

detach(spider)
```



