---
title: "ADA1: Class 15, Hypothesis testing (one- and two-sample)"
author: anonymous
date: "10/11/2016"
output:
html_document:
toc: true
---
Include your answers in this document in the sections below the rubric.
# Rubric
Answer the questions with the two data examples.
---
# Mechanics of a hypothesis test (review)
1. Set up the __null and alternative hypotheses__ in words and notation.
* In words: ``The population mean for [what is being studied] is different from [value of $\mu_0$].''
(Note that the statement in words is in terms of the alternative hypothesis.)
* In notation: $H_0: \mu=\mu_0$ versus $H_A: \mu \ne \mu_0$
(where $\mu_0$ is specified by the context of the problem).
2. Choose the __significance level__ of the test, such as $\alpha=0.05$.
3. Compute the __test statistic__, such as $t_{s} = \frac{\bar{Y}-\mu_0}{SE_{\bar{Y}}}$, where $SE_{\bar{Y}}=s/\sqrt{n}$ is the standard error.
4. Determine the __tail(s)__ of the sampling distribution where the __$p$-value__ from the test statistic will be calculated
(for example, both tails, right tail, or left tail).
(Historically, we would compare the observed test statistic, $t_{s}$,
with the __critical value__ $t_{\textrm{crit}}=t_{\alpha/2}$
in the direction of the alternative hypothesis from the
$t$-distribution table with degrees of freedom $df = n-1$.)
5. State the __conclusion__ in terms of the problem.
* Reject $H_0$ in favor of $H_A$ if $p\textrm{-value} < \alpha$.
* Fail to reject $H_0$ if $p\textrm{-value} \ge \alpha$.
(Note: We DO NOT _accept_ $H_0$.)
6. __Check assumptions__ of the test (next week).
---
# Height data for our class
Is the population mean height of UNM students eligible to take Stat 427/527
different from the [US average](https://en.wikipedia.org/wiki/Human_height#Average_height_around_the_world)
for men (5 ft 9 1/2 in) or women (5 ft 4 in)?
```{R, cache = TRUE}
# install.packages("gsheet")
# Height vs Hand Span
library(gsheet)
dat.hand.url <- "docs.google.com/spreadsheets/d/1D22hva-Em1ZMsLpmo5yNJpuwzAHScZbphN281wunw9I"
dat.hand <- gsheet2tbl(dat.hand.url)
dat.hand <- as.data.frame(dat.hand)
dat.hand <- na.omit(dat.hand)
dat.hand$Gender_M_F <- factor(dat.hand$Gender_M_F, levels = c("F", "M"))
str(dat.hand)
```
Plot the estimated mean from our class sample versus the true US mean.
```{R}
## If we create a summary data.frame with a similar structure as our data, then we
## can annotate our plot with those summaries.
# calculate the estimated mean and order M then F
est.mean <- as.numeric(by(dat.hand$Height_in, dat.hand$Gender_M_F, mean))
# combine true US mean with our estimated mean
height.true <- data.frame(Gender_M_F = unique(dat.hand$Gender_M_F)
, Height_in = c(64, 69.5, est.mean)
, TrueEst = c(rep("True", 2), rep("Est", 2)))
height.true
```
Here's two ways to plot our data, annotating the observed and hypothesized means.
```{R}
#$
library(ggplot2)
p <- ggplot(data = dat.hand, aes(x = Gender_M_F, y = Height_in))
p <- p + geom_boxplot(alpha = 1/4)
p <- p + geom_jitter(position = position_jitter(width = 0.1))
p <- p + geom_point(data = height.true, aes(colour = TrueEst, shape = TrueEst), size = 4, alpha = 3/4)
print(p)
library(ggplot2)
p <- ggplot(data = dat.hand, aes(x = Height_in))
p <- p + geom_histogram(binwidth = 1)
p <- p + geom_vline(data = height.true, aes(xintercept = Height_in, colour = TrueEst, linetype = TrueEst))
p <- p + facet_grid(Gender_M_F ~ .)
print(p)
```
## Conduct the hypothesis tests
```{R}
# look at help for t.test
# ?t.test
# defaults include: alternative = "two.sided", conf.level = 0.95
```
_(I've hidden some code from Chapter 02 defining a function to plot the t-distribution with shaded p-value.)_
```{R, echo=FALSE}
## From ADA1 Chapter 02
# Function to plot t-distribution with shaded p-value
t.dist.pval <- function(t.summary) {
par(mfrow=c(1,1))
lim.extreme <- max(4, abs(t.summary$statistic) + 0.5)
lim.lower <- -lim.extreme;
lim.upper <- lim.extreme;
x.curve <- seq(lim.lower, lim.upper, length=200)
y.curve <- dt(x.curve, df = t.summary$parameter)
plot(x.curve, y.curve, type = "n"
, ylab = paste("t-dist( df =", signif(t.summary$parameter, 3), ")")
, xlab = paste("t-stat =", signif(t.summary$statistic, 5)
, ", Shaded area is p-value =", signif(t.summary$p.value, 5)))
if ((t.summary$alternative == "less")
| (t.summary$alternative == "two.sided")) {
x.pval.l <- seq(lim.lower, -abs(t.summary$statistic), length=200)
y.pval.l <- dt(x.pval.l, df = t.summary$parameter)
polygon(c(lim.lower, x.pval.l, -abs(t.summary$statistic))
, c(0, y.pval.l, 0), col="gray")
}
if ((t.summary$alternative == "greater")
| (t.summary$alternative == "two.sided")) {
x.pval.u <- seq(abs(t.summary$statistic), lim.upper, length=200)
y.pval.u <- dt(x.pval.u, df = t.summary$parameter)
polygon(c(abs(t.summary$statistic), x.pval.u, lim.upper)
, c(0, y.pval.u, 0), col="gray")
}
points(x.curve, y.curve, type = "l", lwd = 2, col = "black")
}
```
### Test female height equal to US (model example).
```{R, fig.width=6, fig.height=4}
# test females
t.summary.F <- t.test(subset(dat.hand, Gender_M_F == "F", Height_in)
, mu = 64)
t.summary.F
names(t.summary.F)
t.dist.pval(t.summary.F)
```
__Hypothesis test__
1. ``The population mean height for females at UNM eligible to take Stat 427/527 is different from the US population value of $\mu_0=64$ inches.''
* $H_0: \mu=64$ versus $H_A: \mu \ne 64$
2. Let $\alpha=0.05$, the significance level of the test and the Type-I error probability if the null hypothesis is true.
3. $t_{s} = `r signif(t.summary.F$statistic, 4)`$.
4. $p=`r signif(t.summary.F$p.value, 3)`$, this is the observed significance of the test.
5. Because $p=`r signif(t.summary.F$p.value, 3)` < 0.05$,
we have sufficient evidence to reject $H_0$, concluding that the
observed mean height is different than the US population mean.
## Questions to answer
1. (3 p) As above, set up the hypothesis test for males, but whether UNM males are taller on average than the US population.
### Test male height greater than US (your turn).
```{R, fig.width=6, fig.height=4}
## You'll need to modify the statement below to correspond
## to the hypothesis you wish to test
# test males
t.summary.M <- t.test(subset(dat.hand, Gender_M_F == "M", Height_in)
, mu = 0
, alternative = "two.sided")
t.summary.M
t.dist.pval(t.summary.M)
```
__Hypothesis test__
1. ``''
* $H_0: $ versus $H_A: $
2. Let $\alpha=0.05$, the significance level of the test and the Type-I error probability if the null hypothesis is true.
3. $t_{s} = $.
4. $p = $, this is the observed significance of the test.
5. Because $p = $, ...
Also, given your conclusion, state whether you could have made a Type-I or Type-II error.
---
# Earth's water example
## Questions to answer
1. (3 p) In class we sampled the beach ball and observed 20 of 30 observations were water.
Conduct a hypothesis test to determine whether the proportion of water on the beach ball
is different from the amount of water on the earth's surface (71%).
```{R}
## notes for prop.test() and binom.test()
# x = number of "successes"
# n = total sample size
#n = 2
#x = 1
n = 30
x = 20
dat.globe <- data.frame(type = c("Water", "Land"), freq = c(x, n - x), prop = c(x, n - x) / n)
dat.globe
# binom.test() is an exact test for a binomial random variable
b.summary <- binom.test(x = x, n = n, p = 0.5, conf.level = 0.95)
b.summary
```
```{R, fig.width=6, fig.height=2}
library(ggplot2)
p <- ggplot(data = subset(dat.globe, type == "Water"), aes(x = type, y = prop))
p <- p + geom_hline(yintercept = c(0, 1), alpha = 1/4)
p <- p + geom_bar(stat = "identity")
p <- p + geom_errorbar(aes(min = b.summary$conf.int[1], max = b.summary$conf.int[2]), width=0.25)
p <- p + geom_hline(yintercept = 0.71, colour = "red")
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + coord_flip() # flip the x and y axes for horizontal plot
print(p)
```
__Hypothesis test__
1. ``''
* $H_0: $ versus $H_A: $
2. Let $\alpha=0.05$, the significance level of the test and the Type-I error probability if the null hypothesis is true.
3. $t_{s} = $.
4. $p = $, this is the observed significance of the test.
5. Because $p = $, ...
Also, given your conclusion, state whether you could have made a Type-I or Type-II error.
---
# African countries in the UN example
Previously in class we collected data using a randomized experiment.
We provided a priming number (X = 10 or 65, not actually a random number) then asked you two questions:
1. Do you think the percentage of countries, among all those in the United
Nations, that are in Africa is higher or lower than X?
2. Give your best estimate of the percentage of countries, among all those in
the United Nations, that are in Africa.
The data were compiled into a google doc which we read in below.
```{R, cache = TRUE}
# install.packages("gsheet")
# Height vs Hand Span
library(gsheet)
dat.UN.Africa.url <- "docs.google.com/spreadsheets/d/1PkjIygZdtlopH6kvyle-KeoECsycHWyCaya92erlvaE"
dat.UN.Africa <- gsheet2tbl(dat.UN.Africa.url)
dat.UN.Africa <- as.data.frame(dat.UN.Africa)
dat.UN.Africa <- na.omit(dat.UN.Africa)
dat.UN.Africa$PrimingNumber <- factor(dat.UN.Africa$PrimingNumber)
dat.UN.Africa$HighLow <- factor(dat.UN.Africa$HighLow)
str(dat.UN.Africa)
```
Here are some summaries and plots.
```{R}
## If we create a summary data.frame with a similar structure as our data, then we
## can annotate our plot with those summaries.
# calculate the estimated mean and order M then F
est.mean <- as.numeric(by(dat.UN.Africa$UN_Percentage, dat.UN.Africa$PrimingNumber, mean))
# combine true US mean with our estimated mean
mean.UN.Africa <- data.frame(PrimingNumber = unique(dat.UN.Africa$PrimingNumber)
, UN_Percentage = est.mean)
mean.UN.Africa
# histogram using ggplot
p <- ggplot(dat.UN.Africa, aes(x = UN_Percentage))
p <- p + geom_histogram(binwidth = 4)
p <- p + geom_rug()
p <- p + geom_vline(data = mean.UN.Africa, aes(xintercept = UN_Percentage, colour = "red"))
p <- p + facet_grid(PrimingNumber ~ .)
print(p)
p <- ggplot(dat.UN.Africa, aes(x = UN_Percentage, fill=PrimingNumber))
p <- p + geom_histogram(binwidth = 4, alpha = 0.5, position="identity")
p <- p + geom_rug()
p <- p + geom_vline(data = mean.UN.Africa, aes(xintercept = UN_Percentage, colour = PrimingNumber, linetype = PrimingNumber))
p <- p + geom_rug(aes(colour = PrimingNumber), alpha = 1/2)
print(p)
```
## Questions to answer
_A priori_, before we observed the data, we hypothesized that those who were primed
with a larger number (65) would provide a higher percentage (`UN_Percentage`)
than those with the lower number (10).
Therefore, this is a one-sided test.
1. (4 p) Set up the two-sample t-test and state the conclusions.
```{R, fig.width=6, fig.height=4}
# two-sample t-test
t.summary.UN <- t.test(UN_Percentage ~ PrimingNumber, data = dat.UN.Africa
, alternative = "less")
t.summary.UN
t.dist.pval(t.summary.UN)
```
__Hypothesis test__
1. ``''
* $H_0: $ versus $H_A: $
2. Let $\alpha=0.05$, the significance level of the test and the Type-I error probability if the null hypothesis is true.
3. $t_{s} = $.
4. $p = $, this is the observed significance of the test.
5. Because $p = $, ...
Also, given your conclusion, state whether you could have made a Type-I or Type-II error.