--- title: 'STAT 3202: Practice 07' author: "Autumn 2018, OSU" date: '' output: html_document: theme: simplex pdf_document: default urlcolor: BrickRed --- *** # Exercise 1 Suppose that a researcher is interested in the effect of caffeine on [typing speed](https://en.wikipedia.org/wiki/Words_per_minute). A group of nine individuals are administered a typing test. The following day, they repeat the typing test, this time after taking 400 mg of caffeine. (Note: This is not recommended.) The data gathered, measured in words per minute, is ```{r} decaf = c(98, 124, 107, 105, 80, 43, 73, 68, 69) caff = c(104, 128, 110, 108, 86, 53, 72, 73, 72) ``` ```{r, echo = FALSE} library(kableExtra) kable(data.frame(decaf, caff)) %>% kable_styling(bootstrap_options = "striped", full_width = FALSE) ``` Note that these are paired observations. Use the **sign test** with a significance level of 0.05 to assess whether or not caffeine has an effect on typing speed. That is, test $$ H_0\colon \ m_D = m_C - m_N = 0 \quad \text{vs} \quad H_A\colon \ m_D = m_C - m_N \neq 0 $$ where - $m_C$ is the median typing speed in words per minute of individuals using caffeine - $m_N$ is the median typing speed in words per minute of individuals not using caffeine Since it is possible that the caffeine makes typing speed worse, use a two-sided test. You may use the following probabilities calculated in `R`. ```{r} round(dbinom(x = 0:9, size = 9, prob = 0.5), 3) ``` Report: - The **p-value** of the test - A **decision** when $\alpha = 0.05$. ### Solution ```{r} # the "test statistic" for the sign test sum(caff - decaf > 0) # the expected value of the test stat under the null # this is used to determine "extreme" values of the test statistic length(caff - decaf) / 2 # add up the probabilities for test stat values that are as extreme or more extreme sum(dbinom(c(0, 1, 8, 9), size = 9, prob = 0.5)) ``` - **P-Value:** `r sum(dbinom(c(0, 1, 8, 9), size = 9, prob = 0.5))` - **Decision:** Reject $H_0$. *** # Exercise 2 Does meditation have an effect on [blood pressure](https://en.wikipedia.org/wiki/Blood_pressure). A group of six college aged individuals were given a routine physical examination including a measurement of their [systolic](https://en.wikipedia.org/wiki/Systole) blood pressure. (Measured in millimeters of mercury.) A week after their physicals, the same six individuals returned for a guided [meditation session](https://en.wikipedia.org/wiki/Guided_meditation). Immediately afterwords there (systolic) blood pressure was measured. The data gathered is ```{r} physical = c(125, 108, 185, 135, 112, 133) meditation = c(120, 114, 160, 131, 124, 125) ``` ```{r, echo = FALSE} library(kableExtra) kable(data.frame(physical, meditation)) %>% kable_styling(bootstrap_options = "striped", full_width = FALSE) ``` Note that these are paired observations. Use the **sign test** with a significance level of 0.10 to assess whether or not meditation has an effect on blood pressure. That is, test $$ H_0\colon \ m_D = m_M - m_P = 0 \quad \text{vs} \quad H_A\colon \ m_D = m_M - m_P \neq 0 $$ where - $m_P$ is the median systolic blood pressure in millimeters of mercury measured without meditation - $m_M$ is the median systolic blood pressure in millimeters of mercury measured with meditation Since it is possible that the meditation makes blood pressure worse, use a two-sided test. You may use the following probabilities calculated in `R`. ```{r} round(dbinom(x = 0:6, size = 6, prob = 0.5), 3) ``` Report: - The **p-value** of the test - A **decision** when $\alpha = 0.10$. ### Solution ```{r} # the "test statistic" for the sign test sum(meditation - physical > 0) # the expected value of the test stat under the null # this is used to determine "extreme" values of the test statistic length(meditation - physical) / 2 # add up the probabilities for test stat values that are as extreme or more extreme sum(dbinom(c(0, 1, 2, 4, 5, 6), size = 6, prob = 0.5)) ``` - **P-Value:** `r sum(dbinom(c(0, 1, 2, 4, 5, 6), size = 6, prob = 0.5))` - **Decision:** Fail to reject $H_0$. *** # Exercise 3 Return to the sleep data in Exercise 2. This time test - $H_0$: The distribution of systolic blood pressure is **the same** with and without meditation - $H_A$: The distribution of systolic blood pressure is **different** with and without meditation To do so, use a **permutation test** that permutes the *statistic* $$ \bar{x}_D $$ where $\bar{x}_D$ is the sample mean difference. Assume that the distribution of blood pressure with and without meditation has the same shape, but may have different locations. Use at least 10000 permutations. ```{r} physical = c(125, 108, 185, 135, 112, 133) meditation = c(120, 114, 160, 131, 124, 125) ``` - Create a histogram that illustrates the distribution of the statistic used. - Report the p-value of the test. ### Solution ```{r} # create difference data bp_diff = meditation - physical # function to shuffle data and calculate statistic permute_x_bar = function(data) { sample_size = length(data) permuted_data = sample(c(-1, 1), size = sample_size, replace = TRUE) * data mean(permuted_data) } # generate t statistics for sleep data set.seed(42) bp_x_bars = replicate(n = 10000, permute_x_bar(data = bp_diff)) # calculate t statistic on observed data bp_x_bar_obs = mean(bp_diff) ``` ```{r} hist(bp_x_bars, col = "darkgrey", xlab = "t", probability = TRUE, main = "Permutation Test, Sample Mean, Blood Pressure Data") box() grid() abline(v = c(-1, 1) * bp_x_bar_obs, col = "firebrick", lwd = 2) ``` ```{r} mean(bp_x_bars > abs(bp_x_bar_obs)) + mean(bp_x_bars < -abs(bp_x_bar_obs)) ``` - This histogram looks real weird... This is mostly due to the choice of statistic to permute. Had we used some like the usual $t$ statistic, the histogram would be more normal looking. *** # Example 4 Which profession pays more? Data Scientist of Actuary? A (far too small) survey of junior (less than three years experience) data scientist and actuaries resulted in the following data: ```{r} data_sci = c(88000, 121000, 91000, 50000, 78000, 95000) actuary = c(63000, 75000, 81000, 75000, 85000) ``` Use a **permutation test** that permutes the *statistic* $$ t = \frac{(\bar{x} - \bar{y}) - 0}{s_p\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}} $$ to test - $H_0$: The distribution of salaries is **the same** for junior data scientists and actuaries - $H_A$: The distribution of salaries is **different** for junior data scientists and actuaries Assume that the distribution of salaries for both has the same shape, but may have different locations. Use at least 10000 permutations. - Create a histogram that illustrates the distribution of the statistic used. - Report the p-value of the test. ### Solution ```{r} # function to shuffle data and calculate statistic permute_two_t_stat = function(data_1, data_2) { # determine samples sizes of both groups sample_size_1 = length(data_1) sample_size_2 = length(data_2) # create variable for group structure groups = c(rep(TRUE, sample_size_1), rep(FALSE, sample_size_2)) # shuffle the groups shuffled_groups = sample(groups) # merge the data into a single group (null hypothesis) all_data = c(data_1, data_2) # create new groups shuffled_data_1 = all_data[shuffled_groups] shuffled_data_2 = all_data[!shuffled_groups] # calculate statistics on permuted data t.test(x = shuffled_data_1, y = shuffled_data_2, var.equal = TRUE)$statistic } # generate t statistics for exam data set.seed(42) salary_t_stats = replicate(n = 10000, permute_two_t_stat(data_1 = data_sci, data_2 = actuary)) # calculate t statistic on observed data salary_t_obs = t.test(x = data_sci, y = actuary, var.equal = TRUE)$statistic ``` ```{r} hist(salary_t_stats, col = "darkgrey", xlab = "t", probability = TRUE, main = "Permutation t-Test, Salary Data") box() grid() abline(v = c(-1, 1) * salary_t_obs, col = "firebrick", lwd = 2) ``` ```{r} mean(salary_t_stats > abs(salary_t_obs)) + mean(salary_t_stats < -abs(salary_t_obs)) ``` *** # Exercise 5 Repeat exercise 3, but use an appropriate test available in the `R` function `wilcox.test()`. Report: - The **p-value** of the test - A **decision** when $\alpha = 0.05$. ### Solution ```{r} wilcox.test(x = meditation, y = physical, paired = TRUE) ``` - **P-Value:** `r wilcox.test(x = meditation, y = physical, paired = TRUE)$p.value` - **Decision:** Fail to reject $H_0$. Here we used the Wilcoxon **signed rank** test as this is paired data. ***