---
title: "BIOS:4120 - Lab 11: Fisher Test, CI for OR, and Quiz 3 Review"
output:
  html_document:
    toc: TRUE
    toc_float: TRUE
date: "April 2-3, 2024"
editor_options: 
  markdown: 
    wrap: 72
---

## Objectives

1.  Revisit Fisher's Exact Test
2.  Discuss the confidence intervals for the odds ratio
3.  Quiz 3 Review

#### Note: The material in the first part of this lab will not be on quiz 3.

## Fisher Test and Confidence Interval

### Contingency tables with Lister data

Let's start by reviewing the lister dataset.

```{r}
lister <- read.delim('https://github.com/IowaBiostat/data-sets/raw/main/lister/lister.txt')
table(lister) # make table of the data
```

Do not fret that this table seems like it's in reverse. It turns out
that the table being in this "reverse" order doesn't actually matter for
when we do our analysis, but to look at the data as it was displayed
during lecture, we can switch the row and column order with a bit of
manipulation to the dataset.

```{r}
lister$Group <- factor(lister$Group, levels = c("Sterile", "Control")) # specify order you want in the levels
lister$Outcome <- factor(lister$Outcome, levels = c("Survived", "Died"))
(lister_tab <- table(lister)) # we will use this table in later analysis
```

Alternatively, if we knew the counts of the data themselves, we could
make a table in the following way:

```{r}
lister_manual <- rbind(c(34, 6),
                       c(19, 16))
# Add row and column names
rownames(lister_manual) <- c("Sterile", "Control")
colnames(lister_manual) <- c("Survived", "Died")
lister_manual

# Alternative way to do this:
lister_manual2 <- data.frame("Group" = c(rep("Control", 35),
                                        rep("Sterile", 40)),
                            "Outcome" = c(rep("Survived", 19),
                                          rep("Died", 16),
                                          rep("Survived", 34),
                                          rep("Died", 6)))

# specify order of rows and columns in the levels:
lister_manual2$Group <- factor(lister_manual2$Group, levels = c("Sterile", "Control"))
lister_manual2$Outcome <- factor(lister_manual2$Outcome, levels = c("Survived", "Died"))

table(lister_manual2)

```

**Note: You may need to do something like this for your next homework
assignment**

<br>

### Fisher's exact test and odds ratio

Remember that Fisher's Exact Test can be performed using the
fisher.test() function. In lab last week using the lister dataset, we
saw that this function not only gives you the p-value but also the odds
ratio (OR) and a 95% CI for the odds ratio.

```{r}
fisher.test(lister_tab)
```

Note that this Fisher's Exact Test odds ratio is calculated using the
formula: \newline

\Large

$\frac{a/b}{c/d}=\frac{ad}{bc}$ \normalsize, where a, b, c, and d are
defined by their location according to the following table:

```{r,echo=FALSE}
Table<-matrix(data=c("a","b","c","d"),nrow=2,byrow=TRUE)
rownames(Table)<-c("Option 1","Option 2")
colnames(Table)<-c("Success","Failure")
print(Table, quote = FALSE)
```

The easiest way to think of this odds ratio is to say:

If OR \> 1: The odds of success are (OR-1)\*100 percent higher if we do
Option 1 than if we do Option 2.

If OR \< 1: The odds of success are (1-OR)\*100 percent lower if we do
Option 1 than if we do Option 2.

If OR = 1: The odds of success are equal between Option 1 and Option 2.

**Be careful when calculating the OR that you interpret what you
calculated. For instance, if I were to switch my Success and Failure
column names, I'd be calculating the OR for Failure given Option 1
instead of Success. If you switch both the rows and column
simultaneously (like we did with the very first table), the odds ratio
will be unchanged.**

### Confidence interval for odds ratio

Calculating a confidence interval for an odds ratio is slightly more
complicated than what we've done so far because it's on a different
scale. Here are the steps:

1.  Find the odds ratio: $\frac{ad}{bc}$

2.  Find the log odds ratio (**natural log**): $\log{(OR)}$

3.  Find the error term:
    $SE_{logOR}=\sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}}$

4.  Calculate the CI on this scale:
    $CI_{logOR}=\log{(OR)} \pm z_{\alpha/2} * SE_{logOR}$

5.  Calculate the CI on the original scale: $CI = \exp{(CI_{logOR})}$

Let's find a CI for the OR for survival given sterile procedure by hand.

```{r,echo=FALSE}
print(lister_tab)
table2<-lister_tab
table2[1:2,1:2]<-c("a","c","b","d")
print(table2)
```

1.  Find the OR.

```{r}
OR <- (34*16) / (19*6)

OR
```

2.  Find the log(OR)

```{r}
logOR <- log(OR) # the log function in R is the natural log

logOR
```

3.  Find the error term.

```{r}
SElogOR <- sqrt((1 / 19) + (1 / 6) + (1 / 16) + (1 / 34))

SElogOR
```

4.  Calculate the CI on this scale.

```{r}
logCI <- logOR + qnorm(c(0.025, 0.975)) * SElogOR

logCI
```

5.  Calculate the CI on the original scale.

```{r}
CI <- exp(logCI)
CI
```

**Interpretation:** We can say with 95% confidence that the true odds
ratio for survival with sterile procedure relative to the control
procedure is (`r round(CI,3)`), indicating significantly increased odds
of survival on the sterile procedure compared to the non-sterile
procedure. (Significant because the CI does not include the value 1.)

Note that this is the same as the relative odds of dying with the
control surgery. This happens because the OR is symmetric.

## Quiz 3 Review

A flowchart for when to use each distribution

![](flowchart.png)

### General approach to stats problems

#### **Step 1: Figure out and list the given information.**

Write out the values of $n, \bar{x}, s$, etc. for everything given in
the prompt.

<br>

#### **Step 2: Figure out what the problem is asking. Write the null hypothesis (if applicable).**

What is the population we're looking at? What kind of distributions can
be used with the information given? Which is most appropriate in this
situation? (see chart) Do we want a test or an interval? Most
importantly, what is the question being asked?

**If Defining Hypotheses:** The null case is the one based on no change,
difference, or improvement. The alternative is what we want to show and 
that there is a change or difference. Remember that the hypotheses describe and
apply to a population parameter.

<br>

#### **Step 3: Write out the equation for the test statistic or interval.**

In this class, most test statistics will be similar to the form:
$z=\frac{(\hat{X}-X_0)}{SE_{\hat{X}}}$, and most confidence intervals
will look similar to this: $\hat{X} \pm z^**SE_{\hat{X}}$,

where:

-   $X$ is the parameter of interest\
-   $\hat{X}$ is the estimate of the parameter\
-   $X_0$ is the value of the parameter under the null hypothesis\
-   $SE_{\hat{X}}$ is the standard error of the estimate\
-   $z^*$ is the critical value for the interval (either from the z or t
    distribution)

<br>

#### **Step 4: Plug in the values.**

Calculate the statistic or interval using the given values. Find the
p-value (if applicable). This means that if we're conducting a test,
compare the test statistic to the distribution we picked in step 2, and
find the probability of being as or more extreme than the value we
observed.

<br>

#### **Step 5: Interpret in the context of the problem.**

This might very well be the most important step. Our interpretation is
dependent on our study, so it varies. Below are some examples.

*Situation: Paired study looking at the effect of oat bran consumption
compared to cornflakes levels where the calculated p-value was 0.005*

-   "There is strong evidence that eating oat bran as opposed to corn
    flakes leads to lower cholesterol levels."

*Situation: Weights of U.S. Males (mean believed to be 172.2) where
observed mean was 180, and the calculated p-value was 0.09.*

-   "This study provides borderline evidence that the true
    mean weight of males in the United States is greater than 172.2
    pounds."

**Table 1. Modifiers based on p-value:**

| p     | Evidence against null |
|-------|-----------------------|
| \> .1 | insufficient          |
| .1    | borderline            |
| .05   | moderate              |
| .01   | strong                |
| .001  | overwhelming          |

<br>

### Concepts to know

#### Effect size

Recall the Nexium/Prilosec example. If we have a large enough n, we can
find statistical significance in the smallest true difference in means.
Here, the difference, while real, was only 3%. That's effectively REALLY
close to being the same. For practical purposes, this difference doesn't
matter, but the p-value doesn't tell us that the actual difference is
little, just that it is significant.

#### Paired data

Paired data can be analyzed using the binomial distribution (proportion
of patients that saw any improvement) or using continuous methods (using
1-sample methods on a set of the differences between the two).

#### Power

Given that the null hypothesis is false, the power of a test is the
probability that we will correctly reject the null.

Properties of power:

-   Power = (1 - $\beta$), where $\beta$ is the type II error rate
-   Power increases as sample size increases
-   Power increases as type I error ($\alpha$) increases
-   Power decreases as standard deviation increases
-   Power increases as effect size increases

## Practice Problems

#### 1. We are interested in testing whether a certain at-risk population for diabetes has a daily sugar intake that is equal to the general population, which is equal to 77 grams/day. A sample of size 37 was taken from this at-risk population, and we obtained a sample mean of 80 and sample standard deviation of 11 grams. Perform a hypothesis test to test whether this population has a significantly different mean sugar intake from 77 grams.

<details>

<summary>Answer</summary>

We have continuous data and we have the sample standard deviation(s),
so we should perform a t-test. The null hypothesis is that the true mean
daily sugar intake is 77 grams/day.

$$
\bar{x} = 80\\
\mu_0 = 77\\
n = 37\\
s = 11\\
df = 37-1=36
$$

$$
SE = \frac{s}{\sqrt{n}} = \frac{11}{\sqrt{37}} = 1.808\\
t = \frac{\bar{x}}{SE} = \frac{80-77}{1.808} = 1.659
$$

The p-value is the probability of observing a value as extreme or more
extreme than the one we observed under the null hypothesis. To compute
it, find 2\*Pr(T \> 1.659). Using the t table, look at the row
corresponding to our df. Since 36 is not on the table, you can just
round to the nearest df you see on the table (or always round down if
you like).

In that case, let's look at df = 35, our test statistic t = 1.66 is
between 1.47 and 1.69 which corresponds to the alpha values 0.15 and 0.1
so we can report the p-value in the following way:

0.1 \< p-value \< 0.15

**Recall that the t table Patrick provides gives the area OUTSIDE the
t-value so the p-value directly corresponds to the alpha values at the
top of the table**

We conclude that we do not have significant evidence to conclude that
the at-risk mean sugar intake is different from the general population
mean sugar intake of 77 grams/day (0.1 \< p \< 0.15).

```{r}
mu <- 77
xbar <- 80
s <- 11
n <- 37

# Use a t-test for the sample

t <- (xbar-mu) / (s / sqrt(n))

2*pt(abs(t), n-1, lower.tail=FALSE)
```

</details>

<br>

#### 2. The distribution of LDL cholesterol levels in a certain population is approximately normal with mean 90 mg/dl and standard deviation 8 mg/dl.

<!-- -->

##### a) What is the probability an individual will have a LDL cholesterol level above 95 mg/dl? 

<details>

<summary>Answer</summary>

The data are approximately normal and we are looking at INDIVIDUALS.

$$
x = 95 \\
\mu = 90 \\
\sigma = 8 \\
P(Z > \frac{x - \mu}{\sigma}) \\
= P(Z > \frac{95-90}{8}) = P(Z > 0.625) \\
= 1 - P(Z < 0.625) = P(Z <-0.625) \\
= 1 - 0.736 = 0.264
$$

```{r}
mu <- 90
sigma <- 8

(p <- pnorm(95, mu, sigma, lower.tail=FALSE))
```

</details>

<br>

##### b) Suppose we have a sample of 10 people from this population. What is the probability of exactly 3 of them being above 95 mg/dl? 

<details>

<summary>Answer</summary>

We now can use the binomial distribution to find this probability where
$\pi$ is the probability of an INDIVIDUAL being above 95 mg/dl. We found
exactly this probability in part a.

$$
X \sim Binom(n = 10, \pi = 0.264) \\
p(x) = \frac{n!}{x!(n-x)!}\pi^x(1-\pi)^{n-x} \\
p(3) = \frac{10!}{3!(7)!}0.264^3(1-0.264)^{7} \\
p(3) = 0.258
$$

```{r}
dbinom(3, 10, p)
```

</details>

<br>

##### c) Take the sample of size 10, as in part b. What is the probability that the sample mean will be above 95 mg/dl? 

<details>

<summary>Answer</summary>

The data are still approximately normal and we are interested in the
SAMPLE MEAN being above 95 mg/dl. This means we need to use the Z
formula but use SE instead of the population sd ($\sigma$)

$$
x = 95 \\
\mu = 90 \\
\sigma = 8 \\
n = 10 \\
SE = \frac{\sigma}{\sqrt{n}} = \frac{8}{\sqrt{10}} = 2.529822 \\
P(Z > \frac{x - \mu}{SE}) \\
= P(Z > \frac{95-90}{2.529822}) = P(Z > 1.976) \\
= 1 - P(Z < 1.976) = P(Z < -1.976) \\
= 1 - 0.976 = 0.024
$$

```{r}
z <- (95-90)/(sigma/sqrt(10))

(p <- pnorm(z, lower.tail=FALSE))
```

</details>

<br>

##### d) Suppose we take 5 samples of size 10 from the population. What is the probability that at least one of the sample means will be greater than 95 mg/dl?

<details>

<summary>Answer</summary>

Similarly to part b, we can use the Binomial distribution. Now the
probability of succes ($\pi$) corresponds to the probability that a
SAMPLE MEAN of sample size 10 will be greater than 95 mg/dl. Hopefully,
you see that we calculated this exact probability in part c! Be careful
here, the n for the Binomial distribution is how many samples we take,
not the size of the samples.

$$
X \sim Binom(n = 5, \pi = 0.024) \\
p(x) = \frac{n!}{x!(n-x)!}\pi^x(1-\pi)^{n-x} \\
P(X \geq 1) = 1 - P(X \leq 0) = 1 - P(X = 0) = 1-p(0) \\
1 - p(0) = 1- (\frac{5!}{0!(5)!}0.024^0(1-0.024)^{5}) \\
1-p(0) = 1 - (1-0.024)^5 = 1 - 0.886 = 0.114
$$

Theh probability that AT LEAST ONE of the SAMPLE MEANS will be greater
than 95 mg/dl is 0.114.

```{r}
1-dbinom(0, 5, p)
```

</details>

<br>

#### 3. In the city of Chicago, 235 baseball fans were sampled and asked if they cheer for the White Sox or the Cubs. 155 of the questioned people preferred the Cubs, while the remaining 80 fans preferred the White Sox. Use normal approximation to test the hypothesis that the Cubs and White Sox have an equal number of fans in the city. Construct a 95% Confidence Interval for the true proportion of Chicago fans who cheer for the Cubs, still using normal approximation.

<details>

<summary>Answer</summary>

We have binary/proportion data and we want to use a normal
approximation. The null hypothesis is that each team has equal number of
fans in the city, this corresponds to a proportion of 0.5.

$$
\pi_0 = 0.5 \\
n = 235 \\
\hat\pi = \frac{155}{235} = 0.6595745 \\
SE = \sqrt{\frac{\pi_0(1-\pi_0)}{n}} = \sqrt{\frac{0.5(1-0.5)}{235}} = 0.0326164
$$

$$
z = \frac{\hat\pi - \pi_0}{SE}  = \frac{0.6595745 - 0.5} {0.0326164} = 4.89 \\
p = 2* P(Z > 4.89) = 2 * P(Z < -4.89) \approx 0.000
$$

```{r}
# Writing down our known values
n <- 235
p <- 155 / n
p_0 <- 0.5
alpha <- 0.05

# Finding our z-statistic
z <- (p - p_0) / (sqrt(p_0*(1 - p_0) / n))

# Using our z-statistic to get our p-value
pvalue <- 2*(1 - pnorm(z, mean = 0, sd = 1))
pvalue
```

There is overwhelming evidence to suggest that the proportion of Cubs
fans is not equal to the proportion of White Sox fans in Chicago.

Now we use the corresponding confidence interval for proportions.

**NOTE: SE for the confidence interval is different than the SE for the
hypothesis test when working with the normal approximation for
proportions!!**

$$
n = 235 \\
\hat\pi = \frac{155}{235} = 0.6595745\\
SE = \sqrt{\frac{\hat\pi(1-\hat\pi)}{n}} = \sqrt{\frac{0.6596(1-0.6596)}{235}} = 0.03091015 \\
z_{\alpha/2} = z_{0.05/2} =  z_{0.025} = -1.96
$$

95% CI

$$
\hat\pi \pm z_{\alpha/2} * SE \\
0.6595745 \pm (1.96 * 0.03091015) \\
(0.599, 0.720)
$$

```{r}
SE <- sqrt(p*(1-p) / n)
z_alpha_2 <- qnorm(0.025)
p + c(1, -1)*z_alpha_2*SE
```

Our null hypothesis value $\pi_0=0.5$ is not contained in the interval.
This is further confirmation that we can reject the null hypothesis, as it is
not in the range of plausible values.

</details>

<br>

#### 4. A pharmaceutical company is testing the effectiveness of a new drug for lowering blood pressure in people with prehypertension (systolic BP between 120 mmHg and 139 mmHg). They recruit 9 participants and measure their systolic blood pressure before and after taking the drug for one month. The mean systolic blood pressure for participants before taking the drug was 125 mmHg. The mean systolic blood pressure one month after taking the drug was 120.8 mmHg. The standard deviation for the difference (after - before) was 3.4 mmHg. This information is summarized in the table below.

| **Table 2. Systolic Blood Pressures (mmHg) of prehypertension patients (n = 9)**

|      | Before drug | After drug | Difference |
|------|:-----------:|:----------:|:----------:|
| Mean |     125     |   120.8    |            |
| Sd   |     3.9     |    3.1     |    3.4     |

##### a) Test the hypothesis that the drug has no effect on lowering systolic blood pressure and state your conclusion.

<details>

<summary>Answer</summary>

The data are paired here because each participant's blood pressure is
measured both before and after taking the drug. Since we are given the
means of the two groups and the sample standard deviation of the
difference, we can use the t-test on the difference (use the formula
just like we usually do but make $\bar{x}$ the difference in means). The
difference in means (after - before) is -4.2 mmHg. The null hypothesis
is that the drug has no effect on lowering systolic blood pressure which
is equivalent to $\mu_0 = 0$

$$
\bar{x} = -4.2 \\
\mu_0 = 0 \\
s = 3.4 \\
n = 9 \\
df = 8 \\
SE = \frac{s}{\sqrt{n}} =  \frac{3.4}{\sqrt{9}} = 1.133\\
t = \frac{\bar{x} - \mu_0}{SE} = \frac{-4.2}{1.133} = -3.706
$$

To find the p-value, we can look for the +3.706 in the row with 8
degrees of freedom on the t table. This value falls between 3.36 and
3.83 which corresponds to a p-value between 0.01 and 0.005.

There is strong evidence to suggest that the difference in blood
pressures before and after treatment is not equal to 0 (0.005 \< p \<
0.01). From the data we observed, we can conclude that the treatment is
helpful in lowering blood pressure.

```{r}
n <- 9
s <- 3.4
SE <- 3.4 / sqrt(n)
(t <- (-4.2 - 0) / SE)
(p_value <- 2*pt(t, 8))
```

</details>

<br>

##### b) Calculate the 95% confidence interval for the mean difference in blood pressure before and after taking the drug.

<details>

<summary>Answer</summary>

We have all the parts we need except for $t_{\alpha/2, n-1}$. To find
this using the table, go to the df = 8 row and $\alpha = 0.05$ column.
We don't look up $\alpha = 0.025$ because the table Patrick provides
accounts for both tails but the value we get still corresponds to the
0.025 quantile. When we do that, we get the corresponding t-value is
2.31.

$$
\bar{x} = -4.2 \\
s = 3.4 \\
n = 9 \\
df = 8 \\
SE = \frac{s}{\sqrt{n}} =  \frac{3.4}{\sqrt{9}} = 1.133\\
t_{0.025, 8} = 2.31 \\
95\% CI: \bar{x} \pm t_{0.025, 8}*SE \\
-4.2 \pm 2.31 * 1.133 \\
(-6.82, -1.58)
$$

```{r}
t_0.025 <- qt(0.975, 8)
-4.2 + c(-1,1)*t_0.025*SE
```

</details>

<br>

#### 5. An article in the New England Journal of Medicine reported that among adults living in the United States, the levels of albumin in cerebrospinal fluid is normally distributed with a mean of 29.5 mg/dl, and a standard deviation of 9.25 mg/dl. We are going to select a sample of size 20.

<!-- -->

##### a) How does the variability of our sample mean compare with the variability of albumin levels in the population?

<details>

<summary>Answer</summary>

The standard deviation of the sample mean is the standard error.

$SE = \frac{\sigma}{\sqrt{n}} = \frac{9.25}{\sqrt{20}} = 2.07$

This is LESS than the standard deviation of the albumin levels (9.25)

</details>

<br>

##### b) What is the probability that our sample mean will lie between 29 and 31 mg/dl?

<details>

<summary>Answer</summary>

When we are looking at the probability that our SAMPLE mean will lie
betwen 2 values, we need to use the SE. Since we can assume the data are
normally distributed we can still use the Z formula.

$$
P(29 < \bar{X} < 31) = P(\frac{29 - 29.5}{9.25/\sqrt{20}} < Z < \frac{31 - 29.5}{9.25/\sqrt{20}}) \\
= P(-0.242 < Z < 0.725) \\
= P(Z < 0.725) - P(Z < -0.242)\\
= 0.766 - 0.404 = 0.362
$$

</details>

<br>

##### c) What two values will contain the middle 50% of our sample means?

<details>

<summary>Answer</summary>

For the middle 50% of the sample means, we want to find $z_{0.25}$ and
$z_{0.75}$ in the table. These 2 values are $\pm0.675$. Since we are
finding the middle 50% of the sample means rather than the individuals,
we will use the standard error rather than the standard deviation.
Therefore,

$\mu \pm z *SE$ = $29.5 \pm 0.675 * \frac{9.25}{\sqrt{20}}$ =
$(28.104, 30.896)$

</details>

<br>

##### d) Now assume we don't know the standard deviation or mean of the population and our sample of 20 has a mean albumin level of 30.1 mg/dl and a standard deviation of 8.95 mg/dl. Construct a 95% confidence interval for the population mean and interpret it.

<details>

<summary>Answer</summary>

Since we no longer know the standard deviation of the population, we
will use the T distribution:

$95\% \space CI = \bar{x} \pm t_{0.025, \space 19}{\frac{s}{\sqrt{n}}}$
= $30.1 \pm 2.09*\frac{8.95}{\sqrt{20}}$ = $(25.92, 34.28)$\

(note that if the course website table were used, if 19 degrees of
freedom were on the table, 2.09 would be under $\alpha = 0.05$)

We are 95% confident that the true mean level of albumnin in
cerebrospinal fluid in adults living in the US is between 25.92 mg/dl
and 34.28 mg/dl.

</details>

<br>

##### e) Why is the confidence interval constructed in (d) wider than the confidence interval that would be constructed if we used the normal distribution? What variables affect the width of the confidence interval?

<details>

<summary>Answer</summary>

Three values affect the width of the interval: the standard deviation,
the sample size, and the critical value. In this case, since we are
using the t-distribution, the tails are fatter and the critical value is
larger, causing the interval to be wider.

</details>

<br>

#### 6. In the following scenarios, identify what will happen to the power of a hypothesis test:

<!-- -->

##### a) We increase the sample size.

<details>

<summary>Answer</summary>

Power increases

</details>

<br>

##### b) The standard deviation of the sample is larger than what we expected.

<details>

<summary>Answer</summary>

Power decreases

</details>

<br>

##### c) Our effect size moves from 5 units to 10 units.

<details>

<summary>Answer</summary>

Power increases

</details>

<br>
