```{r}
#| echo: false
#| message: false
#| warning: false

library(dplyr)
library(ggplot2)
library(broom)
```


# Inference for Difference of Two Means

## Hypothesis Testing:

**Example:**

Prices of diamonds are determined by what is known as the 4 Cs: cut, clarity, color, and carat weight. The prices of diamonds go up as the carat weight increases, but the increase is not smooth. For example, the difference between the size of a 0.99 carat diamond and a 1 carat diamond is undetectable to the naked human eye, but the price of a 1 carat diamond tends to be much higher than the price of a 0.99 diamond. In this question we use two random samples of diamonds, 0.99 carats and 1 carat, each sample of size 33, and compare the average prices of the diamonds. In order to be able to compare equivalent units, we first divide the price for each diamond by 100 times its weight in carats. That is, for a 0.99 carat diamond, we divide the price by 99. For a 1 carat diamond, we divide the price by 100. The distributions and some sample statistics are shown below.

|      | 0.99 carats | 1 carat |
|:-----|------------:|--------:|
| Mean |     \$44.51 | \$56.81 |
| SD   |     \$13.32 | \$16.13 |
| n    |          33 |      33 |

We need to perform a hypothesis test to evaluate if the average prices differ between the two groups.

### Hypotheses:

**Null Hypothesis (**$H_0$): There is no difference in the average standardized prices of 0.99 and 1 carat diamonds.

$H_0: \mu_.99 = \mu_1$

**Alternative Hypothesis (**$H_A$): There is a difference in the average standardized prices of 0.99 and 1 carat diamonds.

$H_A: \mu_.99 \neq \mu_1$

## Checking Assumptions:

Before we conduct the hypothesis test, we need to check some conditions for the t-test.

**Independence:** Both samples are random and represent less than 10% of their respective populations. Also, we have no reason to think that the 0.99 carats are not independent of the 1 carat diamonds since they are both sampled randomly.

**Normality**: n > 30 for both samples, hence we can assume that the distribution of the average differences will be nearly normal as well.


xbar0.99 - xbar1 ~ approx N(mu0.99 - mu1, ((sd0.99^2) / n0.99) + ((sd1^2) / n1))

## Calculations

### Sample statistics for 0.99 carat and 1 carat diamonds:

```{r}
mean_099 <- 44.51
sd_099 <- 13.32
n_099 <- 33

mean_1 <- 56.81
sd_1 <- 16.13
n_1 <- 33

```

### Difference in means:

```{r}
diff_means <- mean_099 - mean_1
diff_means
```

### Standard error:

```{r}
se <- sqrt((sd_099^2 / n_099) + (sd_1^2 / n_1))
se
```

### t-statistic (t-score):

t-score (test statistic) = ((xbar0.99 - xbar1) - (mu0.99 - mu1)) / se
= ((xbar0.99 - xbar1) - (0)) / se

```{r}
t_stat <- diff_means / se
t_stat
```

t-score ~ t distribution (df = min(n0.99, n1) - 1)
n1 = n2 = 33
df = 33 - 1 = 32

### Degrees of freedom:

```{r}
df <- 33 - 1
df
```

### Calculate p-value for two-tailed test:

```{r}
2*pt(t_stat, df)
p_value <- 2 * (1 - pt(abs(t_stat), df))
p_value
```

**Interpretation:** 

alpha = 0.05

Because our p-value = 0.0019, which is less than alpha = 0.05, we reject our null hypothesis and conclude that there is a difference in the average standardized prices of 0.99 and 1 carat diamonds. 

## Confidence Intervals:

### Difference in means:

```{r}
diff_means <- mean_099 - mean_1
diff_means
```

### Standard error:

```{r}
se <- sqrt((sd_099^2 / n_099) + (sd_1^2 / n_1))
se
```

### t critical value for 95% CI (degrees of freedom = 32):

```{r}
t_critical <- qt(0.975, df = 32)
qt(0.025, df = 32)
t_critical
```

### Confidence interval:

```{r}
ci_lower <- diff_means - t_critical * se
ci_upper <- diff_means + t_critical * se


ci_lower
ci_upper
```

mu0.99 - mu1

**Interpretation**: We are 95% confident that the true difference in average standardized prices of 0.99 and 1 carat diamonds falls within the interval (-19.72, -4.88). This interval is fully negative, which indicates that 1 carat diamonds are more expensive than 0.99 carat diamonds.

mu0.99 - mu1 > 0, then mu0.99 > mu1
mu0.99 - mu1 < 0, then mu0.99 < mu1

(-0.2, 0.5) <- includes 0, we can't say that there is a difference in average ____. 

(4.88, 19.72) <- indicates that 0.99 carat diamonds are more expensive 1 carate diamonds on average.


Unpaired vs paired:

unpaired: 2 independent samples (basically everything we went over above)
- difference between 2 means hypothesis test

paired: 1 person who does two different treatments/conditions. you want to compare the average difference in _____ between the two treatments/conditions (or before/after). 

```{r}
score_before <- c(80, 85, 92, 78)
score_after <- c(87, 85, 95, 80)
score_diff <- score_after - score_before
score_diff
```

- 1 sample mean hypothesis test



# Simple Linear Regression

## Dataset

The dataset **faithful** describes the waiting time between eruptions _(waiting)_ and the duration of the eruption _(eruptions)_ for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA. Both time are in minutes.

```{r}
glimpse(faithful)
```

## Linear Regression Setup

We would like to explore the linear relationship between the duration of the eruption (response variable) and the waiting time between eruptions (explanatory variable).

Does the waiting time between eruptions have a linear relationship with the duration of the eruption?

_Write out the linear equation between response variable and explanatory variable._

$$
eruptions = \beta_0 + \beta_1*waiting + \epsilon
$$

## Check the Conditions of Least Squares Regression


```{r}
lm_model <- lm(eruptions ~ waiting, data = faithful)

# The following lines stores residuals and predictions (y_hat) in the faithful data frame
faithful <- faithful |>
  modelr::add_residuals(lm_model) |> 
  modelr::add_predictions(lm_model)

glimpse(faithful)
```

### Linearity

```{r}
ggplot(faithful, aes(x = waiting, y = eruptions)) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE)
```

Looks generally linear (no curves). 

we can also see that it looks like there might be two groups of points.

### Constant Variance

```{r}
faithful |> 
  ggplot(aes(x = waiting, y = resid)) +
  geom_point() +
  geom_hline(yintercept = 0)
```

if there was constant variance, we would expect the residuals to not have any patterns in the above plot. here, we can see that there is a slight fanning pattern (like < symbol, going from small to wide in terms of how the points are distributed about the line in the middle), so there is not constant variance.

### Normality of Residuals

```{r}
faithful |> 
  ggplot(aes(x = resid)) +
  geom_density()
```

It is not quite normal. it looks like there are almost 2 peaks.

### Independence

_What about the independence of the observations?_

No, it is not fulfilled because all of the observations are from the same geyser. 

### Interpretation of Results

```{r}
lm_model_summary <- broom::tidy(lm_model)
lm_model_summary
```

$$
\hat{eruptions} = -1.87 + 0.0756*waiting
$$

_What are the point estimates for slope and intercept? How to interpret them?_

$b_0$ = -1.87

expected y value when our x value is 0

We expect the duration of the eruption to be -1.87 minutes when the waiting time is 0 minutes. 

$b_1$ = 0.0756

expected increase in our y variable for every 1 increase in our x variable.

We expect that for every increase of 1 minute in waiting time, the average duration of the eruption will increase by 0.0756 minutes. 

Note: if b_1 < 0: we expect that for every increase in 1 waiting time, the average duration of the eruption will decrease by ____ minutes. 

_What is the estimated eruption duration when the waiting time is 60 minutes? What is the estimated eruption duration when the waiting time is 120 minutes?_

estimating for 120 minutes is extrapolation because our waiting variable only ranges from 40-100. 

```{r}
-1.87 + 0.0756*60
```


### Inference

_Perform the hypothesis testing on whether the linear relationship between the duration of the eruption and waiting time between eruptions exists._

_Do we have significant evidence that the waiting time between eruptions affects the duration of the eruption?_

_Find the $95\%$ confidence interval for the slope. How to interpret this interval?_



```{r}
# point estimate for beta1
b1 <- as.numeric(lm_model_summary[2,2]) # This gets the second row and second column of the lm_model_summary table

# b0 <- as.numeric(lm_model_summary[1, 2]) # you can also use this if needed

critical_val <- qt(0.975, df = 271) 
se_b1 <- as.numeric(lm_model_summary[2,3])

me <- critical_val * se_b1
c(b1 - me, b1 + me)
```

```{r}
confint(lm_model)
confint(lm_model, level = 0.9)
```
