---
title: "Statistics"
output:
  ioslides_presentation:
    css: ../styles.css
    widescreen: yes
---

```{r knit-setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE, 
               message = FALSE, 
               warning = FALSE,
               fig.height = 4,
               fig.width = 7, 
               comment = "")
library(dplyr)
options(scipen=999)
library(readr)
library(ggplot2)
```

## Overview

We will cover how to use R to compute some of basic statistics and fit some basic statistical models. 

* Correlation
* T-test
* Linear Regression
* Logistic Regression

**DISCLAIMER**: We will focus on how to use R software to do these. We will be glossing over the statistical theory and "formulas" for these tests. Moreover, we do not claim the data we use for demonstration meet assumptions of the methods. 

There are plenty of resources online for learning more about these methods, as well as dedicated Biostatistics series (at different advancement levels) at the JHU School of Public Health.

# Correlation

## Correlation 

Function `cor()` computes correlation in R

```
cor(x, y = NULL, use = "everything",
    method = c("pearson", "kendall", "spearman"))
```

To use: 

- provide two numeric vectors (arguments `x`, `y`) to compute correlation between them, or
- provide matrix or data frame (argument `x`) that has at least 2 columns (must be numeric) to compute correlation between all pairs 

By default, Pearson correlation coefficient is computed. 

## Correlation {.small}

https://jhudatascience.org/intro_to_R_class/data/Charm_City_Circulator_Ridership.csv

```{r cor1, comment="", message = FALSE}
circ <- jhur::read_circulator()
head(circ)
```

## Correlation for two vectors

First, we compute correlation by providing two vectors.

Like other functions, if there are `NA`s, you get `NA` as the result.  But if you specify use only the complete observations, then it will give you correlation using the non-missing data. 

```{r}
x <- pull(circ, orangeAverage)
y <- pull(circ, purpleAverage)
```

```{r}
cor(x, y)
cor(x, y, use = "complete.obs")
```

## Correlation for two vectors with plot

Note that you can add the correlation value to a plot, via the `annotate()`. 

```{r, fig.width=3, fig.height=3}
cor_val <- cor(x, y, use = "complete.obs")
cor_val_label <- paste0("r = ", round(cor_val, 3))

circ %>%
  ggplot(aes(x = orangeAverage, y = purpleAverage)) +
  geom_point(size = 0.3) +
  annotate("text", x = 2000, y = 7500, label = cor_val_label, size = 5)
```

## Correlation for data frame columns

We can compute correlation for all pairs of columns of a data frame / matrix. We typically just say, *"compute correlation matrix"*.

Columns must be all numeric! 

```{r}
circ_subset_Average <- circ %>% select(ends_with("Average"))
dim(circ_subset_Average)

cor_mat <- cor(circ_subset_Average, use = "complete.obs")
cor_mat
```

## Correlation for data frame columns with plot

- Google, *"r correlation matrix plot"*

```{r, fig.width=4, fig.height=4}
library(corrplot)
corrplot(cor_mat, type = "upper", order = "hclust")
```

## Lab Part 1

[Website](http://jhudatascience.org/intro_to_R_class/index.html)

# T-test

## T-test

The commonly used are: 

- **one-sample t-test** -- used to test mean of a variable in one group 
- **two-sample t-test** -- used to test difference in means of a variable between two groups (if the "two groups" are data of the *same* individuals collected at 2 time points, we say it is two-sample paired t-test)

The `t.test()` function in R is one to address the above. 

```
t.test(x, y = NULL,
       alternative = c("two.sided", "less", "greater"),
       mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95, ...)
```

## Running one-sample t-test

It tests mean of a variable in one group. By default (i.e., without us explicitly specifying values of other arguments): 

- tests whether a mean of a variable is equal to 0 (`mu=0`) 
- uses "two sided" alternative (`alternative = "two.sided"`) 
- returns result assuming confidence level 0.95 (`conf.level = 0.95`)

```{r}
x <- pull(circ, orangeAverage)
t.test(x)
```

## Running two-sample t-test{.small}

It tests test difference in means of a variable between two groups. By default:

- tests whether difference in means of a variable is equal to 0 (`mu=0`) 
- uses "two sided" alternative (`alternative = "two.sided"`) 
- returns result assuming confidence level 0.95 (`conf.level = 0.95`)
- assumes data are not paired (`paired = FALSE`)
- assumes true variance in the two groups is not equal (`var.equal = FALSE`)

```{r}
x <- pull(circ, orangeAverage)
y <- pull(circ, purpleAverage)
t.test(x, y)
```

## T-test: retrieving information from the result 

Object returned from `t.test()` function is a named list. We can use it to access test result elements. The easiest way to do this is to use base R (`$` notation).

```{r}
result <- t.test(x, y)
is.list(result)
names(result)
result$statistic
result$p.value
```

## T-test: retrieving information from the result with `broom` package 

The `broom` package has a `tidy()` function that can organize results into a data frame so that they are easily manipulated (or nicely printed)

```{r broom, comment=""}
library(broom)

result <- t.test(x, y)
result_tidy <- tidy(result)
result_tidy
```

## Some other statistical tests

- `wilcox.test()` -- Wilcoxon signed rank test, Wilcoxon rank sum test
- `shapiro.test()` -- Shapiro test
- `ks.test()` -- Kolmogorov-Smirnov test 
- `var.test()`-- Fisher’s F-Test
- `chisq.test()` -- Chi-squared test

## Lab Part 2

[Website](http://jhudatascience.org/intro_to_R_class/index.html)

# Regression

## Linear regression

Linear regression is a method to model the relationship between a response and one or more explanatory variables. 

We provide a little notation here so some of the commands are easier to put in the proper context.

$$
y_i = \alpha + \beta x_{i} + \varepsilon_i 
$$
where:

* $y_i$ is the outcome for person i
* $\alpha$ is the intercept
* $\beta$ is the slope 
* $x_i$ is the predictor for person i
* $\varepsilon_i$ is the residual variation for person i

## Linear regression

```{r,echo=FALSE}
coefs=lm(data=iris,Petal.Width~Petal.Length) %>%  coef

line_df=data.frame(x=c(0,7),y=c(coefs[1],(coefs[2]*7+coefs[1])))
riserun_df=data.frame(x=c(2.1,2.9,2.9),y=c((coefs[1]+coefs[2]*2.1),(coefs[1]+coefs[2]*2.1),(coefs[1]+coefs[2]*2.9)))
residual_df=data.frame(x=c(3,3),y=c(1.1,(coefs[1]+coefs[2]*3)))
labels_df=data.frame(x=c(0.5,3.2,2.7),y=c(coefs[1],0.65,1),labels=c("alpha","beta","epsilon[i]"))   
 
print(ggplot(data=iris,aes(x=Petal.Length,y=Petal.Width))+geom_point()+geom_line(data=line_df,aes(x=x,y=y))+
    geom_point(data=line_df[1,],aes(x=x,y=y,color="red",size=1.5))+
    geom_line(data=riserun_df,aes(x=x,y=y,color="orange",size=1.4))+
    geom_line(data=residual_df,aes(x=x,y=y,color="blue",size=1.4))+
    geom_label(data=labels_df,aes(x=x,y=y,label=labels,size=2),parse=TRUE)+
    theme_bw()+ylab(expression(y))+xlab(expression(x))+
    theme(legend.position="none",text=element_text(size=20)))

```

## Linear regression

Linear regression is a method to model the relationship between a response and one or more explanatory variables. 

We provide a little notation here so some of the commands are easier to put in the proper context.

$$
y_i = \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \varepsilon_i 
$$
where:

* $y_i$ is the outcome for person i
* $\alpha$ is the intercept
* $\beta_1$, $\beta_2$, $\beta_2$ are the slopes for variables $x_{i1}$, $x_{i2}$, $x_{i3}$ 
* $x_{i1}$, $x_{i2}$, $x_{i3}$ are the predictors for person i
* $\varepsilon_i$ is the residual variation for person i

## Linear regression fit in R 

To fit linear models in R, we use function `lm()`. 

```
lm(formula, data, subset, weights, na.action,
   method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
   singular.ok = TRUE, contrasts = NULL, offset, ...)
```

We typically provide two arguments: 

- `formula` -- model formula written using names of columns in our data 
- `data` -- our data frame 

## Linear regression fit in R: model formula 

Model formula 
$$
y_i = \alpha + \beta x_{i} + \varepsilon_i 
$$
translates to `y ~ x` in R formula for this example. 

In practice, `y` and `x` are replaced with the **names of columns from our data set**. 

- For example, if we want to fit a regression model where outcome is `income` and predictor is `years_of_education`, our formula would be: 

  `income ~ years_of_education`

## Linear regression fit in R: model formula 

Model formula 
$$
y_i = \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \varepsilon_i 
$$
translates to `y ~ x1 + x2 + x3` in R formula for this example. 

In practice, `y` and `x1`, `x2`, `x3` are replaced with the **names of columns from our data set**. 

- For example, if we want to fit a regression model where outcome is `income` and predictors are `years_of_education`, `age`, `location` then our formula would be: 

  `income ~ years_of_education + age + location`

## Linear regression 

We will use `kaggleCarAuction.csv` dataset from one of the Kaggle competitions. 

https://jhudatascience.org/intro_to_R_class/data/kaggleCarAuction.csv

```{r}
cars <- jhur::read_kaggle()
head(cars)
```

## Linear regression: model fitting 

We fit linear regression model with  vehicles odometer (distance traveled by a vehicle; `VehOdo`) as an outcome and vehicle (`VehicleAge`) age as a predictor. 

```{r}
fit <- lm(VehOdo ~ VehicleAge, data = cars)
print(fit)
```

## Linear regression: model summary

The `summary()` command returns a list that shows us some more detail

```{r}
sfit <- summary(fit)
print(sfit)
```

## Linear regression: retrieving information with `broom` package 

Use `tidy` to create a tibble with the coefficient estimates.

```{r}
# tidy() is a function from the broom package
tidy(sfit)
```

## Linear regression: model summary

Model summary is a named list and we can access its specific elements. Again, we should use base R (`$` notation).

```{r}
names(sfit)
sfit$r.squared
```

## Linear regression: multiple predictors {.smaller}

Let's try adding another explanatory variable to our model, Warranty price (`WarrantyCost`) 

```{r}
fit_2 <- lm(VehOdo ~ VehicleAge + WarrantyCost, data = cars)
summary(fit_2)
```

## Linear regression: factors 

Factors get special treatment in regression models - lowest level of the factor is the comparison group, and all other factors are relative to its values.

`TopThreeAmericanName` states if the manufacturer is one of the top three American manufacturers.

```{r}
top_3 <- pull(cars, TopThreeAmericanName)
table(top_3)
```

## Linear regression: factors {.smaller}

```{r regress6, comment="", fig.height=4,fig.width=8}
fit_3 <- lm(VehOdo ~ factor(TopThreeAmericanName), data = cars)
summary(fit_3)
```

## Linear regression: retrieving information with `broom` package 

```{r, comment="", fig.height=4,fig.width=8}
tidy(fit_3)
```

## Generalized Linear Models (GLMs)

Generalized Linear Models (GLMs) allow for fitting regressions for  non-continuous/normal outcomes. Examples include: logistic regression, Poisson regression. 

We fit GLM with a `glm()` function that has a very similar syntax to the `lm()` function. 

The primary difference is in `glm()`, we additionally specify the `family` argument -- a description of the error distribution and link function to be used in the model. These include: 

- `binomial(link = "logit")`
- `poisson(link = "log")`, and other.

See `?family` documentation for details of family functions. 

## Logistic regression {.smaller}

`IsBadBuy` is a 0/1-valued variable stating "if the kicked vehicle was an avoidable purchase". 

```{r regress7, comment="", fig.height=4,fig.width=8}
glm_fit <- glm(IsBadBuy ~ VehOdo + VehicleAge, data = cars, family = binomial())
summary(glm_fit)  
```

## Final note 

Some final notes: 

- Researcher's responsibility to **understand the statistical method**  they use -- underlying assumptions, correct interpretation of method results 

- Researcher's responsibility to **understand the R software**  they use -- meaning of function's arguments and meaning of  function's output elements

## Lab Part 3

[Website](http://jhudatascience.org/intro_to_R_class/index.html)
