---
-title: "Week 7 R Lecture"
-author: "Jeremy Foote"
-date: "April 4, 2019"
+title: "Week 6 R lecture"
+subtitle: "Statistics and statistical programming \nNorthwestern University \nMTS 525"
+author: "Aaron Shaw"
+date: "May 3, 2019"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
+# Using built-in functions to conduct hypothesis tests
-## Categorical Data
-
-The goal of this script is to help you think about analyzing categorical data, including proportions, tables, chi-squared tests, and simulation.
-
-### Estimating proportions
-
-If a survey of 50 randomly sampled Chicagoans found that 45% of them thought that Giordano's made the best deep dish pizza, what would be the 95% confidence interval for the true proportion of Chicagoans who prefer Giordano's?
-
-Can we reject the hypothesis that 50% of Chicagoans prefer Giordano's?
+## T-tests
+You learned the theory/concepts behind t-tests last week, so here's a brief run-down on how to use built-in functions in R to conduct them and interpret the results. As a reminder, t-tests are used to estimate the difference-in-means between two distributions. You might be working with a one-sample test, two (independent samples, pooled or paired samples, or difference scores. R can handle most of this stuff and, if not, remember that you know how to write your own functions!
+The usual procedure for conducting t-tests in R is pretty straightforward: you use the `t.test()` function. For my example, I'll simulate two groups by randomly sampling the `rivers` dataset. I will then conduct a two-sample t-test for difference of means under a null hypothesis of no difference.
```{r}
-est = .45
-sample_size = 50
-SE = sqrt(est*(1-est)/sample_size)
+set.seed(20190503)
+group.a <- sample(rivers, 35, replace=TRUE)
+group.b <- sample(rivers, 42, replace=TRUE)
-conf_int = c(est - 1.96 * SE, est + 1.96 * SE)
-conf_int
+t.test(group.a, group.b)
```
+Notice everything that is contained in the output: It tells me what test I just ran (a two sample t-test). It also tells me the names of the two groups. Then it reports the test statistic (t) value, the degrees of freedom, and a p-value. It states my alternative hypothesis explicitly and includes a 95% confidence interval before also reporting the sample means for each group.
-What if we had the same result but had sampled 500 people?
-
+Go read the documentation for `t.test()` in R as the function has many possible arguments to help you conduct exactly the sort of t-test you want given the particulars of your data, hypotheses, etc. The documentation will also tell you quite a lot about the output or "values" returned by the function. As with most tests and operations in R, you can store the output of a `t.test()` call as an object. In the code below, notice what values the object has and how you might use the object to call specific values later.
```{r}
-est = .45
-sample_size = 500
-SE = sqrt(est*(1-est)/sample_size)
-
-conf_int = c(est - 1.96 * SE, est + 1.96 * SE)
-conf_int
-```
-
-### Tabular Data
-
-The Iris dataset is composed of measurements of flower dimensions. It comes packaged with R and is often used in examples. Here we make a table of how often each species in the dataset has a sepal width greater than 3.
+t.result.ab <- t.test(group.a, group.b)
-```{r}
+t.result.ab ## same as before
-table(iris$Species, iris$Sepal.Width > 3)
+names(t.result.ab)
+t.result.ab$conf.int
```
+You never need to calculate a confidence interval for a difference in means by hand again. You can even set other $\alpha$ values in the call to t-test if you want to estimate another interval (see the documentation).
+
+## ANOVAs
+Here's a very brief introduction to how ANOVAs work in R. Let's use the `iris` dataset (look up `help(iris)` to learn about it) to build an example in which we test whether irises of different species have different average sepal width. Remember that the purpose of an ANOVA is to compare some continuous outcome across two or more categories to test the global hypothesis of any differences between groups.
-The chi-squared test is a test of how much the frequencies we see in a table differ from what we would expect if there was no difference between the groups.
+In the code block below, I explicitly store the `iris` dataset and then call `aov()` to conduct the test. In the call to `aov()` I indicate the variables and dataset I want to use. Notice the use of the tilde character ("~"). This is formula notation and tells R that I'm using the variables before and after as the dimensions of my ANOVA. The `aov()` function follows a convention that the dependent or outcome variable has to go first in the formula (just try it the other way if you want to see what happens — since ANOVAs assume a continuous outcome it won't work to try and use the categorical measure as the DV).
```{r}
+iris <- iris
-chisq.test(table(iris$Species, iris$Sepal.Width > 3))
+aov(Sepal.Width ~ Species, data=iris)
```
+Wait a minute. That doesn't look quite right. Where's the F statistic? Where's the p-value? Why doesn't this table look like the ANOVA tables the book showed us?
-The incredibly low p-value means that it is very unlikely that these came from the same distribution and that sepal width differs by species.
-
-
-
-## Using Simulation
-
-When the assumptions of Chi-squared tests aren't met, we can use simulation to approximate how likely a given result is.
-
-The book uses the example of a medical practitioner who has 3 complications out of 62 procedures, while the typical rate is 10%.
-
-The null hypothesis is that this practitioner's true rate is also 10%, so we're trying to figure out how rare it would be to have 3 or fewer complications, if the true rate is 10%.
-
+Not to worry. As with many statistical tests in R, the `aov()` function did everything it promised, you just need to ask it explicitly to show you the "summary" of the test to get the information you're looking for. You can also store the results and call `summary()` as a separate command.
```{r}
-# We write a function that we are going to replicate
-simulation <- function(rate = .1, n = 62){
- # Draw n random numbers from a uniform distribution from 0 to 1
- draws = runif(n)
- # If rate = .4, on average, .4 of the draws will be less than .4
- # So, we consider those draws where the value is less than `rate` as complications
- complication_count = sum(draws < rate)
- # Then, we return the total count
- return(complication_count)
-}
+summary(aov(Sepal.Width ~ Species, data=iris))
-# The replicate function runs a function many times
-
-simulated_complications <- replicate(5000, simulation())
+aov1.result <- aov(Sepal.Width ~ Species, data=iris)
+summary(aov1.result)
```
+That should look more familiar. It's good to get accustomed to storing your model results and then calling `summary()` on the output. Such a routine is typical for many other types of statistical tests in R
-We can look at our simulated complications
+Note that the object returned by `aov()` has values that you can also call directly just like we did with the output of `t.test()` above. In the example below, I've asked R to show me the residuals (deviation from the sample mean across all groups) for each row of the dataset. I'm not sure what you might use this for, but it's available.
```{r}
+names(aov1.result)
-hist(simulated_complications)
-```
-
-And determine how many of them are as extreme or more extreme than the value we saw. This is the p-value.
+aov1.result$residuals
-```{r}
-
-sum(simulated_complications <= 3)/length(simulated_complications)
```
+As with all functions in R, I *strongly* encourage you to read the documentation before you set out to conduct your own ANOVAs in the wild. There are details and settings that it may be important to consider. In addition, you might notice in the documentation for ANOVAs that R actually conducts the ANOVA by performing a linear regression (whaaaat??!?!?). Turns out that ANOVA is just a special case of linear regression with some idiosyncratic assumptions and standards for the type of summary information that needs to be reported. No need to worry about the details of that for now, but it can be good to have in the back of your mind for future reference in case you find yourself in a situation where you had planned to run an ANOVA, but your data just doesn't fit the assumptions necessary to identify the model.
+
projects / stats_class_2019.git / blobdiff