Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Statistical inference is the process of drawing formal conclusions from -data.
- -In our class, we wil define formal statistical inference as settings where one wants to infer facts about a population using noisy -statistical data where uncertainty must be accounted for.
- -In every major election, pollsters would like to know, ahead of the -actual election, who's going to win. Here, the target of -estimation (the estimand) is clear, the percentage of people in -a particular group (city, state, county, country or other electoral -grouping) who will vote for each candidate.
- -We can not poll everyone. Even if we could, some polled -may change their vote by the time the election occurs. -How do we collect a reasonable subset of data and quantify the -uncertainty in the process to produce a good guess at who will win?
- -A large clinical trial (the Women’s Health Initiative) published results in 2002 that contradicted prior evidence on the efficacy of hormone replacement therapy for post menopausal women and suggested a negative impact of HRT for several key health outcomes. Based on a statistically based protocol, the study was stopped early due an excess number of negative events.
- -Here's there's two inferential problems.
- -See WHI writing group paper JAMA 2002, Vol 288:321 - 333. for the paper and Steinkellner et al. Menopause 2012, Vol 19:616 621 for adiscussion of the long term impacts
- -In 1985 a group at a major neonatal intensive care center published the results of a trial comparing a standard treatment and a promising new extracorporeal membrane oxygenation treatment (ECMO) for newborn infants with severe respiratory failure. Ethical considerations lead to a statistical randomization scheme whereby one infant received the control therapy, thereby opening the study to sample-size based criticisms.
- -For a review and statistical discussion, see Royall Statistical Science 1991, Vol 6, No. 1, 52-88
- -We won't spend too much time talking about this, but there are several different -styles of inference. Two broad categories that get discussed a lot are:
- -Data scientists tend to fall within shades of gray of these and various other schools of inference.
- -Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Statistical inference is the process of drawing formal conclusions from +data.
+ +In our class, we wil define formal statistical inference as settings where one wants to infer facts about a population using noisy +statistical data where uncertainty must be accounted for.
+ +In every major election, pollsters would like to know, ahead of the +actual election, who's going to win. Here, the target of +estimation (the estimand) is clear, the percentage of people in +a particular group (city, state, county, country or other electoral +grouping) who will vote for each candidate.
+ +We can not poll everyone. Even if we could, some polled +may change their vote by the time the election occurs. +How do we collect a reasonable subset of data and quantify the +uncertainty in the process to produce a good guess at who will win?
+ +A large clinical trial (the Women’s Health Initiative) published results in 2002 that contradicted prior evidence on the efficacy of hormone replacement therapy for post menopausal women and suggested a negative impact of HRT for several key health outcomes. Based on a statistically based protocol, the study was stopped early due an excess number of negative events.
+ +Here's there's two inferential problems.
+ +See WHI writing group paper JAMA 2002, Vol 288:321 - 333. for the paper and Steinkellner et al. Menopause 2012, Vol 19:616 621 for adiscussion of the long term impacts
+ +In 1985 a group at a major neonatal intensive care center published the results of a trial comparing a standard treatment and a promising new extracorporeal membrane oxygenation treatment (ECMO) for newborn infants with severe respiratory failure. Ethical considerations lead to a statistical randomization scheme whereby one infant received the control therapy, thereby opening the study to sample-size based criticisms.
+ +For a review and statistical discussion, see Royall Statistical Science 1991, Vol 6, No. 1, 52-88
+ +We won't spend too much time talking about this, but there are several different +styles of inference. Two broad categories that get discussed a lot are:
+ +Data scientists tend to fall within shades of gray of these and various other schools of inference.
+ +Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Normal set operations have particular interpretations in this setting
- -A probability measure, \(P\), is a function from the collection of possible events so that the following hold
- -Part 3 of the definition implies finite additivity
- -\[ -P(\cup_{i=1}^n A_i) = \sum_{i=1}^n P(A_i) -\] -where the \(\{A_i\}\) are mutually exclusive. (Note a more general version of -additivity is used in advanced classes.)
- -The National Sleep Foundation (www.sleepfoundation.org) reports that around 3% of the American population has sleep apnea. They also report that around 10% of the North American and European population has restless leg syndrome. Does this imply that 13% of people will have at least one sleep problems of these sorts?
- -Answer: No, the events are not mutually exclusive. To elaborate let:
- -\[ -\begin{eqnarray*} - A_1 & = & \{\mbox{Person has sleep apnea}\} \\ - A_2 & = & \{\mbox{Person has RLS}\} - \end{eqnarray*} -\]
- -Then
- -\[ -\begin{eqnarray*} - P(A_1 \cup A_2 ) & = & P(A_1) + P(A_2) - P(A_1 \cap A_2) \\ - & = & 0.13 - \mbox{Probability of having both} - \end{eqnarray*} -\] -Likely, some fraction of the population has both.
- -A probability mass function evaluated at a value corresponds to the -probability that a random variable takes that value. To be a valid -pmf a function, \(p\), must satisfy
- -The sum is taken over all of the possible values for \(x\).
- -Let \(X\) be the result of a coin flip where \(X=0\) represents -tails and \(X = 1\) represents heads. -\[ -p(x) = (1/2)^{x} (1/2)^{1-x} ~~\mbox{ for }~~x = 0,1 -\] -Suppose that we do not know whether or not the coin is fair; Let -\(\theta\) be the probability of a head expressed as a proportion -(between 0 and 1). -\[ -p(x) = \theta^{x} (1 - \theta)^{1-x} ~~\mbox{ for }~~x = 0,1 -\]
- -A probability density function (pdf), is a function associated with -a continuous random variable
- -Areas under pdfs correspond to probabilities for that random variable
- -To be a valid pdf, a function \(f\) must satisfy
- -\(f(x) \geq 0\) for all \(x\)
The area under \(f(x)\) is one.
Suppose that the proportion of help calls that get addressed in -a random day by a help line is given by -\[ -f(x) = \left\{\begin{array}{ll} - 2 x & \mbox{ for } 1 > x > 0 \\ - 0 & \mbox{ otherwise} -\end{array} \right. -\]
- -Is this a mathematically valid density?
- -x <- c(-0.5, 0, 1, 1, 1.5)
-y <- c(0, 0, 2, 0, 0)
-plot(x, y, lwd = 3, frame = FALSE, type = "l")
-
What is the probability that 75% or fewer of calls get addressed?
- -
1.5 * 0.75/2
-## [1] 0.5625
-pbeta(0.75, 2, 1)
-## [1] 0.5625
-What are the survival function and CDF from the density considered before?
- -For \(1 \geq x \geq 0\) -\[ -F(x) = P(X \leq x) = \frac{1}{2} Base \times Height = \frac{1}{2} (x) \times (2 x) = x^2 -\]
- -\[ -S(x) = 1 - x^2 -\]
- -pbeta(c(0.4, 0.5, 0.6), 2, 1)
-## [1] 0.16 0.25 0.36
-sqrt(0.5)
-## [1] 0.7071
-qbeta(0.5, 2, 1)
-## [1] 0.7071
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Normal set operations have particular interpretations in this setting
+ +A probability measure, \(P\), is a function from the collection of possible events so that the following hold
+ +Part 3 of the definition implies finite additivity
+ +\[ +P(\cup_{i=1}^n A_i) = \sum_{i=1}^n P(A_i) +\] +where the \(\{A_i\}\) are mutually exclusive. (Note a more general version of +additivity is used in advanced classes.)
+ +The National Sleep Foundation (www.sleepfoundation.org) reports that around 3% of the American population has sleep apnea. They also report that around 10% of the North American and European population has restless leg syndrome. Does this imply that 13% of people will have at least one sleep problems of these sorts?
+ +Answer: No, the events are not mutually exclusive. To elaborate let:
+ +\[ +\begin{eqnarray*} + A_1 & = & \{\mbox{Person has sleep apnea}\} \\ + A_2 & = & \{\mbox{Person has RLS}\} + \end{eqnarray*} +\]
+ +Then
+ +\[ +\begin{eqnarray*} + P(A_1 \cup A_2 ) & = & P(A_1) + P(A_2) - P(A_1 \cap A_2) \\ + & = & 0.13 - \mbox{Probability of having both} + \end{eqnarray*} +\] +Likely, some fraction of the population has both.
+ +A probability mass function evaluated at a value corresponds to the +probability that a random variable takes that value. To be a valid +pmf a function, \(p\), must satisfy
+ +The sum is taken over all of the possible values for \(x\).
+ +Let \(X\) be the result of a coin flip where \(X=0\) represents +tails and \(X = 1\) represents heads. +\[ +p(x) = (1/2)^{x} (1/2)^{1-x} ~~\mbox{ for }~~x = 0,1 +\] +Suppose that we do not know whether or not the coin is fair; Let +\(\theta\) be the probability of a head expressed as a proportion +(between 0 and 1). +\[ +p(x) = \theta^{x} (1 - \theta)^{1-x} ~~\mbox{ for }~~x = 0,1 +\]
+ +A probability density function (pdf), is a function associated with +a continuous random variable
+ +Areas under pdfs correspond to probabilities for that random variable
+ +To be a valid pdf, a function \(f\) must satisfy
+ +\(f(x) \geq 0\) for all \(x\)
The area under \(f(x)\) is one.
Suppose that the proportion of help calls that get addressed in +a random day by a help line is given by +\[ +f(x) = \left\{\begin{array}{ll} + 2 x & \mbox{ for } 1 > x > 0 \\ + 0 & \mbox{ otherwise} +\end{array} \right. +\]
+ +Is this a mathematically valid density?
+ +x <- c(-0.5, 0, 1, 1, 1.5)
+y <- c(0, 0, 2, 0, 0)
+plot(x, y, lwd = 3, frame = FALSE, type = "l")
+
What is the probability that 75% or fewer of calls get addressed?
+ +
1.5 * 0.75/2
+## [1] 0.5625
+pbeta(0.75, 2, 1)
+## [1] 0.5625
+What are the survival function and CDF from the density considered before?
+ +For \(1 \geq x \geq 0\) +\[ +F(x) = P(X \leq x) = \frac{1}{2} Base \times Height = \frac{1}{2} (x) \times (2 x) = x^2 +\]
+ +\[ +S(x) = 1 - x^2 +\]
+ +pbeta(c(0.4, 0.5, 0.6), 2, 1)
+## [1] 0.16 0.25 0.36
+sqrt(0.5)
+## [1] 0.7071
+qbeta(0.5, 2, 1)
+## [1] 0.7071
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
library(manipulate)
-myHist <- function(mu){
- hist(galton$child,col="blue",breaks=100)
- lines(c(mu, mu), c(0, 150),col="red",lwd=5)
- mse <- mean((galton$child - mu)^2)
- text(63, 150, paste("mu = ", mu))
- text(63, 140, paste("Imbalance = ", round(mse, 2)))
-}
-manipulate(myHist(mu), mu = slider(62, 74, step = 0.5))
-hist(galton$child, col = "blue", breaks = 100)
-meanChild <- mean(galton$child)
-lines(rep(meanChild, 100), seq(0, 150, length = 100), col = "red", lwd = 5)
-
\[ -Var(X) = E[(X - \mu)^2] -\]
- -the expected (squared) distance from the mean
- -What's the sample variance from the result of a toss of a die?
- -\(Var(X) = E[X^2] - E[X]^2 \approx 2.92\)
What's the sample variance from the result of the toss of a coin with probability of heads (1) of \(p\)?
- -\(Var(X) = E[X^2] - E[X]^2 = p - p^2 = p(1 - p)\)
Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
## Loading required package: MASS
+
library(manipulate)
+myHist <- function(mu){
+ hist(galton$child,col="blue",breaks=100)
+ lines(c(mu, mu), c(0, 150),col="red",lwd=5)
+ mse <- mean((galton$child - mu)^2)
+ text(63, 150, paste("mu = ", mu))
+ text(63, 140, paste("Imbalance = ", round(mse, 2)))
+}
+manipulate(myHist(mu), mu = slider(62, 74, step = 0.5))
+hist(galton$child, col = "blue", breaks = 100)
+meanChild <- mean(galton$child)
+lines(rep(meanChild, 100), seq(0, 150, length = 100), col = "red", lwd = 5)
+
\[ +Var(X) = E[(X - \mu)^2] +\]
+ +the expected (squared) distance from the mean
+ +What's the sample variance from the result of a toss of a die?
+ +\(Var(X) = E[X^2] - E[X]^2 \approx 2.92\)
What's the sample variance from the result of the toss of a coin with probability of heads (1) of \(p\)?
+ +\(Var(X) = E[X^2] - E[X]^2 = p - p^2 = p(1 - p)\)
+
+
+---
+## Using manipulate
+```
+library(manipulate)
+myHist <- function(mu){
+ hist(galton$child,col="blue",breaks=100)
+ lines(c(mu, mu), c(0, 150),col="red",lwd=5)
+ mse <- mean((galton$child - mu)^2)
+ text(63, 150, paste("mu = ", mu))
+ text(63, 140, paste("Imbalance = ", round(mse, 2)))
+}
+manipulate(myHist(mu), mu = slider(62, 74, step = 0.5))
+```
+
+---
+## The center of mass is the empirical mean
+
+```r
+hist(galton$child, col = "blue", breaks = 100)
+meanChild <- mean(galton$child)
+lines(rep(meanChild, 100), seq(0, 150, length = 100), col = "red", lwd = 5)
+```
+
+
+
+
+---
+## Example
+
+- Suppose a coin is flipped and $X$ is declared $0$ or $1$ corresponding to a head or a tail, respectively
+- What is the expected value of $X$?
+ $$
+ E[X] = .5 \times 0 + .5 \times 1 = .5
+ $$
+- Note, if thought about geometrically, this answer is obvious; if two equal weights are spaced at 0 and 1, the center of mass will be $.5$
+
+
+
+---
+
+## Example
+
+- Suppose that a die is rolled and $X$ is the number face up
+- What is the expected value of $X$?
+ $$
+ E[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} +
+ 3 \times \frac{1}{6} + 4 \times \frac{1}{6} +
+ 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5
+ $$
+- Again, the geometric argument makes this answer obvious without calculation.
+
+---
+
+## Continuous random variables
+
+- For a continuous random variable, $X$, with density, $f$, the expected
+ value is defined as follows
+ $$
+ E[X] = \mbox{the area under the function}~~~ t f(t)
+ $$
+- This definition borrows from the definition of center of mass for a continuous body
+
+---
+
+## Example
+
+- Consider a density where $f(x) = 1$ for $x$ between zero and one
+- (Is this a valid density?)
+- Suppose that $X$ follows this density; what is its expected value?
+
+
+
+---
+
+## Rules about expected values
+
+- The expected value is a linear operator
+- If $a$ and $b$ are not random and $X$ and $Y$ are two random variables then
+ - $E[aX + b] = a E[X] + b$
+ - $E[X + Y] = E[X] + E[Y]$
+
+---
+
+## Example
+
+- You flip a coin, $X$ and simulate a uniform random number $Y$, what is the expected value of their sum?
+ $$
+ E[X + Y] = E[X] + E[Y] = .5 + .5 = 1
+ $$
+- Another example, you roll a die twice. What is the expected value of the average?
+- Let $X_1$ and $X_2$ be the results of the two rolls
+ $$
+ E[(X_1 + X_2) / 2] = \frac{1}{2}(E[X_1] + E[X_2])
+ = \frac{1}{2}(3.5 + 3.5) = 3.5
+ $$
+
+---
+
+## Example
+
+1. Let $X_i$ for $i=1,\ldots,n$ be a collection of random variables, each from a distribution with mean $\mu$
+2. Calculate the expected value of the sample average of the $X_i$
+$$
+ \begin{eqnarray*}
+ E\left[ \frac{1}{n}\sum_{i=1}^n X_i\right]
+ & = & \frac{1}{n} E\left[\sum_{i=1}^n X_i\right] \\
+ & = & \frac{1}{n} \sum_{i=1}^n E\left[X_i\right] \\
+ & = & \frac{1}{n} \sum_{i=1}^n \mu = \mu.
+ \end{eqnarray*}
+$$
+
+---
+
+## Remark
+
+- Therefore, the expected value of the **sample mean** is the population mean that it's trying to estimate
+- When the expected value of an estimator is what its trying to estimate, we say that the estimator is **unbiased**
+
+---
+
+## The variance
+
+- The variance of a random variable is a measure of *spread*
+- If $X$ is a random variable with mean $\mu$, the variance of $X$ is defined as
+
+$$
+Var(X) = E[(X - \mu)^2]
+$$
+
+the expected (squared) distance from the mean
+- Densities with a higher variance are more spread out than densities with a lower variance
+
+---
+
+- Convenient computational form
+$$
+Var(X) = E[X^2] - E[X]^2
+$$
+- If $a$ is constant then $Var(aX) = a^2 Var(X)$
+- The square root of the variance is called the **standard deviation**
+- The standard deviation has the same units as $X$
+
+---
+
+## Example
+
+- What's the sample variance from the result of a toss of a die?
+
+ - $E[X] = 3.5$
+ - $E[X^2] = 1 ^ 2 \times \frac{1}{6} + 2 ^ 2 \times \frac{1}{6} + 3 ^ 2 \times \frac{1}{6} + 4 ^ 2 \times \frac{1}{6} + 5 ^ 2 \times \frac{1}{6} + 6 ^ 2 \times \frac{1}{6} = 15.17$
+
+- $Var(X) = E[X^2] - E[X]^2 \approx 2.92$
+
+---
+
+## Example
+
+- What's the sample variance from the result of the toss of a coin with probability of heads (1) of $p$?
+
+ - $E[X] = 0 \times (1 - p) + 1 \times p = p$
+ - $E[X^2] = E[X] = p$
+
+- $Var(X) = E[X^2] - E[X]^2 = p - p^2 = p(1 - p)$
+
+---
+
+## Interpreting variances
+
+- Chebyshev's inequality is useful for interpreting variances
+- This inequality states that
+$$
+P(|X - \mu| \geq k\sigma) \leq \frac{1}{k^2}
+$$
+- For example, the probability that a random variable lies beyond $k$ standard deviations from its mean is less than $1/k^2$
+$$
+\begin{eqnarray*}
+ 2\sigma & \rightarrow & 25\% \\
+ 3\sigma & \rightarrow & 11\% \\
+ 4\sigma & \rightarrow & 6\%
+\end{eqnarray*}
+$$
+- Note this is only a bound; the actual probability might be quite a bit smaller
+
+---
+
+## Example
+
+- IQs are often said to be distributed with a mean of $100$ and a sd of $15$
+- What is the probability of a randomly drawn person having an IQ higher than $160$ or below $40$?
+- Thus we want to know the probability of a person being more than $4$ standard deviations from the mean
+- Thus Chebyshev's inequality suggests that this will be no larger than 6\%
+- IQs distributions are often cited as being bell shaped, in which case this bound is very conservative
+- The probability of a random draw from a bell curve being $4$ standard deviations from the mean is on the order of $10^{-5}$ (one thousandth of one percent)
+
+---
+
+## Example
+
+- A former buzz phrase in industrial quality control is Motorola's "Six Sigma" whereby businesses are suggested to control extreme events or rare defective parts
+- Chebyshev's inequality states that the probability of a "Six Sigma" event is less than $1/6^2 \approx 3\%$
+- If a bell curve is assumed, the probability of a "six sigma" event is on the order of $10^{-9}$ (one ten millionth of a percent)
+
diff --git a/06_StatisticalInference/01_03_Expectations/index.pdf b/06_StatisticalInference/01_03_Expectations/index.pdf
index c9c43b5a3..aa71b7bd4 100644
Binary files a/06_StatisticalInference/01_03_Expectations/index.pdf and b/06_StatisticalInference/01_03_Expectations/index.pdf differ
diff --git a/06_StatisticalInference/01_04_Independence/assets/fig/unnamed-chunk-2.png b/06_StatisticalInference/01_04_Independence/assets/fig/unnamed-chunk-2.png
new file mode 100644
index 000000000..eded2a301
Binary files /dev/null and b/06_StatisticalInference/01_04_Independence/assets/fig/unnamed-chunk-2.png differ
diff --git a/06_StatisticalInference/01_04_Independence/index.html b/06_StatisticalInference/01_04_Independence/index.html
index 57c94e8b1..3ce8e98d5 100644
--- a/06_StatisticalInference/01_04_Independence/index.html
+++ b/06_StatisticalInference/01_04_Independence/index.html
@@ -1,496 +1,496 @@
-
-
-
- Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
If \(A\) is independent of \(B\) then
- -We will use the following fact extensively in this class:
- -If a collection of random variables \(X_1, X_2, \ldots, X_n\) are independent, then their joint distribution is the product of their individual densities or mass functions
- -That is, if \(f_i\) is the density for random variable \(X_i\) we have that -\[ -f(x_1,\ldots, x_n) = \prod_{i=1}^n f_i(x_i) -\]
- -Let \(\{X_i\}_{i=1}^n\) be a collection of random variables
- -A commonly used subcase from these properties is that if a collection of random variables \(\{X_i\}\) are uncorrelated, then the variance of the sum is the sum of the variances -\[ -Var\left(\sum_{i=1}^n X_i \right) = \sum_{i=1}^n Var(X_i) -\]
Therefore, it is sums of variances that tend to be useful, not sums of standard deviations; that is, the standard deviation of the sum of bunch of independent random variables is the square root of the sum of the variances, not the sum of the standard deviations
Suppose \(X_i\) are iid with variance \(\sigma^2\)
- -\[ -\begin{eqnarray*} - Var(\bar X) & = & Var \left( \frac{1}{n}\sum_{i=1}^n X_i \right)\\ \\ - & = & \frac{1}{n^2} Var\left(\sum_{i=1}^n X_i \right)\\ \\ - & = & \frac{1}{n^2} \sum_{i=1}^n Var(X_i) \\ \\ - & = & \frac{1}{n^2} \times n\sigma^2 \\ \\ - & = & \frac{\sigma^2}{n} - \end{eqnarray*} -\]
- -\[ - \begin{eqnarray*} - E\left[\sum_{i=1}^n (X_i - \bar X)^2\right] & = & \sum_{i=1}^n E\left[X_i^2\right] - n E\left[\bar X^2\right] \\ \\ - & = & \sum_{i=1}^n \left\{Var(X_i) + \mu^2\right\} - n \left\{Var(\bar X) + \mu^2\right\} \\ \\ - & = & \sum_{i=1}^n \left\{\sigma^2 + \mu^2\right\} - n \left\{\sigma^2 / n + \mu^2\right\} \\ \\ - & = & n \sigma^2 + n \mu ^ 2 - \sigma^2 - n \mu^2 \\ \\ - & = & (n - 1) \sigma^2 - \end{eqnarray*} -\]
- -data(father.son)
-x <- father.son$sheight
-n <- length(x)
-
round(c(sum((x - mean(x))^2)/(n - 1), var(x), var(x)/n, sd(x), sd(x)/sqrt(n)),
- 2)
-## [1] 7.92 7.92 0.01 2.81 0.09
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
If \(A\) is independent of \(B\) then
+ +We will use the following fact extensively in this class:
+ +If a collection of random variables \(X_1, X_2, \ldots, X_n\) are independent, then their joint distribution is the product of their individual densities or mass functions
+ +That is, if \(f_i\) is the density for random variable \(X_i\) we have that +\[ +f(x_1,\ldots, x_n) = \prod_{i=1}^n f_i(x_i) +\]
+ +Let \(\{X_i\}_{i=1}^n\) be a collection of random variables
+ +A commonly used subcase from these properties is that if a collection of random variables \(\{X_i\}\) are uncorrelated, then the variance of the sum is the sum of the variances +\[ +Var\left(\sum_{i=1}^n X_i \right) = \sum_{i=1}^n Var(X_i) +\]
Therefore, it is sums of variances that tend to be useful, not sums of standard deviations; that is, the standard deviation of the sum of bunch of independent random variables is the square root of the sum of the variances, not the sum of the standard deviations
Suppose \(X_i\) are iid with variance \(\sigma^2\)
+ +\[ +\begin{eqnarray*} + Var(\bar X) & = & Var \left( \frac{1}{n}\sum_{i=1}^n X_i \right)\\ \\ + & = & \frac{1}{n^2} Var\left(\sum_{i=1}^n X_i \right)\\ \\ + & = & \frac{1}{n^2} \sum_{i=1}^n Var(X_i) \\ \\ + & = & \frac{1}{n^2} \times n\sigma^2 \\ \\ + & = & \frac{\sigma^2}{n} + \end{eqnarray*} +\]
+ +\[ + \begin{eqnarray*} + E\left[\sum_{i=1}^n (X_i - \bar X)^2\right] & = & \sum_{i=1}^n E\left[X_i^2\right] - n E\left[\bar X^2\right] \\ \\ + & = & \sum_{i=1}^n \left\{Var(X_i) + \mu^2\right\} - n \left\{Var(\bar X) + \mu^2\right\} \\ \\ + & = & \sum_{i=1}^n \left\{\sigma^2 + \mu^2\right\} - n \left\{\sigma^2 / n + \mu^2\right\} \\ \\ + & = & n \sigma^2 + n \mu ^ 2 - \sigma^2 - n \mu^2 \\ \\ + & = & (n - 1) \sigma^2 + \end{eqnarray*} +\]
+ +data(father.son)
+x <- father.son$sheight
+n <- length(x)
+
round(c(sum((x - mean(x))^2)/(n - 1), var(x), var(x)/n, sd(x), sd(x)/sqrt(n)),
+ 2)
+## [1] 7.92 7.92 0.01 2.81 0.09
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
\[ -P(B ~|~ A) = \frac{P(A ~|~ B) P(B)}{P(A ~|~ B) P(B) + P(A ~|~ B^c)P(B^c)}. -\]
- -\[ -\begin{eqnarray*} - P(D ~|~ +) & = &\frac{P(+~|~D)P(D)}{P(+~|~D)P(D) + P(+~|~D^c)P(D^c)}\\ \\ - & = & \frac{P(+~|~D)P(D)}{P(+~|~D)P(D) + \{1-P(-~|~D^c)\}\{1 - P(D)\}} \\ \\ - & = & \frac{.997\times .001}{.997 \times .001 + .015 \times .999}\\ \\ - & = & .062 -\end{eqnarray*} -\]
- -Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
\[ +P(B ~|~ A) = \frac{P(A ~|~ B) P(B)}{P(A ~|~ B) P(B) + P(A ~|~ B^c)P(B^c)}. +\]
+ +\[ +\begin{eqnarray*} + P(D ~|~ +) & = &\frac{P(+~|~D)P(D)}{P(+~|~D)P(D) + P(+~|~D^c)P(D^c)}\\ \\ + & = & \frac{P(+~|~D)P(D)}{P(+~|~D)P(D) + \{1-P(-~|~D^c)\}\{1 - P(D)\}} \\ \\ + & = & \frac{.997\times .001}{.997 \times .001 + .015 \times .999}\\ \\ + & = & .062 +\end{eqnarray*} +\]
+ +Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
n <- 5
-pvals <- seq(0, 1, length = 1000)
-plot(c(0, 1), c(0, 1.2), type = "n", frame = FALSE, xlab = "p", ylab = "likelihood")
-text((0 : n) /n, 1.1, as.character(0 : n))
-sapply(0 : n, function(x) {
- phat <- x / n
- if (x == 0) lines(pvals, ( (1 - pvals) / (1 - phat) )^(n-x), lwd = 3)
- else if (x == n) lines(pvals, (pvals / phat) ^ x, lwd = 3)
- else lines(pvals, (pvals / phat ) ^ x * ( (1 - pvals) / (1 - phat) ) ^ (n-x), lwd = 3)
- }
-)
-title(paste("Likelihoods for n = ", n))
-
\[\left( - \begin{array}{c} - n \\ 0 - \end{array} - \right) = -\left( - \begin{array}{c} - n \\ n - \end{array} - \right) = 1 - \]
- -choose(8, 7) * .5 ^ 8 + choose(8, 8) * .5 ^ 8
-[1] 0.03516
-pbinom(6, size = 8, prob = .5, lower.tail = FALSE)
-[1] 0.03516
-plot(pvals, dbinom(7, 8, pvals) / dbinom(7, 8, 7/8) ,
- lwd = 3, frame = FALSE, type = "l", xlab = "p", ylab = "likelihood")
-
zvals <- seq(-3, 3, length = 1000)
-plot(zvals, dnorm(zvals),
- type = "l", lwd = 3, frame = FALSE, xlab = "z", ylab = "Density")
-sapply(-3 : 3, function(k) abline(v = k))
-
[[1]]
-NULL
-
-[[2]]
-NULL
-
-[[3]]
-NULL
-
-[[4]]
-NULL
-
-[[5]]
-NULL
-
-[[6]]
-NULL
-
-[[7]]
-NULL
-qnorm(.95, mean = mu, sd = sd)The number of people that show up at a bus stop is Poisson with -a mean of \(2.5\) per hour.
- -If watching the bus stop for 4 hours, what is the probability that \(3\) -or fewer people show up for the whole time?
- -ppois(3, lambda = 2.5 * 4)
-[1] 0.01034
-We flip a coin with success probablity \(0.01\) five hundred times.
- -What's the probability of 2 or fewer successes?
- -pbinom(2, size = 500, prob = .01)
-[1] 0.1234
-ppois(2, lambda=500 * .01)
-[1] 0.1247
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
n <- 5
+pvals <- seq(0, 1, length = 1000)
+plot(c(0, 1), c(0, 1.2), type = "n", frame = FALSE, xlab = "p", ylab = "likelihood")
+text((0 : n) /n, 1.1, as.character(0 : n))
+sapply(0 : n, function(x) {
+ phat <- x / n
+ if (x == 0) lines(pvals, ( (1 - pvals) / (1 - phat) )^(n-x), lwd = 3)
+ else if (x == n) lines(pvals, (pvals / phat) ^ x, lwd = 3)
+ else lines(pvals, (pvals / phat ) ^ x * ( (1 - pvals) / (1 - phat) ) ^ (n-x), lwd = 3)
+ }
+)
+title(paste("Likelihoods for n = ", n))
+
\[\left( + \begin{array}{c} + n \\ 0 + \end{array} + \right) = +\left( + \begin{array}{c} + n \\ n + \end{array} + \right) = 1 + \]
+ +choose(8, 7) * 0.5^8 + choose(8, 8) * 0.5^8
+## [1] 0.03516
+pbinom(6, size = 8, prob = 0.5, lower.tail = FALSE)
+## [1] 0.03516
+plot(pvals, dbinom(7, 8, pvals)/dbinom(7, 8, 7/8), lwd = 3, frame = FALSE, type = "l",
+ xlab = "p", ylab = "likelihood")
+
zvals <- seq(-3, 3, length = 1000)
+plot(zvals, dnorm(zvals), type = "l", lwd = 3, frame = FALSE, xlab = "z", ylab = "Density")
+sapply(-3:3, function(k) abline(v = k))
+
qnorm(.95, mean = mu, sd = sd)The number of people that show up at a bus stop is Poisson with +a mean of \(2.5\) per hour.
+ +If watching the bus stop for 4 hours, what is the probability that \(3\) +or fewer people show up for the whole time?
+ +ppois(3, lambda = 2.5 * 4)
+## [1] 0.01034
+We flip a coin with success probablity \(0.01\) five hundred times.
+ +What's the probability of 2 or fewer successes?
+ +pbinom(2, size = 500, prob = 0.01)
+## [1] 0.1234
+ppois(2, lambda = 500 * 0.01)
+## [1] 0.1247
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Imagine a sequence
- -Clearly this sequence converges to \(1\)
Definition of a limit: For any fixed distance we can find a point in the sequence so that the sequence is closer to the limit than that distance from that point on
Consider \(\bar X_n\) the sample average of the first \(n\) of a collection of \(iid\) observations
- -We say that \(\bar X_n\) converges in probability to a limit if for any fixed distance the probability of \(\bar X_n\) being closer (further away) than that distance from the limit converges to one (zero)
n <- 10000; means <- cumsum(rnorm(n)) / (1 : n)
-plot(1 : n, means, type = "l", lwd = 2,
- frame = FALSE, ylab = "cumulative means", xlab = "sample size")
-abline(h = 0)
-
in Galton's data
- -library(UsingR);data(father.son); x <- father.son$sheight
-(mean(x) + c(-1, 1) * qnorm(.975) * sd(x) / sqrt(length(x))) / 12
-[1] 5.710 5.738
-1/sqrt(100)=.1 so a back of the envelope calculation gives an approximate 95% interval of (0.46, 0.66)
-
-round(1 / sqrt(10 ^ (1 : 6)), 3)
-[1] 0.316 0.100 0.032 0.010 0.003 0.001
-x <- 5; t <- 94.32; lambda <- x / t
-round(lambda + c(-1, 1) * qnorm(.975) * sqrt(lambda / t), 3)
-[1] 0.007 0.099
-poisson.test(x, T = 94.32)$conf
-[1] 0.01721 0.12371
-attr(,"conf.level")
-[1] 0.95
-exp(confint(glm(x ~ 1 + offset(log(t)), family = poisson(link = log))))
- 2.5 % 97.5 %
-0.01901 0.11393
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Imagine a sequence
+ +Clearly this sequence converges to \(1\)
Definition of a limit: For any fixed distance we can find a point in the sequence so that the sequence is closer to the limit than that distance from that point on
Consider \(\bar X_n\) the sample average of the first \(n\) of a collection of \(iid\) observations
+ +We say that \(\bar X_n\) converges in probability to a limit if for any fixed distance the probability of \(\bar X_n\) being closer (further away) than that distance from the limit converges to one (zero)
n <- 10000
+means <- cumsum(rnorm(n))/(1:n)
+plot(1:n, means, type = "l", lwd = 2, frame = FALSE, ylab = "cumulative means",
+ xlab = "sample size")
+abline(h = 0)
+
in Galton's data
+ +library(UsingR)
+data(father.son)
+x <- father.son$sheight
+(mean(x) + c(-1, 1) * qnorm(0.975) * sd(x)/sqrt(length(x)))/12
+## [1] 5.710 5.738
+1/sqrt(100)=.1 so a back of the envelope calculation gives an approximate 95% interval of (0.46, 0.66)
+
+round(1/sqrt(10^(1:6)), 3)
+## [1] 0.316 0.100 0.032 0.010 0.003 0.001
+x <- 5
+t <- 94.32
+lambda <- x/t
+round(lambda + c(-1, 1) * qnorm(0.975) * sqrt(lambda/t), 3)
+## [1] 0.007 0.099
+poisson.test(x, T = 94.32)$conf
+## [1] 0.01721 0.12371
+## attr(,"conf.level")
+## [1] 0.95
+exp(confint(glm(x ~ 1 + offset(log(t)), family = poisson(link = log))))
+## Waiting for profiling to be done...
+## 2.5 % 97.5 %
+## 0.01901 0.11393
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Throughout we use the following general procedure for creating CIs
- -a. Create a Pivot or statistic that does not depend on the parameter of interest
- -b. Solve the probability that the pivot lies between bounds for the parameter
Note that if \(\chi^2_{n-1, \alpha}\) is the \(\alpha\) quantile of the -Chi-squared distribution then
- -\[ -\begin{eqnarray*} - 1 - \alpha & = & P \left( \chi^2_{n-1, \alpha/2} \leq \frac{(n - 1) S^2}{\sigma^2} \leq \chi^2_{n-1,1 - \alpha/2} \right) \\ \\ -& = & P\left(\frac{(n-1)S^2}{\chi^2_{n-1,1-\alpha/2}} \leq \sigma^2 \leq -\frac{(n-1)S^2}{\chi^2_{n-1,\alpha/2}} \right) \\ -\end{eqnarray*} -\] -So that -\[ -\left[\frac{(n-1)S^2}{\chi^2_{n-1,1-\alpha/2}}, \frac{(n-1)S^2}{\chi^2_{n-1,\alpha/2}}\right] -\] -is a \(100(1-\alpha)\%\) confidence interval for \(\sigma^2\)
- -library(UsingR)
-data(father.son)
-x <- father.son$sheight
-s <- sd(x)
-n <- length(x)
-round(sqrt((n - 1) * s^2/qchisq(c(0.975, 0.025), n - 1)), 3)
-## [1] 2.701 2.939
-Suppose that \((X_1,\ldots,X_n)\) are iid \(N(\mu,\sigma^2)\), then: -a. \(\frac{\bar X - \mu}{\sigma / \sqrt{n}}\) is standard normal -b. \(\sqrt{\frac{(n - 1) S^2}{\sigma^2 (n - 1)}} = S / \sigma\) is the square root of a Chi-squared divided by its df
Therefore -\[ -\frac{\frac{\bar X - \mu}{\sigma /\sqrt{n}}}{S/\sigma} -= \frac{\bar X - \mu}{S/\sqrt{n}} -\] -follows Gosset's \(t\) distribution with \(n-1\) degrees of freedom
In R typing data(sleep) brings up the sleep data originally
-analyzed in Gosset's Biometrika paper, which shows the increase in
-hours for 10 patients on two soporific drugs. R treats the data as two
-groups rather than paired.
data(sleep)
-head(sleep)
-## extra group ID
-## 1 0.7 1 1
-## 2 -1.6 1 2
-## 3 -0.2 1 3
-## 4 -1.2 1 4
-## 5 -0.1 1 5
-## 6 3.4 1 6
-g1 <- sleep$extra[1:10]
-g2 <- sleep$extra[11:20]
-difference <- g2 - g1
-mn <- mean(difference)
-s <- sd(difference)
-n <- 10
-mn + c(-1, 1) * qt(0.975, n - 1) * s/sqrt(n)
-## [1] 0.7001 2.4599
-t.test(difference)$conf.int
-## [1] 0.7001 2.4599
-## attr(,"conf.level")
-## [1] 0.95
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Throughout we use the following general procedure for creating CIs
+ +a. Create a Pivot or statistic that does not depend on the parameter of interest
+ +b. Solve the probability that the pivot lies between bounds for the parameter
Note that if \(\chi^2_{n-1, \alpha}\) is the \(\alpha\) quantile of the +Chi-squared distribution then
+ +\[ +\begin{eqnarray*} + 1 - \alpha & = & P \left( \chi^2_{n-1, \alpha/2} \leq \frac{(n - 1) S^2}{\sigma^2} \leq \chi^2_{n-1,1 - \alpha/2} \right) \\ \\ +& = & P\left(\frac{(n-1)S^2}{\chi^2_{n-1,1-\alpha/2}} \leq \sigma^2 \leq +\frac{(n-1)S^2}{\chi^2_{n-1,\alpha/2}} \right) \\ +\end{eqnarray*} +\] +So that +\[ +\left[\frac{(n-1)S^2}{\chi^2_{n-1,1-\alpha/2}}, \frac{(n-1)S^2}{\chi^2_{n-1,\alpha/2}}\right] +\] +is a \(100(1-\alpha)\%\) confidence interval for \(\sigma^2\)
+ +library(UsingR)
+data(father.son)
+x <- father.son$sheight
+s <- sd(x)
+n <- length(x)
+round(sqrt((n - 1) * s^2/qchisq(c(0.975, 0.025), n - 1)), 3)
+## [1] 2.701 2.939
+Suppose that \((X_1,\ldots,X_n)\) are iid \(N(\mu,\sigma^2)\), then: +a. \(\frac{\bar X - \mu}{\sigma / \sqrt{n}}\) is standard normal +b. \(\sqrt{\frac{(n - 1) S^2}{\sigma^2 (n - 1)}} = S / \sigma\) is the square root of a Chi-squared divided by its df
Therefore +\[ +\frac{\frac{\bar X - \mu}{\sigma /\sqrt{n}}}{S/\sigma} += \frac{\bar X - \mu}{S/\sqrt{n}} +\] +follows Gosset's \(t\) distribution with \(n-1\) degrees of freedom
In R typing data(sleep) brings up the sleep data originally
+analyzed in Gosset's Biometrika paper, which shows the increase in
+hours for 10 patients on two soporific drugs. R treats the data as two
+groups rather than paired.
data(sleep)
+head(sleep)
+## extra group ID
+## 1 0.7 1 1
+## 2 -1.6 1 2
+## 3 -0.2 1 3
+## 4 -1.2 1 4
+## 5 -0.1 1 5
+## 6 3.4 1 6
+g1 <- sleep$extra[1:10]
+g2 <- sleep$extra[11:20]
+difference <- g2 - g1
+mn <- mean(difference)
+s <- sd(difference)
+n <- 10
+mn + c(-1, 1) * qt(0.975, n - 1) * s/sqrt(n)
+## [1] 0.7001 2.4599
+t.test(difference)$conf.int
+## [1] 0.7001 2.4599
+## attr(,"conf.level")
+## [1] 0.95
+Brian Caffo, Roger Peng, Jeff Leek
Johns Hopkins Bloomberg School of Public Health
Given a statistical probability mass function or density, say \(f(x, \theta)\), where \(\theta\) is an unknown parameter, the likelihood is \(f\) viewed as a function of \(\theta\) for a fixed, observed value of \(x\).
- -The likelihood has the following properties:
- -pvals <- seq(0, 1, length = 1000)
-plot(pvals, dbinom(3, 4, pvals) / dbinom(3, 4, 3/4), type = "l", frame = FALSE, lwd = 3, xlab = "p", ylab = "likelihood / max likelihood")
-lambda <- seq(0, .2, length = 1000)
-likelihood <- dpois(5, 94 * lambda) / dpois(5, 5)
-plot(lambda, likelihood, frame = FALSE, lwd = 3, type = "l", xlab = expression(lambda))
-lines(rep(5/94, 2), 0 : 1, col = "red", lwd = 3)
-lines(range(lambda[likelihood > 1/16]), rep(1/16, 2), lwd = 2)
-lines(range(lambda[likelihood > 1/8]), rep(1/8, 2), lwd = 2)
-Brian Caffo, Roger Peng, Jeff Leek
Johns Hopkins Bloomberg School of Public Health
Given a statistical probability mass function or density, say \(f(x, \theta)\), where \(\theta\) is an unknown parameter, the likelihood is \(f\) viewed as a function of \(\theta\) for a fixed, observed value of \(x\).
+ +The likelihood has the following properties:
+ +pvals <- seq(0, 1, length = 1000)
+plot(pvals, dbinom(3, 4, pvals)/dbinom(3, 4, 3/4), type = "l", frame = FALSE,
+ lwd = 3, xlab = "p", ylab = "likelihood / max likelihood")
+
lambda <- seq(0, 0.2, length = 1000)
+likelihood <- dpois(5, 94 * lambda)/dpois(5, 5)
+plot(lambda, likelihood, frame = FALSE, lwd = 3, type = "l", xlab = expression(lambda))
+lines(rep(5/94, 2), 0:1, col = "red", lwd = 3)
+lines(range(lambda[likelihood > 1/16]), rep(1/16, 2), lwd = 2)
+lines(range(lambda[likelihood > 1/8]), rep(1/8, 2), lwd = 2)
+
Brian Caffo, Roger Peng, Jeff Leek
Johns Hopkins Bloomberg School of Public Health
## Exploring the beta density
-library(manipulate)
-pvals <- seq(0.01, 0.99, length = 1000)
-manipulate(
- plot(pvals, dbeta(pvals, alpha, beta), type = "l", lwd = 3, frame = FALSE),
- alpha = slider(0.01, 10, initial = 1, step = .5),
- beta = slider(0.01, 10, initial = 1, step = .5)
- )
-\[ -\begin{align} -E[p ~|~ X] & = \frac{\tilde \alpha}{\tilde \alpha + \tilde \beta}\\ \\ -& = \frac{x + \alpha}{x + \alpha + n - x + \beta}\\ \\ -& = \frac{x + \alpha}{n + \alpha + \beta} \\ \\ -& = \frac{x}{n} \times \frac{n}{n + \alpha + \beta} + \frac{\alpha}{\alpha + \beta} \times \frac{\alpha + \beta}{n + \alpha + \beta} \\ \\ -& = \mbox{MLE} \times \pi + \mbox{Prior Mean} \times (1 - \pi) -\end{align} -\]
- -pvals <- seq(0.01, 0.99, length = 1000)
-x <- 13; n <- 20
-myPlot <- function(alpha, beta){
- plot(0 : 1, 0 : 1, type = "n", xlab = "p", ylab = "", frame = FALSE)
- lines(pvals, dbeta(pvals, alpha, beta) / max(dbeta(pvals, alpha, beta)),
- lwd = 3, col = "darkred")
- lines(pvals, dbinom(x,n,pvals) / dbinom(x,n,x/n), lwd = 3, col = "darkblue")
- lines(pvals, dbeta(pvals, alpha+x, beta+(n-x)) / max(dbeta(pvals, alpha+x, beta+(n-x))),
- lwd = 3, col = "darkgreen")
- title("red=prior,green=posterior,blue=likelihood")
-}
-manipulate(
- myPlot(alpha, beta),
- alpha = slider(0.01, 10, initial = 1, step = .5),
- beta = slider(0.01, 10, initial = 1, step = .5)
- )
-library(binom)
-binom.bayes(13, 20, type = "highest")
- method x n shape1 shape2 mean lower upper sig
-1 bayes 13 20 13.5 7.5 0.6429 0.4423 0.8361 0.05
-gives the HPD interval.
- -pvals <- seq(0.01, 0.99, length = 1000)
-x <- 13; n <- 20
-myPlot2 <- function(alpha, beta, cl){
- plot(pvals, dbeta(pvals, alpha+x, beta+(n-x)), type = "l", lwd = 3,
- xlab = "p", ylab = "", frame = FALSE)
- out <- binom.bayes(x, n, type = "highest",
- prior.shape1 = alpha,
- prior.shape2 = beta,
- conf.level = cl)
- p1 <- out$lower; p2 <- out$upper
- lines(c(p1, p1, p2, p2), c(0, dbeta(c(p1, p2), alpha+x, beta+(n-x)), 0),
- type = "l", lwd = 3, col = "darkred")
-}
-manipulate(
- myPlot2(alpha, beta, cl),
- alpha = slider(0.01, 10, initial = 1, step = .5),
- beta = slider(0.01, 10, initial = 1, step = .5),
- cl = slider(0.01, 0.99, initial = 0.95, step = .01)
- )
-Brian Caffo, Roger Peng, Jeff Leek
Johns Hopkins Bloomberg School of Public Health
## Exploring the beta density
+library(manipulate)
+pvals <- seq(0.01, 0.99, length = 1000)
+manipulate(
+ plot(pvals, dbeta(pvals, alpha, beta), type = "l", lwd = 3, frame = FALSE),
+ alpha = slider(0.01, 10, initial = 1, step = .5),
+ beta = slider(0.01, 10, initial = 1, step = .5)
+ )
+\[ +\begin{align} +E[p ~|~ X] & = \frac{\tilde \alpha}{\tilde \alpha + \tilde \beta}\\ \\ +& = \frac{x + \alpha}{x + \alpha + n - x + \beta}\\ \\ +& = \frac{x + \alpha}{n + \alpha + \beta} \\ \\ +& = \frac{x}{n} \times \frac{n}{n + \alpha + \beta} + \frac{\alpha}{\alpha + \beta} \times \frac{\alpha + \beta}{n + \alpha + \beta} \\ \\ +& = \mbox{MLE} \times \pi + \mbox{Prior Mean} \times (1 - \pi) +\end{align} +\]
+ +pvals <- seq(0.01, 0.99, length = 1000)
+x <- 13; n <- 20
+myPlot <- function(alpha, beta){
+ plot(0 : 1, 0 : 1, type = "n", xlab = "p", ylab = "", frame = FALSE)
+ lines(pvals, dbeta(pvals, alpha, beta) / max(dbeta(pvals, alpha, beta)),
+ lwd = 3, col = "darkred")
+ lines(pvals, dbinom(x,n,pvals) / dbinom(x,n,x/n), lwd = 3, col = "darkblue")
+ lines(pvals, dbeta(pvals, alpha+x, beta+(n-x)) / max(dbeta(pvals, alpha+x, beta+(n-x))),
+ lwd = 3, col = "darkgreen")
+ title("red=prior,green=posterior,blue=likelihood")
+}
+manipulate(
+ myPlot(alpha, beta),
+ alpha = slider(0.01, 100, initial = 1, step = .5),
+ beta = slider(0.01, 100, initial = 1, step = .5)
+ )
+library(binom)
+## Error: there is no package called 'binom'
+binom.bayes(13, 20, type = "highest")
+## Error: could not find function "binom.bayes"
+gives the HPD interval.
+ +pvals <- seq(0.01, 0.99, length = 1000)
+x <- 13; n <- 20
+myPlot2 <- function(alpha, beta, cl){
+ plot(pvals, dbeta(pvals, alpha+x, beta+(n-x)), type = "l", lwd = 3,
+ xlab = "p", ylab = "", frame = FALSE)
+ out <- binom.bayes(x, n, type = "highest",
+ prior.shape1 = alpha,
+ prior.shape2 = beta,
+ conf.level = cl)
+ p1 <- out$lower; p2 <- out$upper
+ lines(c(p1, p1, p2, p2), c(0, dbeta(c(p1, p2), alpha+x, beta+(n-x)), 0),
+ type = "l", lwd = 3, col = "darkred")
+}
+manipulate(
+ myPlot2(alpha, beta, cl),
+ alpha = slider(0.01, 10, initial = 1, step = .5),
+ beta = slider(0.01, 10, initial = 1, step = .5),
+ cl = slider(0.01, 0.99, initial = 0.95, step = .01)
+ )
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
sp <- sqrt((7 * 15.34^2 + 20 * 18.23^2) / (8 + 21 - 2))
-132.86 - 127.44 + c(-1, 1) * qt(.975, 27) * sp * (1 / 8 + 1 / 21)^.5
-[1] -9.521 20.361
-data(sleep)
-x1 <- sleep$extra[sleep$group == 1]
-x2 <- sleep$extra[sleep$group == 2]
-n1 <- length(x1)
-n2 <- length(x2)
-sp <- sqrt( ((n1 - 1) * sd(x1)^2 + (n2-1) * sd(x2)^2) / (n1 + n2-2))
-md <- mean(x1) - mean(x2)
-semd <- sp * sqrt(1 / n1 + 1/n2)
-md + c(-1, 1) * qt(.975, n1 + n2 - 2) * semd
-[1] -3.3639 0.2039
-t.test(x1, x2, paired = FALSE, var.equal = TRUE)$conf
-[1] -3.3639 0.2039
-attr(,"conf.level")
-[1] 0.95
-t.test(x1, x2, paired = TRUE)$conf
-[1] -2.4599 -0.7001
-attr(,"conf.level")
-[1] 0.95
-t.test(..., var.equal = FALSE)Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
sp <- sqrt((7 * 15.34^2 + 20 * 18.23^2)/(8 + 21 - 2))
+132.86 - 127.44 + c(-1, 1) * qt(0.975, 27) * sp * (1/8 + 1/21)^0.5
+## [1] -9.521 20.361
+data(sleep)
+x1 <- sleep$extra[sleep$group == 1]
+x2 <- sleep$extra[sleep$group == 2]
+n1 <- length(x1)
+n2 <- length(x2)
+sp <- sqrt(((n1 - 1) * sd(x1)^2 + (n2 - 1) * sd(x2)^2)/(n1 + n2 - 2))
+md <- mean(x1) - mean(x2)
+semd <- sp * sqrt(1/n1 + 1/n2)
+md + c(-1, 1) * qt(0.975, n1 + n2 - 2) * semd
+## [1] -3.3639 0.2039
+t.test(x1, x2, paired = FALSE, var.equal = TRUE)$conf
+## [1] -3.3639 0.2039
+## attr(,"conf.level")
+## [1] 0.95
+t.test(x1, x2, paired = TRUE)$conf
+## [1] -2.4599 -0.7001
+## attr(,"conf.level")
+## [1] 0.95
+
t.test(..., var.equal = FALSE)Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
| Truth | -Decide | -Result | -
|---|---|---|
| \(H_0\) | -\(H_0\) | -Correctly accept null | -
| \(H_0\) | -\(H_a\) | -Type I error | -
| \(H_a\) | -\(H_a\) | -Correctly reject null | -
| \(H_a\) | -\(H_0\) | -Type II error | -
\[ -\begin{align} -0.05 & = P\left(\bar X \geq C ~|~ \mu = 30 \right) \\ - & = P\left(\frac{\bar X - 30}{10 / \sqrt{100}} \geq \frac{C - 30}{10/\sqrt{100}} ~|~ \mu = 30\right) \\ - & = P\left(Z \geq \frac{C - 30}{1}\right) \\ -\end{align} -\]
- -library(UsingR); data(father.son)
-t.test(father.son$sheight - father.son$fheight)
-
- One Sample t-test
-
-data: father.son$sheight - father.son$fheight
-t = 11.79, df = 1077, p-value < 2.2e-16
-alternative hypothesis: true mean is not equal to 0
-95 percent confidence interval:
- 0.831 1.163
-sample estimates:
-mean of x
- 0.997
-| Rejection region | -Type I error rate | -
|---|---|
| [0 : 8] | -1 | -
| [1 : 8] | -0.9961 | -
| [2 : 8] | -0.9648 | -
| [3 : 8] | -0.8555 | -
| [4 : 8] | -0.6367 | -
| [5 : 8] | -0.3633 | -
| [6 : 8] | -0.1445 | -
| [7 : 8] | -0.0352 | -
| [8 : 8] | -0.0039 | -
Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
| Truth | +Decide | +Result | +
|---|---|---|
| \(H_0\) | +\(H_0\) | +Correctly accept null | +
| \(H_0\) | +\(H_a\) | +Type I error | +
| \(H_a\) | +\(H_a\) | +Correctly reject null | +
| \(H_a\) | +\(H_0\) | +Type II error | +
\[ +\begin{align} +0.05 & = P\left(\bar X \geq C ~|~ \mu = 30 \right) \\ + & = P\left(\frac{\bar X - 30}{10 / \sqrt{100}} \geq \frac{C - 30}{10/\sqrt{100}} ~|~ \mu = 30\right) \\ + & = P\left(Z \geq \frac{C - 30}{1}\right) \\ +\end{align} +\]
+ +library(UsingR)
+data(father.son)
+t.test(father.son$sheight - father.son$fheight)
+##
+## One Sample t-test
+##
+## data: father.son$sheight - father.son$fheight
+## t = 11.79, df = 1077, p-value < 2.2e-16
+## alternative hypothesis: true mean is not equal to 0
+## 95 percent confidence interval:
+## 0.831 1.163
+## sample estimates:
+## mean of x
+## 0.997
+| Rejection region | +Type I error rate | +
|---|---|
| [0 : 8] | +1 | +
| [1 : 8] | +0.9961 | +
| [2 : 8] | +0.9648 | +
| [3 : 8] | +0.8555 | +
| [4 : 8] | +0.6367 | +
| [5 : 8] | +0.3633 | +
| [6 : 8] | +0.1445 | +
| [7 : 8] | +0.0352 | +
| [8 : 8] | +0.0039 | +
Brian Caffo, Jeffrey Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Idea: Suppose nothing is going on - how unusual is it to see the estimate we got?
- -Approach:
- -pt(0.8, 15, lower.tail = FALSE)
-[1] 0.2181
-choose(8, 7) * .5 ^ 8 + choose(8, 8) * .5 ^ 8
-[1] 0.03516
-pbinom(6, size = 8, prob = .5, lower.tail = FALSE)
-[1] 0.03516
-ppois(9, 5, lower.tail = FALSE)
-[1] 0.03183
-Brian Caffo, Jeffrey Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Idea: Suppose nothing is going on - how unusual is it to see the estimate we got?
+ +Approach:
+ +pt(0.8, 15, lower.tail = FALSE)
+## [1] 0.2181
+choose(8, 7) * 0.5^8 + choose(8, 8) * 0.5^8
+## [1] 0.03516
+pbinom(6, size = 8, prob = 0.5, lower.tail = FALSE)
+## [1] 0.03516
+ppois(9, 5, lower.tail = FALSE)
+## [1] 0.03183
+Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Assume that \(n\) is large and that we know \(\sigma\) -\[ -\begin{align} -1 -\beta & = -P\left(\frac{\bar X - 30}{\sigma /\sqrt{n}} > z_{1-\alpha} ~|~ \mu = \mu_a \right)\\ -& = P\left(\frac{\bar X - \mu_a + \mu_a - 30}{\sigma /\sqrt{n}} > z_{1-\alpha} ~|~ \mu = \mu_a \right)\\ \\ -& = P\left(\frac{\bar X - \mu_a}{\sigma /\sqrt{n}} > z_{1-\alpha} - \frac{\mu_a - 30}{\sigma /\sqrt{n}} ~|~ \mu = \mu_a \right)\\ \\ -& = P\left(Z > z_{1-\alpha} - \frac{\mu_a - 30}{\sigma /\sqrt{n}} ~|~ \mu = \mu_a \right)\\ \\ -\end{align} -\]
- -pnorm(-0.355, lower.tail = FALSE)
-[1] 0.6387
-sigma <- 10; mu_0 = 0; mu_a = 2; n <- 100; alpha = .05
-plot(c(-3, 6),c(0, dnorm(0)), type = "n", frame = false, xlab = "Z value", ylab = "")
-xvals <- seq(-3, 6, length = 1000)
-lines(xvals, dnorm(xvals), type = "l", lwd = 3)
-lines(xvals, dnorm(xvals, mean = sqrt(n) * (mu_a - mu_0) / sigma), lwd =3)
-abline(v = qnorm(1 - alpha))
-power.t.test does this very well
-
-power.t.test(n = 16, delta = 2 / 4, sd=1, type = "one.sample", alt = "one.sided")$power
-[1] 0.604
-power.t.test(n = 16, delta = 2, sd=4, type = "one.sample", alt = "one.sided")$power
-[1] 0.604
-power.t.test(n = 16, delta = 100, sd=200, type = "one.sample", alt = "one.sided")$power
-[1] 0.604
-power.t.test(power = .8, delta = 2 / 4, sd=1, type = "one.sample", alt = "one.sided")$n
-[1] 26.14
-power.t.test(power = .8, delta = 2, sd=4, type = "one.sample", alt = "one.sided")$n
-[1] 26.14
-power.t.test(power = .8, delta = 100, sd=200, type = "one.sample", alt = "one.sided")$n
-[1] 26.14
-Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
Assume that \(n\) is large and that we know \(\sigma\) +\[ +\begin{align} +1 -\beta & = +P\left(\frac{\bar X - 30}{\sigma /\sqrt{n}} > z_{1-\alpha} ~|~ \mu = \mu_a \right)\\ +& = P\left(\frac{\bar X - \mu_a + \mu_a - 30}{\sigma /\sqrt{n}} > z_{1-\alpha} ~|~ \mu = \mu_a \right)\\ \\ +& = P\left(\frac{\bar X - \mu_a}{\sigma /\sqrt{n}} > z_{1-\alpha} - \frac{\mu_a - 30}{\sigma /\sqrt{n}} ~|~ \mu = \mu_a \right)\\ \\ +& = P\left(Z > z_{1-\alpha} - \frac{\mu_a - 30}{\sigma /\sqrt{n}} ~|~ \mu = \mu_a \right)\\ \\ +\end{align} +\]
+ +pnorm(-0.355, lower.tail = FALSE)
+## [1] 0.6387
+sigma <- 10; mu_0 = 0; mu_a = 2; n <- 100; alpha = .05
+plot(c(-3, 6),c(0, dnorm(0)), type = "n", frame = FALSE, xlab = "Z value", ylab = "")
+xvals <- seq(-3, 6, length = 1000)
+lines(xvals, dnorm(xvals), type = "l", lwd = 3)
+lines(xvals, dnorm(xvals, mean = sqrt(n) * (mu_a - mu_0) / sigma), lwd =3)
+abline(v = qnorm(1 - alpha))
+
power.t.test does this very well
+
+power.t.test(n = 16, delta = 2/4, sd = 1, type = "one.sample", alt = "one.sided")$power
+## [1] 0.604
+power.t.test(n = 16, delta = 2, sd = 4, type = "one.sample", alt = "one.sided")$power
+## [1] 0.604
+power.t.test(n = 16, delta = 100, sd = 200, type = "one.sample", alt = "one.sided")$power
+## [1] 0.604
+power.t.test(power = 0.8, delta = 2/4, sd = 1, type = "one.sample", alt = "one.sided")$n
+## [1] 26.14
+power.t.test(power = 0.8, delta = 2, sd = 4, type = "one.sample", alt = "one.sided")$n
+## [1] 26.14
+power.t.test(power = 0.8, delta = 100, sd = 200, type = "one.sample", alt = "one.sided")$n
+## [1] 26.14
+
-
-
----
-
-## Why correct for multiple tests?
-
-
-
-
-[http://xkcd.com/882/](http://xkcd.com/882/)
-
----
-
-## Why correct for multiple tests?
-
-
-
-[http://xkcd.com/882/](http://xkcd.com/882/)
-
-
----
-
-## Types of errors
-
-Suppose you are testing a hypothesis that a parameter $\beta$ equals zero versus the alternative that it does not equal zero. These are the possible outcomes.
-
-
- | $\beta=0$ | $\beta\neq0$ | Hypotheses
---------------------|-------------|----------------|---------
-Claim $\beta=0$ | $U$ | $T$ | $m-R$
-Claim $\beta\neq 0$ | $V$ | $S$ | $R$
- Claims | $m_0$ | $m-m_0$ | $m$
-
-
-
-__Type I error or false positive ($V$)__ Say that the parameter does not equal zero when it does
-
-__Type II error or false negative ($T$)__ Say that the parameter equals zero when it doesn't
-
-
----
-
-## Error rates
-
-__False positive rate__ - The rate at which false results ($\beta = 0$) are called significant: $E\left[\frac{V}{m_0}\right]$*
-
-__Family wise error rate (FWER)__ - The probability of at least one false positive ${\rm Pr}(V \geq 1)$
-
-__False discovery rate (FDR)__ - The rate at which claims of significance are false $E\left[\frac{V}{R}\right]$
-
-* The false positive rate is closely related to the type I error rate [http://en.wikipedia.org/wiki/False_positive_rate](http://en.wikipedia.org/wiki/False_positive_rate)
-
----
-
-## Controlling the false positive rate
-
-If P-values are correctly calculated calling all $P < \alpha$ significant will control the false positive rate at level $\alpha$ on average.
-
-
-
-Controlling all error rates at $\alpha = 0.20$
-
----
-
-## Adjusted P-values
-
-* One approach is to adjust the threshold $\alpha$
-* A different approach is to calculate "adjusted p-values"
-* They _are not p-values_ anymore
-* But they can be used directly without adjusting $\alpha$
-
-__Example__:
-* Suppose P-values are $P_1,\ldots,P_m$
-* You could adjust them by taking $P_i^{fwer} = \max{m \times P_i,1}$ for each P-value.
-* Then if you call all $P_i^{fwer} < \alpha$ significant you will control the FWER.
-
----
-
-## Case study I: no true positives
-
-```{r createPvals,cache=TRUE}
-set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- y <- rnorm(20)
- x <- rnorm(20)
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-
-# Controls false positive rate
-sum(pValues < 0.05)
-```
-
----
-
-## Case study I: no true positives
-
-```{r, dependson="createPvals"}
-# Controls FWER
-sum(p.adjust(pValues,method="bonferroni") < 0.05)
-# Controls FDR
-sum(p.adjust(pValues,method="BH") < 0.05)
-```
-
-
----
-
-## Case study II: 50% true positives
-
-```{r createPvals2,cache=TRUE}
-set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- x <- rnorm(20)
- # First 500 beta=0, last 500 beta=2
- if(i <= 500){y <- rnorm(20)}else{ y <- rnorm(20,mean=2*x)}
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-trueStatus <- rep(c("zero","not zero"),each=500)
-table(pValues < 0.05, trueStatus)
-```
-
----
-
-
-## Case study II: 50% true positives
-
-```{r, dependson="createPvals2"}
-# Controls FWER
-table(p.adjust(pValues,method="bonferroni") < 0.05,trueStatus)
-# Controls FDR
-table(p.adjust(pValues,method="BH") < 0.05,trueStatus)
-```
-
-
----
-
-
-## Case study II: 50% true positives
-
-__P-values versus adjusted P-values__
-```{r, dependson="createPvals2",fig.height=4,fig.width=8}
-par(mfrow=c(1,2))
-plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
-plot(pValues,p.adjust(pValues,method="BH"),pch=19)
-```
-
-
----
-
-
-## Notes and resources
-
-__Notes__:
-* Multiple testing is an entire subfield
-* A basic Bonferroni/BH correction is usually enough
-* If there is strong dependence between tests there may be problems
- * Consider method="BY"
-
-__Further resources__:
-* [Multiple testing procedures with applications to genomics](http://www.amazon.com/Multiple-Procedures-Applications-Genomics-Statistics/dp/0387493166/ref=sr_1_2/102-3292576-129059?ie=UTF8&s=books&qid=1187394873&sr=1-2)
-* [Statistical significance for genome-wide studies](http://www.pnas.org/content/100/16/9440.full)
-* [Introduction to multiple testing](http://ies.ed.gov/ncee/pubs/20084018/app_b.asp)
-
+---
+title : Multiple testing
+subtitle : Statistical Inference
+author : Brian Caffo, Jeffrey Leek, Roger Peng
+job : Johns Hopkins Bloomberg School of Public Health
+logo : bloomberg_shield.png
+framework : io2012 # {io2012, html5slides, shower, dzslides, ...}
+highlighter : highlight.js # {highlight.js, prettify, highlight}
+hitheme : tomorrow #
+url:
+ lib: ../../librariesNew
+ assets: ../../assets
+widgets : [mathjax] # {mathjax, quiz, bootstrap}
+mode : selfcontained # {standalone, draft}
+---
+## Key ideas
+
+* Hypothesis testing/significance analysis is commonly overused
+* Correcting for multiple testing avoids false positives or discoveries
+* Two key components
+ * Error measure
+ * Correction
+
+
+---
+
+## Three eras of statistics
+
+__The age of Quetelet and his successors, in which huge census-level data sets were brought to bear on simple but important questions__: Are there more male than female births? Is the rate of insanity rising?
+
+The classical period of Pearson, Fisher, Neyman, Hotelling, and their successors, intellectual giants who __developed a theory of optimal inference capable of wringing every drop of information out of a scientific experiment__. The questions dealt with still tended to be simple Is treatment A better than treatment B?
+
+__The era of scientific mass production__, in which new technologies typified by the microarray allow a single team of scientists to produce data sets of a size Quetelet would envy. But now the flood of data is accompanied by a deluge of questions, perhaps thousands of estimates or hypothesis tests that the statistician is charged with answering together; not at all what the classical masters had in mind. Which variables matter among the thousands measured? How do you relate unrelated information?
+
+[http://www-stat.stanford.edu/~ckirby/brad/papers/2010LSIexcerpt.pdf](http://www-stat.stanford.edu/~ckirby/brad/papers/2010LSIexcerpt.pdf)
+
+---
+
+## Reasons for multiple testing
+
+
+
+
+---
+
+## Why correct for multiple tests?
+
+
+
+
+[http://xkcd.com/882/](http://xkcd.com/882/)
+
+---
+
+## Why correct for multiple tests?
+
+
+
+[http://xkcd.com/882/](http://xkcd.com/882/)
+
+
+---
+
+## Types of errors
+
+Suppose you are testing a hypothesis that a parameter $\beta$ equals zero versus the alternative that it does not equal zero. These are the possible outcomes.
+
+
+ | $\beta=0$ | $\beta\neq0$ | Hypotheses
+--------------------|-------------|----------------|---------
+Claim $\beta=0$ | $U$ | $T$ | $m-R$
+Claim $\beta\neq 0$ | $V$ | $S$ | $R$
+ Claims | $m_0$ | $m-m_0$ | $m$
+
+
+
+__Type I error or false positive ($V$)__ Say that the parameter does not equal zero when it does
+
+__Type II error or false negative ($T$)__ Say that the parameter equals zero when it doesn't
+
+
+---
+
+## Error rates
+
+__False positive rate__ - The rate at which false results ($\beta = 0$) are called significant: $E\left[\frac{V}{m_0}\right]$*
+
+__Family wise error rate (FWER)__ - The probability of at least one false positive ${\rm Pr}(V \geq 1)$
+
+__False discovery rate (FDR)__ - The rate at which claims of significance are false $E\left[\frac{V}{R}\right]$
+
+* The false positive rate is closely related to the type I error rate [http://en.wikipedia.org/wiki/False_positive_rate](http://en.wikipedia.org/wiki/False_positive_rate)
+
+---
+
+## Controlling the false positive rate
+
+If P-values are correctly calculated calling all $P < \alpha$ significant will control the false positive rate at level $\alpha$ on average.
+
+
+
+Controlling all error rates at $\alpha = 0.20$
+
+---
+
+## Adjusted P-values
+
+* One approach is to adjust the threshold $\alpha$
+* A different approach is to calculate "adjusted p-values"
+* They _are not p-values_ anymore
+* But they can be used directly without adjusting $\alpha$
+
+__Example__:
+* Suppose P-values are $P_1,\ldots,P_m$
+* You could adjust them by taking $P_i^{fwer} = \max{m \times P_i,1}$ for each P-value.
+* Then if you call all $P_i^{fwer} < \alpha$ significant you will control the FWER.
+
+---
+
+## Case study I: no true positives
+
+```{r createPvals,cache=TRUE}
+set.seed(1010093)
+pValues <- rep(NA,1000)
+for(i in 1:1000){
+ y <- rnorm(20)
+ x <- rnorm(20)
+ pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
+}
+
+# Controls false positive rate
+sum(pValues < 0.05)
+```
+
+---
+
+## Case study I: no true positives
+
+```{r, dependson="createPvals"}
+# Controls FWER
+sum(p.adjust(pValues,method="bonferroni") < 0.05)
+# Controls FDR
+sum(p.adjust(pValues,method="BH") < 0.05)
+```
+
+
+---
+
+## Case study II: 50% true positives
+
+```{r createPvals2,cache=TRUE}
+set.seed(1010093)
+pValues <- rep(NA,1000)
+for(i in 1:1000){
+ x <- rnorm(20)
+ # First 500 beta=0, last 500 beta=2
+ if(i <= 500){y <- rnorm(20)}else{ y <- rnorm(20,mean=2*x)}
+ pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
+}
+trueStatus <- rep(c("zero","not zero"),each=500)
+table(pValues < 0.05, trueStatus)
+```
+
+---
+
+
+## Case study II: 50% true positives
+
+```{r, dependson="createPvals2"}
+# Controls FWER
+table(p.adjust(pValues,method="bonferroni") < 0.05,trueStatus)
+# Controls FDR
+table(p.adjust(pValues,method="BH") < 0.05,trueStatus)
+```
+
+
+---
+
+
+## Case study II: 50% true positives
+
+__P-values versus adjusted P-values__
+```{r, dependson="createPvals2",fig.height=4,fig.width=8}
+par(mfrow=c(1,2))
+plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
+plot(pValues,p.adjust(pValues,method="BH"),pch=19)
+```
+
+
+---
+
+
+## Notes and resources
+
+__Notes__:
+* Multiple testing is an entire subfield
+* A basic Bonferroni/BH correction is usually enough
+* If there is strong dependence between tests there may be problems
+ * Consider method="BY"
+
+__Further resources__:
+* [Multiple testing procedures with applications to genomics](http://www.amazon.com/Multiple-Procedures-Applications-Genomics-Statistics/dp/0387493166/ref=sr_1_2/102-3292576-129059?ie=UTF8&s=books&qid=1187394873&sr=1-2)
+* [Statistical significance for genome-wide studies](http://www.pnas.org/content/100/16/9440.full)
+* [Introduction to multiple testing](http://ies.ed.gov/ncee/pubs/20084018/app_b.asp)
+
diff --git a/06_StatisticalInference/03_05_MultipleTesting/index.html b/06_StatisticalInference/03_05_MultipleTesting/index.html
index bb3a271db..dfc498eeb 100644
--- a/06_StatisticalInference/03_05_MultipleTesting/index.html
+++ b/06_StatisticalInference/03_05_MultipleTesting/index.html
@@ -1,581 +1,585 @@
-
-
-
- Brian Caffo, Jeffrey Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
The age of Quetelet and his successors, in which huge census-level data sets were brought to bear on simple but important questions: Are there more male than female births? Is the rate of insanity rising?
- -The classical period of Pearson, Fisher, Neyman, Hotelling, and their successors, intellectual giants who developed a theory of optimal inference capable of wringing every drop of information out of a scientific experiment. The questions dealt with still tended to be simple Is treatment A better than treatment B?
- -The era of scientific mass production, in which new technologies typified by the microarray allow a single team of scientists to produce data sets of a size Quetelet would envy. But now the flood of data is accompanied by a deluge of questions, perhaps thousands of estimates or hypothesis tests that the statistician is charged with answering together; not at all what the classical masters had in mind. Which variables matter among the thousands measured? How do you relate unrelated information?
- -http://www-stat.stanford.edu/~ckirby/brad/papers/2010LSIexcerpt.pdf
- -
Suppose you are testing a hypothesis that a parameter \(\beta\) equals zero versus the alternative that it does not equal zero. These are the possible outcomes. -
- -| - | \(\beta=0\) | -\(\beta\neq0\) | -Hypotheses | -
|---|---|---|---|
| Claim \(\beta=0\) | -\(U\) | -\(T\) | -\(m-R\) | -
| Claim \(\beta\neq 0\) | -\(V\) | -\(S\) | -\(R\) | -
| Claims | -\(m_0\) | -\(m-m_0\) | -\(m\) | -
Type I error or false positive (\(V\)) Say that the parameter does not equal zero when it does
- -Type II error or false negative (\(T\)) Say that the parameter equals zero when it doesn't
- -False positive rate - The rate at which false results (\(\beta = 0\)) are called significant: \(E\left[\frac{V}{m_0}\right]\)*
- -Family wise error rate (FWER) - The probability of at least one false positive \({\rm Pr}(V \geq 1)\)
- -False discovery rate (FDR) - The rate at which claims of significance are false \(E\left[\frac{V}{R}\right]\)
- -If P-values are correctly calculated calling all \(P < \alpha\) significant will control the false positive rate at level \(\alpha\) on average.
- -Suppose that you call all \(P < 0.05\) significant.
- -The expected number of false positives is: \(10,000 \times 0.05 = 500\) false positives.
- -How do we avoid so many false positives?
- -The Bonferroni correction is the oldest multiple testing correction.
- -Basic idea:
- -Pros: Easy to calculate, conservative -Cons: May be very conservative
- -This is the most popular correction when performing lots of tests say in genomics, imaging, astronomy, or other signal-processing disciplines.
- -Basic idea:
- -Pros: Still pretty easy to calculate, less conservative (maybe much less)
- -Cons: Allows for more false positives, may behave strangely under dependence
- -
Controlling all error rates at \(\alpha = 0.20\)
- -Example:
- -set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- y <- rnorm(20)
- x <- rnorm(20)
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-
-# Controls false positive rate
-sum(pValues < 0.05)
-[1] 51
-# Controls FWER
-sum(p.adjust(pValues,method="bonferroni") < 0.05)
-[1] 0
-# Controls FDR
-sum(p.adjust(pValues,method="BH") < 0.05)
-[1] 0
-set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- x <- rnorm(20)
- # First 500 beta=0, last 500 beta=2
- if(i <= 500){y <- rnorm(20)}else{ y <- rnorm(20,mean=2*x)}
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-trueStatus <- rep(c("zero","not zero"),each=500)
-table(pValues < 0.05, trueStatus)
- trueStatus
- not zero zero
- FALSE 0 476
- TRUE 500 24
-# Controls FWER
-table(p.adjust(pValues,method="bonferroni") < 0.05,trueStatus)
- trueStatus
- not zero zero
- FALSE 23 500
- TRUE 477 0
-# Controls FDR
-table(p.adjust(pValues,method="BH") < 0.05,trueStatus)
- trueStatus
- not zero zero
- FALSE 0 487
- TRUE 500 13
-P-values versus adjusted P-values
- -par(mfrow=c(1,2))
-plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
-plot(pValues,p.adjust(pValues,method="BH"),pch=19)
-Notes:
- -Further resources:
- - - -Brian Caffo, Jeffrey Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
The age of Quetelet and his successors, in which huge census-level data sets were brought to bear on simple but important questions: Are there more male than female births? Is the rate of insanity rising?
+ +The classical period of Pearson, Fisher, Neyman, Hotelling, and their successors, intellectual giants who developed a theory of optimal inference capable of wringing every drop of information out of a scientific experiment. The questions dealt with still tended to be simple Is treatment A better than treatment B?
+ +The era of scientific mass production, in which new technologies typified by the microarray allow a single team of scientists to produce data sets of a size Quetelet would envy. But now the flood of data is accompanied by a deluge of questions, perhaps thousands of estimates or hypothesis tests that the statistician is charged with answering together; not at all what the classical masters had in mind. Which variables matter among the thousands measured? How do you relate unrelated information?
+ +http://www-stat.stanford.edu/~ckirby/brad/papers/2010LSIexcerpt.pdf
+ +
Suppose you are testing a hypothesis that a parameter \(\beta\) equals zero versus the alternative that it does not equal zero. These are the possible outcomes. +
+ +| + | \(\beta=0\) | +\(\beta\neq0\) | +Hypotheses | +
|---|---|---|---|
| Claim \(\beta=0\) | +\(U\) | +\(T\) | +\(m-R\) | +
| Claim \(\beta\neq 0\) | +\(V\) | +\(S\) | +\(R\) | +
| Claims | +\(m_0\) | +\(m-m_0\) | +\(m\) | +
Type I error or false positive (\(V\)) Say that the parameter does not equal zero when it does
+ +Type II error or false negative (\(T\)) Say that the parameter equals zero when it doesn't
+ +False positive rate - The rate at which false results (\(\beta = 0\)) are called significant: \(E\left[\frac{V}{m_0}\right]\)*
+ +Family wise error rate (FWER) - The probability of at least one false positive \({\rm Pr}(V \geq 1)\)
+ +False discovery rate (FDR) - The rate at which claims of significance are false \(E\left[\frac{V}{R}\right]\)
+ +If P-values are correctly calculated calling all \(P < \alpha\) significant will control the false positive rate at level \(\alpha\) on average.
+ +Suppose that you call all \(P < 0.05\) significant.
+ +The expected number of false positives is: \(10,000 \times 0.05 = 500\) false positives.
+ +How do we avoid so many false positives?
+ +The Bonferroni correction is the oldest multiple testing correction.
+ +Basic idea:
+ +Pros: Easy to calculate, conservative +Cons: May be very conservative
+ +This is the most popular correction when performing lots of tests say in genomics, imaging, astronomy, or other signal-processing disciplines.
+ +Basic idea:
+ +Pros: Still pretty easy to calculate, less conservative (maybe much less)
+ +Cons: Allows for more false positives, may behave strangely under dependence
+ +
Controlling all error rates at \(\alpha = 0.20\)
+ +Example:
+ +set.seed(1010093)
+pValues <- rep(NA, 1000)
+for (i in 1:1000) {
+ y <- rnorm(20)
+ x <- rnorm(20)
+ pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
+}
+
+# Controls false positive rate
+sum(pValues < 0.05)
+## [1] 51
+# Controls FWER
+sum(p.adjust(pValues, method = "bonferroni") < 0.05)
+## [1] 0
+# Controls FDR
+sum(p.adjust(pValues, method = "BH") < 0.05)
+## [1] 0
+set.seed(1010093)
+pValues <- rep(NA, 1000)
+for (i in 1:1000) {
+ x <- rnorm(20)
+ # First 500 beta=0, last 500 beta=2
+ if (i <= 500) {
+ y <- rnorm(20)
+ } else {
+ y <- rnorm(20, mean = 2 * x)
+ }
+ pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
+}
+trueStatus <- rep(c("zero", "not zero"), each = 500)
+table(pValues < 0.05, trueStatus)
+## trueStatus
+## not zero zero
+## FALSE 0 476
+## TRUE 500 24
+# Controls FWER
+table(p.adjust(pValues, method = "bonferroni") < 0.05, trueStatus)
+## trueStatus
+## not zero zero
+## FALSE 23 500
+## TRUE 477 0
+# Controls FDR
+table(p.adjust(pValues, method = "BH") < 0.05, trueStatus)
+## trueStatus
+## not zero zero
+## FALSE 0 487
+## TRUE 500 13
+P-values versus adjusted P-values
+ +par(mfrow = c(1, 2))
+plot(pValues, p.adjust(pValues, method = "bonferroni"), pch = 19)
+plot(pValues, p.adjust(pValues, method = "BH"), pch = 19)
+
Notes:
+ +Further resources:
+ + + +
-
-
----
-
-## Why correct for multiple tests?
-
-
-
-
-[http://xkcd.com/882/](http://xkcd.com/882/)
-
----
-
-## Why correct for multiple tests?
-
-
-
-[http://xkcd.com/882/](http://xkcd.com/882/)
-
-
----
-
-## Types of errors
-
-Suppose you are testing a hypothesis that a parameter $\beta$ equals zero versus the alternative that it does not equal zero. These are the possible outcomes.
-
-
- | $\beta=0$ | $\beta\neq0$ | Hypotheses
---------------------|-------------|----------------|---------
-Claim $\beta=0$ | $U$ | $T$ | $m-R$
-Claim $\beta\neq 0$ | $V$ | $S$ | $R$
- Claims | $m_0$ | $m-m_0$ | $m$
-
-
-
-__Type I error or false positive ($V$)__ Say that the parameter does not equal zero when it does
-
-__Type II error or false negative ($T$)__ Say that the parameter equals zero when it doesn't
-
-
----
-
-## Error rates
-
-__False positive rate__ - The rate at which false results ($\beta = 0$) are called significant: $E\left[\frac{V}{m_0}\right]$*
-
-__Family wise error rate (FWER)__ - The probability of at least one false positive ${\rm Pr}(V \geq 1)$
-
-__False discovery rate (FDR)__ - The rate at which claims of significance are false $E\left[\frac{V}{R}\right]$
-
-* The false positive rate is closely related to the type I error rate [http://en.wikipedia.org/wiki/False_positive_rate](http://en.wikipedia.org/wiki/False_positive_rate)
-
----
-
-## Controlling the false positive rate
-
-If P-values are correctly calculated calling all $P < \alpha$ significant will control the false positive rate at level $\alpha$ on average.
-
-
-
-Controlling all error rates at $\alpha = 0.20$
-
----
-
-## Adjusted P-values
-
-* One approach is to adjust the threshold $\alpha$
-* A different approach is to calculate "adjusted p-values"
-* They _are not p-values_ anymore
-* But they can be used directly without adjusting $\alpha$
-
-__Example__:
-* Suppose P-values are $P_1,\ldots,P_m$
-* You could adjust them by taking $P_i^{fwer} = \max{m \times P_i,1}$ for each P-value.
-* Then if you call all $P_i^{fwer} < \alpha$ significant you will control the FWER.
-
----
-
-## Case study I: no true positives
-
-
-```r
-set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- y <- rnorm(20)
- x <- rnorm(20)
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-
-# Controls false positive rate
-sum(pValues < 0.05)
-```
-
-```
-[1] 51
-```
-
-
----
-
-## Case study I: no true positives
-
-
-```r
-# Controls FWER
-sum(p.adjust(pValues,method="bonferroni") < 0.05)
-```
-
-```
-[1] 0
-```
-
-```r
-# Controls FDR
-sum(p.adjust(pValues,method="BH") < 0.05)
-```
-
-```
-[1] 0
-```
-
-
-
----
-
-## Case study II: 50% true positives
-
-
-```r
-set.seed(1010093)
-pValues <- rep(NA,1000)
-for(i in 1:1000){
- x <- rnorm(20)
- # First 500 beta=0, last 500 beta=2
- if(i <= 500){y <- rnorm(20)}else{ y <- rnorm(20,mean=2*x)}
- pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
-}
-trueStatus <- rep(c("zero","not zero"),each=500)
-table(pValues < 0.05, trueStatus)
-```
-
-```
- trueStatus
- not zero zero
- FALSE 0 476
- TRUE 500 24
-```
-
-
----
-
-
-## Case study II: 50% true positives
-
-
-```r
-# Controls FWER
-table(p.adjust(pValues,method="bonferroni") < 0.05,trueStatus)
-```
-
-```
- trueStatus
- not zero zero
- FALSE 23 500
- TRUE 477 0
-```
-
-```r
-# Controls FDR
-table(p.adjust(pValues,method="BH") < 0.05,trueStatus)
-```
-
-```
- trueStatus
- not zero zero
- FALSE 0 487
- TRUE 500 13
-```
-
-
-
----
-
-
-## Case study II: 50% true positives
-
-__P-values versus adjusted P-values__
-
-```r
-par(mfrow=c(1,2))
-plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
-plot(pValues,p.adjust(pValues,method="BH"),pch=19)
-```
-
-
+
+
+---
+
+## Why correct for multiple tests?
+
+
+
+
+[http://xkcd.com/882/](http://xkcd.com/882/)
+
+---
+
+## Why correct for multiple tests?
+
+
+
+[http://xkcd.com/882/](http://xkcd.com/882/)
+
+
+---
+
+## Types of errors
+
+Suppose you are testing a hypothesis that a parameter $\beta$ equals zero versus the alternative that it does not equal zero. These are the possible outcomes.
+
+
+ | $\beta=0$ | $\beta\neq0$ | Hypotheses
+--------------------|-------------|----------------|---------
+Claim $\beta=0$ | $U$ | $T$ | $m-R$
+Claim $\beta\neq 0$ | $V$ | $S$ | $R$
+ Claims | $m_0$ | $m-m_0$ | $m$
+
+
+
+__Type I error or false positive ($V$)__ Say that the parameter does not equal zero when it does
+
+__Type II error or false negative ($T$)__ Say that the parameter equals zero when it doesn't
+
+
+---
+
+## Error rates
+
+__False positive rate__ - The rate at which false results ($\beta = 0$) are called significant: $E\left[\frac{V}{m_0}\right]$*
+
+__Family wise error rate (FWER)__ - The probability of at least one false positive ${\rm Pr}(V \geq 1)$
+
+__False discovery rate (FDR)__ - The rate at which claims of significance are false $E\left[\frac{V}{R}\right]$
+
+* The false positive rate is closely related to the type I error rate [http://en.wikipedia.org/wiki/False_positive_rate](http://en.wikipedia.org/wiki/False_positive_rate)
+
+---
+
+## Controlling the false positive rate
+
+If P-values are correctly calculated calling all $P < \alpha$ significant will control the false positive rate at level $\alpha$ on average.
+
+
+
+Controlling all error rates at $\alpha = 0.20$
+
+---
+
+## Adjusted P-values
+
+* One approach is to adjust the threshold $\alpha$
+* A different approach is to calculate "adjusted p-values"
+* They _are not p-values_ anymore
+* But they can be used directly without adjusting $\alpha$
+
+__Example__:
+* Suppose P-values are $P_1,\ldots,P_m$
+* You could adjust them by taking $P_i^{fwer} = \max{m \times P_i,1}$ for each P-value.
+* Then if you call all $P_i^{fwer} < \alpha$ significant you will control the FWER.
+
+---
+
+## Case study I: no true positives
+
+
+```r
+set.seed(1010093)
+pValues <- rep(NA, 1000)
+for (i in 1:1000) {
+ y <- rnorm(20)
+ x <- rnorm(20)
+ pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
+}
+
+# Controls false positive rate
+sum(pValues < 0.05)
+```
+
+```
+## [1] 51
+```
+
+
+---
+
+## Case study I: no true positives
+
+
+```r
+# Controls FWER
+sum(p.adjust(pValues, method = "bonferroni") < 0.05)
+```
+
+```
+## [1] 0
+```
+
+```r
+# Controls FDR
+sum(p.adjust(pValues, method = "BH") < 0.05)
+```
+
+```
+## [1] 0
+```
+
+
+
+---
+
+## Case study II: 50% true positives
+
+
+```r
+set.seed(1010093)
+pValues <- rep(NA, 1000)
+for (i in 1:1000) {
+ x <- rnorm(20)
+ # First 500 beta=0, last 500 beta=2
+ if (i <= 500) {
+ y <- rnorm(20)
+ } else {
+ y <- rnorm(20, mean = 2 * x)
+ }
+ pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
+}
+trueStatus <- rep(c("zero", "not zero"), each = 500)
+table(pValues < 0.05, trueStatus)
+```
+
+```
+## trueStatus
+## not zero zero
+## FALSE 0 476
+## TRUE 500 24
+```
+
+
+---
+
+
+## Case study II: 50% true positives
+
+
+```r
+# Controls FWER
+table(p.adjust(pValues, method = "bonferroni") < 0.05, trueStatus)
+```
+
+```
+## trueStatus
+## not zero zero
+## FALSE 23 500
+## TRUE 477 0
+```
+
+```r
+# Controls FDR
+table(p.adjust(pValues, method = "BH") < 0.05, trueStatus)
+```
+
+```
+## trueStatus
+## not zero zero
+## FALSE 0 487
+## TRUE 500 13
+```
+
+
+
+---
+
+
+## Case study II: 50% true positives
+
+__P-values versus adjusted P-values__
+
+```r
+par(mfrow = c(1, 2))
+plot(pValues, p.adjust(pValues, method = "bonferroni"), pch = 19)
+plot(pValues, p.adjust(pValues, method = "BH"), pch = 19)
+```
+
+
+
+
+
+---
+
+
+## Notes and resources
+
+__Notes__:
+* Multiple testing is an entire subfield
+* A basic Bonferroni/BH correction is usually enough
+* If there is strong dependence between tests there may be problems
+ * Consider method="BY"
+
+__Further resources__:
+* [Multiple testing procedures with applications to genomics](http://www.amazon.com/Multiple-Procedures-Applications-Genomics-Statistics/dp/0387493166/ref=sr_1_2/102-3292576-129059?ie=UTF8&s=books&qid=1187394873&sr=1-2)
+* [Statistical significance for genome-wide studies](http://www.pnas.org/content/100/16/9440.full)
+* [Introduction to multiple testing](http://ies.ed.gov/ncee/pubs/20084018/app_b.asp)
+
diff --git a/06_StatisticalInference/03_05_MultipleTesting/index.pdf b/06_StatisticalInference/03_05_MultipleTesting/index.pdf
index 190c24c34..88d17ad14 100644
Binary files a/06_StatisticalInference/03_05_MultipleTesting/index.pdf and b/06_StatisticalInference/03_05_MultipleTesting/index.pdf differ
diff --git a/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-4.png b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-4.png
new file mode 100644
index 000000000..d88aaaddd
Binary files /dev/null and b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-4.png differ
diff --git a/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-5.png b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-5.png
new file mode 100644
index 000000000..57f92bdae
Binary files /dev/null and b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-5.png differ
diff --git a/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-7.png b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-7.png
new file mode 100644
index 000000000..8542ea3a4
Binary files /dev/null and b/06_StatisticalInference/03_06_resampledInference/assets/fig/unnamed-chunk-7.png differ
diff --git a/06_StatisticalInference/03_06_resampledInference/index.Rmd b/06_StatisticalInference/03_06_resampledInference/index.Rmd
index 775099c76..3f25e8c38 100644
--- a/06_StatisticalInference/03_06_resampledInference/index.Rmd
+++ b/06_StatisticalInference/03_06_resampledInference/index.Rmd
@@ -1,259 +1,243 @@
----
-title : Resampled inference
-subtitle : Statistical Inference
-author : Brian Caffo, Jeff Leek, Roger Peng
-job : Johns Hopkins Bloomberg School of Public Health
-logo : bloomberg_shield.png
-framework : io2012 # {io2012, html5slides, shower, dzslides, ...}
-highlighter : highlight.js # {highlight.js, prettify, highlight}
-hitheme : tomorrow #
-url:
- lib: ../../librariesNew
- assets: ../../assets
-widgets : [mathjax] # {mathjax, quiz, bootstrap}
-mode : selfcontained # {standalone, draft}
-
----
-```{r setup, cache = F, echo = F, message = F, warning = F, tidy = F, results='hide'}
-# make this an external chunk that can be included in any file
-options(width = 100)
-opts_chunk$set(message = F, error = F, warning = F, comment = NA, fig.align = 'center', dpi = 100, tidy = F, cache.path = '.cache/', fig.path = 'fig/')
-
-options(xtable.type = 'html')
-knit_hooks$set(inline = function(x) {
- if(is.numeric(x)) {
- round(x, getOption('digits'))
- } else {
- paste(as.character(x), collapse = ', ')
- }
-})
-knit_hooks$set(plot = knitr:::hook_plot_html)
-runif(1)
-```
-
-## The jackknife
-
-- The jackknife is a tool for estimating standard errors and the bias of estimators
-- As its name suggests, the jackknife is a small, handy tool; in contrast to the bootstrap, which is then the moral equivalent of a giant workshop full of tools
-- Both the jackknife and the bootstrap involve *resampling* data; that is, repeatedly creating new data sets from the original data
-
----
-
-## The jackknife
-
-- The jackknife deletes each observation and calculates an estimate based on the remaining $n-1$ of them
-- It uses this collection of estimates to do things like estimate the bias and the standard error
-- Note that estimating the bias and having a standard error are not needed for things like sample means, which we know are unbiased estimates of population means and what their standard errors are
-
----
-
-## The jackknife
-
-- We'll consider the jackknife for univariate data
-- Let $X_1,\ldots,X_n$ be a collection of data used to estimate a parameter $\theta$
-- Let $\hat \theta$ be the estimate based on the full data set
-- Let $\hat \theta_{i}$ be the estimate of $\theta$ obtained by *deleting observation $i$*
-- Let $\bar \theta = \frac{1}{n}\sum_{i=1}^n \hat \theta_{i}$
-
----
-
-## Continued
-
-- Then, the jackknife estimate of the bias is
- $$
- (n - 1) \left(\bar \theta - \hat \theta\right)
- $$
- (how far the average delete-one estimate is from the actual estimate)
-- The jackknife estimate of the standard error is
- $$
- \left[\frac{n-1}{n}\sum_{i=1}^n (\hat \theta_i - \bar\theta )^2\right]^{1/2}
- $$
-(the deviance of the delete-one estimates from the average delete-one estimate)
-
----
-
-## Example
-### We want to estimate the bias and standard error of the median
-
-```{r, results='hide'}
-library(UsingR)
-data(father.son)
-x <- father.son$sheight
-n <- length(x)
-theta <- median(x)
-jk <- sapply(1 : n,
- function(i) median(x[-i])
- )
-thetaBar <- mean(jk)
-biasEst <- (n - 1) * (thetaBar - theta)
-seEst <- sqrt((n - 1) * mean((jk - thetaBar)^2))
-```
-
----
-
-## Example
-
-```{r}
-c(biasEst, seEst)
-library(bootstrap)
-temp <- jackknife(x, median)
-c(temp$jack.bias, temp$jack.se)
-```
-
----
-
-## Example
-
-- Both methods (of course) yield an estimated bias of `r temp$jack.bias` and a se of `r temp$jack.se`
-- Odd little fact: the jackknife estimate of the bias for the median is always $0$ when the number of observations is even
-- It has been shown that the jackknife is a linear approximation to the bootstrap
-- Generally do not use the jackknife for sample quantiles like the median; as it has been shown to have some poor properties
-
----
-
-## Pseudo observations
-
-- Another interesting way to think about the jackknife uses pseudo observations
-- Let
-$$
- \mbox{Pseudo Obs} = n \hat \theta - (n - 1) \hat \theta_{i}
-$$
-- Think of these as ``whatever observation $i$ contributes to the estimate of $\theta$''
-- Note when $\hat \theta$ is the sample mean, the pseudo observations are the data themselves
-- Then the sample standard error of these observations is the previous jackknife estimated standard error.
-- The mean of these observations is a bias-corrected estimate of $\theta$
-
----
-
-## The bootstrap
-
-- The bootstrap is a tremendously useful tool for constructing confidence intervals and calculating standard errors for difficult statistics
-- For example, how would one derive a confidence interval for the median?
-- The bootstrap procedure follows from the so called bootstrap principle
-
----
-
-## The bootstrap principle
-
-- Suppose that I have a statistic that estimates some population parameter, but I don't know its sampling distribution
-- The bootstrap principle suggests using the distribution defined by the data to approximate its sampling distribution
-
----
-
-## The bootstrap in practice
-
-- In practice, the bootstrap principle is always carried out using simulation
-- We will cover only a few aspects of bootstrap resampling
-- The general procedure follows by first simulating complete data sets from the observed data with replacement
-
- - This is approximately drawing from the sampling distribution of that statistic, at least as far as the data is able to approximate the true population distribution
-
-- Calculate the statistic for each simulated data set
-- Use the simulated statistics to either define a confidence interval or take the standard deviation to calculate a standard error
-
----
-## Nonparametric bootstrap algorithm example
-
-- Bootstrap procedure for calculating confidence interval for the median from a data set of $n$ observations
-
- i. Sample $n$ observations **with replacement** from the observed data resulting in one simulated complete data set
-
- ii. Take the median of the simulated data set
-
- iii. Repeat these two steps $B$ times, resulting in $B$ simulated medians
-
- iv. These medians are approximately drawn from the sampling distribution of the median of $n$ observations; therefore we can
-
- - Draw a histogram of them
- - Calculate their standard deviation to estimate the standard error of the median
- - Take the $2.5^{th}$ and $97.5^{th}$ percentiles as a confidence interval for the median
-
----
-
-## Example code
-
-```{r}
-B <- 1000
-resamples <- matrix(sample(x,
- n * B,
- replace = TRUE),
- B, n)
-medians <- apply(resamples, 1, median)
-sd(medians)
-quantile(medians, c(.025, .975))
-```
-
----
-## Histogram of bootstrap resamples
-
-```{r, fig.height=5, fig.width=5}
-hist(medians)
-```
-
----
-
-## Notes on the bootstrap
-
-- The bootstrap is non-parametric
-- Better percentile bootstrap confidence intervals correct for bias
-- There are lots of variations on bootstrap procedures; the book "An Introduction to the Bootstrap"" by Efron and Tibshirani is a great place to start for both bootstrap and jackknife information
-
-
----
-## Group comparisons
-- Consider comparing two independent groups.
-- Example, comparing sprays B and C
-
-```{r, fig.height=4, fig.width=4}
-data(InsectSprays)
-boxplot(count ~ spray, data = InsectSprays)
-```
-
----
-## Permutation tests
-- Consider the null hypothesis that the distribution of the observations from each group is the same
-- Then, the group labels are irrelevant
-- We then discard the group levels and permute the combined data
-- Split the permuted data into two groups with $n_A$ and $n_B$
- observations (say by always treating the first $n_A$ observations as
- the first group)
-- Evaluate the probability of getting a statistic as large or
- large than the one observed
-- An example statistic would be the difference in the averages between the two groups;
- one could also use a t-statistic
-
----
-## Variations on permutation testing
-Data type | Statistic | Test name
----|---|---|
-Ranks | rank sum | rank sum test
-Binary | hypergeometric prob | Fisher's exact test
-Raw data | | ordinary permutation test
-
-- Also, so-called *randomization tests* are exactly permutation tests, with a different motivation.
-- For matched data, one can randomize the signs
- - For ranks, this results in the signed rank test
-- Permutation strategies work for regression as well
- - Permuting a regressor of interest
-- Permutation tests work very well in multivariate settings
-
----
-## Permutation test for pesticide data
-```{r}
-subdata <- InsectSprays[InsectSprays$spray %in% c("B", "C"),]
-y <- subdata$count
-group <- as.character(subdata$spray)
-testStat <- function(w, g) mean(w[g == "B"]) - mean(w[g == "C"])
-observedStat <- testStat(y, group)
-permutations <- sapply(1 : 10000, function(i) testStat(y, sample(group)))
-observedStat
-mean(permutations > observedStat)
-```
-
----
-## Histogram of permutations
-```{r, echo= FALSE, fig.width=5, fig.height=5}
-hist(permutations)
-```
-
-
+---
+title : Resampled inference
+subtitle : Statistical Inference
+author : Brian Caffo, Jeff Leek, Roger Peng
+job : Johns Hopkins Bloomberg School of Public Health
+logo : bloomberg_shield.png
+framework : io2012 # {io2012, html5slides, shower, dzslides, ...}
+highlighter : highlight.js # {highlight.js, prettify, highlight}
+hitheme : tomorrow #
+url:
+ lib: ../../librariesNew
+ assets: ../../assets
+widgets : [mathjax] # {mathjax, quiz, bootstrap}
+mode : selfcontained # {standalone, draft}
+
+---
+
+## The jackknife
+
+- The jackknife is a tool for estimating standard errors and the bias of estimators
+- As its name suggests, the jackknife is a small, handy tool; in contrast to the bootstrap, which is then the moral equivalent of a giant workshop full of tools
+- Both the jackknife and the bootstrap involve *resampling* data; that is, repeatedly creating new data sets from the original data
+
+---
+
+## The jackknife
+
+- The jackknife deletes each observation and calculates an estimate based on the remaining $n-1$ of them
+- It uses this collection of estimates to do things like estimate the bias and the standard error
+- Note that estimating the bias and having a standard error are not needed for things like sample means, which we know are unbiased estimates of population means and what their standard errors are
+
+---
+
+## The jackknife
+
+- We'll consider the jackknife for univariate data
+- Let $X_1,\ldots,X_n$ be a collection of data used to estimate a parameter $\theta$
+- Let $\hat \theta$ be the estimate based on the full data set
+- Let $\hat \theta_{i}$ be the estimate of $\theta$ obtained by *deleting observation $i$*
+- Let $\bar \theta = \frac{1}{n}\sum_{i=1}^n \hat \theta_{i}$
+
+---
+
+## Continued
+
+- Then, the jackknife estimate of the bias is
+ $$
+ (n - 1) \left(\bar \theta - \hat \theta\right)
+ $$
+ (how far the average delete-one estimate is from the actual estimate)
+- The jackknife estimate of the standard error is
+ $$
+ \left[\frac{n-1}{n}\sum_{i=1}^n (\hat \theta_i - \bar\theta )^2\right]^{1/2}
+ $$
+(the deviance of the delete-one estimates from the average delete-one estimate)
+
+---
+
+## Example
+### We want to estimate the bias and standard error of the median
+
+```{r, results='hide'}
+library(UsingR)
+data(father.son)
+x <- father.son$sheight
+n <- length(x)
+theta <- median(x)
+jk <- sapply(1 : n,
+ function(i) median(x[-i])
+ )
+thetaBar <- mean(jk)
+biasEst <- (n - 1) * (thetaBar - theta)
+seEst <- sqrt((n - 1) * mean((jk - thetaBar)^2))
+```
+
+---
+
+## Example test
+
+```{r}
+c(biasEst, seEst)
+library(bootstrap)
+temp <- jackknife(x, median)
+c(temp$jack.bias, temp$jack.se)
+```
+
+---
+
+## Example
+
+- Both methods (of course) yield an estimated bias of `r temp$jack.bias` and a se of `r temp$jack.se`
+- Odd little fact: the jackknife estimate of the bias for the median is always $0$ when the number of observations is even
+- It has been shown that the jackknife is a linear approximation to the bootstrap
+- Generally do not use the jackknife for sample quantiles like the median; as it has been shown to have some poor properties
+
+---
+
+## Pseudo observations
+
+- Another interesting way to think about the jackknife uses pseudo observations
+- Let
+$$
+ \mbox{Pseudo Obs} = n \hat \theta - (n - 1) \hat \theta_{i}
+$$
+- Think of these as ``whatever observation $i$ contributes to the estimate of $\theta$''
+- Note when $\hat \theta$ is the sample mean, the pseudo observations are the data themselves
+- Then the sample standard error of these observations is the previous jackknife estimated standard error.
+- The mean of these observations is a bias-corrected estimate of $\theta$
+
+---
+
+## The bootstrap
+
+- The bootstrap is a tremendously useful tool for constructing confidence intervals and calculating standard errors for difficult statistics
+- For example, how would one derive a confidence interval for the median?
+- The bootstrap procedure follows from the so called bootstrap principle
+
+---
+
+## The bootstrap principle
+
+- Suppose that I have a statistic that estimates some population parameter, but I don't know its sampling distribution
+- The bootstrap principle suggests using the distribution defined by the data to approximate its sampling distribution
+
+---
+
+## The bootstrap in practice
+
+- In practice, the bootstrap principle is always carried out using simulation
+- We will cover only a few aspects of bootstrap resampling
+- The general procedure follows by first simulating complete data sets from the observed data with replacement
+
+ - This is approximately drawing from the sampling distribution of that statistic, at least as far as the data is able to approximate the true population distribution
+
+- Calculate the statistic for each simulated data set
+- Use the simulated statistics to either define a confidence interval or take the standard deviation to calculate a standard error
+
+---
+## Nonparametric bootstrap algorithm example
+
+- Bootstrap procedure for calculating confidence interval for the median from a data set of $n$ observations
+
+ i. Sample $n$ observations **with replacement** from the observed data resulting in one simulated complete data set
+
+ ii. Take the median of the simulated data set
+
+ iii. Repeat these two steps $B$ times, resulting in $B$ simulated medians
+
+ iv. These medians are approximately drawn from the sampling distribution of the median of $n$ observations; therefore we can
+
+ - Draw a histogram of them
+ - Calculate their standard deviation to estimate the standard error of the median
+ - Take the $2.5^{th}$ and $97.5^{th}$ percentiles as a confidence interval for the median
+
+---
+
+## Example code
+
+```{r}
+B <- 1000
+resamples <- matrix(sample(x,
+ n * B,
+ replace = TRUE),
+ B, n)
+medians <- apply(resamples, 1, median)
+sd(medians)
+quantile(medians, c(.025, .975))
+```
+
+---
+## Histogram of bootstrap resamples
+
+```{r, fig.height=5, fig.width=5}
+hist(medians)
+```
+
+---
+
+## Notes on the bootstrap
+
+- The bootstrap is non-parametric
+- Better percentile bootstrap confidence intervals correct for bias
+- There are lots of variations on bootstrap procedures; the book "An Introduction to the Bootstrap"" by Efron and Tibshirani is a great place to start for both bootstrap and jackknife information
+
+
+---
+## Group comparisons
+- Consider comparing two independent groups.
+- Example, comparing sprays B and C
+
+```{r, fig.height=4, fig.width=4}
+data(InsectSprays)
+boxplot(count ~ spray, data = InsectSprays)
+```
+
+---
+## Permutation tests
+- Consider the null hypothesis that the distribution of the observations from each group is the same
+- Then, the group labels are irrelevant
+- We then discard the group levels and permute the combined data
+- Split the permuted data into two groups with $n_A$ and $n_B$
+ observations (say by always treating the first $n_A$ observations as
+ the first group)
+- Evaluate the probability of getting a statistic as large or
+ large than the one observed
+- An example statistic would be the difference in the averages between the two groups;
+ one could also use a t-statistic
+
+---
+## Variations on permutation testing
+Data type | Statistic | Test name
+---|---|---|
+Ranks | rank sum | rank sum test
+Binary | hypergeometric prob | Fisher's exact test
+Raw data | | ordinary permutation test
+
+- Also, so-called *randomization tests* are exactly permutation tests, with a different motivation.
+- For matched data, one can randomize the signs
+ - For ranks, this results in the signed rank test
+- Permutation strategies work for regression as well
+ - Permuting a regressor of interest
+- Permutation tests work very well in multivariate settings
+
+---
+## Permutation test for pesticide data
+```{r}
+subdata <- InsectSprays[InsectSprays$spray %in% c("B", "C"),]
+y <- subdata$count
+group <- as.character(subdata$spray)
+testStat <- function(w, g) mean(w[g == "B"]) - mean(w[g == "C"])
+observedStat <- testStat(y, group)
+permutations <- sapply(1 : 10000, function(i) testStat(y, sample(group)))
+observedStat
+mean(permutations > observedStat)
+```
+
+---
+## Histogram of permutations
+```{r, echo= FALSE, fig.width=5, fig.height=5}
+hist(permutations)
+```
+
+
diff --git a/06_StatisticalInference/03_06_resampledInference/index.html b/06_StatisticalInference/03_06_resampledInference/index.html
index 15515b719..ce1da45ea 100644
--- a/06_StatisticalInference/03_06_resampledInference/index.html
+++ b/06_StatisticalInference/03_06_resampledInference/index.html
@@ -1,614 +1,609 @@
-
-
-
- Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
library(UsingR)
-data(father.son)
-x <- father.son$sheight
-n <- length(x)
-theta <- median(x)
-jk <- sapply(1 : n,
- function(i) median(x[-i])
- )
-thetaBar <- mean(jk)
-biasEst <- (n - 1) * (thetaBar - theta)
-seEst <- sqrt((n - 1) * mean((jk - thetaBar)^2))
-c(biasEst, seEst)
-[1] 0.0000 0.1014
-library(bootstrap)
-temp <- jackknife(x, median)
-c(temp$jack.bias, temp$jack.se)
-[1] 0.0000 0.1014
-The general procedure follows by first simulating complete data sets from the observed data with replacement
- -Calculate the statistic for each simulated data set
Use the simulated statistics to either define a confidence interval or take the standard deviation to calculate a standard error
Bootstrap procedure for calculating confidence interval for the median from a data set of \(n\) observations
- -i. Sample \(n\) observations with replacement from the observed data resulting in one simulated complete data set
- -ii. Take the median of the simulated data set
- -iii. Repeat these two steps \(B\) times, resulting in \(B\) simulated medians
- -iv. These medians are approximately drawn from the sampling distribution of the median of \(n\) observations; therefore we can
- -B <- 1000
-resamples <- matrix(sample(x,
- n * B,
- replace = TRUE),
- B, n)
-medians <- apply(resamples, 1, median)
-sd(medians)
-[1] 0.08546
-quantile(medians, c(.025, .975))
- 2.5% 97.5%
-68.43 68.82
-hist(medians)
-data(InsectSprays)
-boxplot(count ~ spray, data = InsectSprays)
-
| Data type | -Statistic | -Test name | -
|---|---|---|
| Ranks | -rank sum | -rank sum test | -
| Binary | -hypergeometric prob | -Fisher's exact test | -
| Raw data | -- | ordinary permutation test | -
subdata <- InsectSprays[InsectSprays$spray %in% c("B", "C"),]
-y <- subdata$count
-group <- as.character(subdata$spray)
-testStat <- function(w, g) mean(w[g == "B"]) - mean(w[g == "C"])
-observedStat <- testStat(y, group)
-permutations <- sapply(1 : 10000, function(i) testStat(y, sample(group)))
-observedStat
-[1] 13.25
-mean(permutations > observedStat)
-[1] 0
-
Brian Caffo, Jeff Leek, Roger Peng
Johns Hopkins Bloomberg School of Public Health
library(UsingR)
+data(father.son)
+x <- father.son$sheight
+n <- length(x)
+theta <- median(x)
+jk <- sapply(1:n, function(i) median(x[-i]))
+thetaBar <- mean(jk)
+biasEst <- (n - 1) * (thetaBar - theta)
+seEst <- sqrt((n - 1) * mean((jk - thetaBar)^2))
+c(biasEst, seEst)
+## [1] 0.0000 0.1014
+library(bootstrap)
+temp <- jackknife(x, median)
+c(temp$jack.bias, temp$jack.se)
+## [1] 0.0000 0.1014
+The general procedure follows by first simulating complete data sets from the observed data with replacement
+ +Calculate the statistic for each simulated data set
Use the simulated statistics to either define a confidence interval or take the standard deviation to calculate a standard error
Bootstrap procedure for calculating confidence interval for the median from a data set of \(n\) observations
+ +i. Sample \(n\) observations with replacement from the observed data resulting in one simulated complete data set
+ +ii. Take the median of the simulated data set
+ +iii. Repeat these two steps \(B\) times, resulting in \(B\) simulated medians
+ +iv. These medians are approximately drawn from the sampling distribution of the median of \(n\) observations; therefore we can
+ +B <- 1000
+resamples <- matrix(sample(x, n * B, replace = TRUE), B, n)
+medians <- apply(resamples, 1, median)
+sd(medians)
+## [1] 0.08834
+quantile(medians, c(0.025, 0.975))
+## 2.5% 97.5%
+## 68.41 68.82
+hist(medians)
+
data(InsectSprays)
+boxplot(count ~ spray, data = InsectSprays)
+
| Data type | +Statistic | +Test name | +
|---|---|---|
| Ranks | +rank sum | +rank sum test | +
| Binary | +hypergeometric prob | +Fisher's exact test | +
| Raw data | ++ | ordinary permutation test | +
subdata <- InsectSprays[InsectSprays$spray %in% c("B", "C"), ]
+y <- subdata$count
+group <- as.character(subdata$spray)
+testStat <- function(w, g) mean(w[g == "B"]) - mean(w[g == "C"])
+observedStat <- testStat(y, group)
+permutations <- sapply(1:10000, function(i) testStat(y, sample(group)))
+observedStat
+## [1] 13.25
+mean(permutations > observedStat)
+## [1] 0
+
Load the data set mtcars in the datasets R package. You want
+to test whether the MPG is \(\mu_0\) or smaller using a one sided
+5% level test. (\(H_0 : \mu = \mu_0\) versus \(H_a : \mu < \mu_0\)).
+Using that data set and a Z test:
Both to two decimal places.
+ + + + + + +This is the inversion of a one sided hypothesis test. It yields confidence +bounds. (Note inverting a two sidded test yields confidence intervals.) +Think about the derivation of the confidence interval.
+ +We want to solve +\[ +\frac{\sqrt{n}(\bar{X} - \mu_0)}{s} = Z_{0.05} +\] +Or \[\mu_0 = \bar{X} - Z_{0.05} s / \sqrt{n} = \bar{X} + Z_{0.95} s / \sqrt{n}\]. Note that the quantile is negative.
+ +mn <- mean(mtcars$mpg); s <- sd(mtcars$mpg); z <- qnorm(.05)
+mu0 <- mn - z * s / sqrt(nrow(mtcars))
+Note, it's easy to get tripped up in this problem on signs. If you get a value +that's less than \(\bar X\), then clearly it's wrong, since you'd never reject for +\(H_0:\mu = \mu_0\) versus \(H_a : \mu < \mu_0\) if \(\mu_0\) was less than your observed mean. +Also note the answer to "What is the largest value for which you would reject for?" is +infinity.
+ +21.84
+ +Consider again the mtcars dataset. Use a two group t-test to test
+the hypothesis that the 4 and 6 cyl cars have the same mpg. Use
+a two sided test with unequal variances.
Use t.test with the options var.equal=FALSE, paired=FALSE, altnernative as two.sided.
m4 <- mtcars$mpg[mtcars$cyl == 4]
+m6 <- mtcars$mpg[mtcars$cyl == 6]
+p <- t.test(m4, m6, paired = FALSE, alternative="two.sided", var.equal=FALSE)$p.value
+The answer to 1. is 1
+The answer to 2. is 4e-04
A sample of 100 men yielded an average PSA level of 3.0 with a sd of 1.1. What +are the complete set of values that a 5% two sided Z test of \(H_0 : \mu = \mu_0\) +would reject the null hypothesis for?
+ +This is equivalent to the confidence interval.
+ +The answer to 1 is
+ 2.78
+The answer to 2 is 2.78
You believe the coin that you're flipping is biased towards heads. You get 55 heads out of +100 flips.
+ +Use pbinom for a hypothesis that \(p=.5\) veruss \(p>.5\) where \(p\) is the binomial success
+probability.
Note you have to start at 54 as it lower.tail = FALSE gives the strictly greater than
+probabilities
ans <- round(pbinom(54, prob = .5, size = 100, lower.tail = FALSE),4)
+The answer to 1 is 0.1841
+The answer to 2 is 0
A web site was monitored for a year and it received 520 hits per day. In the first +30 days in the next year, the site received 15,800 hits. Assuming that web hits +are Poisson.
+ +Consider using ppois with \(\lambda=520 * 30\). Note this is nearly exactly Gaussian,
+so one could get away with the Gaussian calculation.
This test comes with the important assumption that web hits are a Poisson process.
+ +pv <- ppois(15800 - 1, lambda = 520 * 30, lower.tail = FALSE)
+The answer to 1 is 0.0553
+The answer to 2 is 0
Also, compare with the Gaussian approximation where \(\hat \lambda \sim N(\lambda, \lambda / t)\)
+ +pnorm(15800 / 30, mean = 520, sd = sqrt(520 / 30), lower.tail = FALSE)
+[1] 0.05466
+As \(t\rightarrow \infty\) this becomes more Gaussian. The approximation is pretty good for this +setting.
+ +Suppose that in an AB test, one advertising scheme led to an average of 10 purchases per day for a sample of 100 days, while the other led to 11 purchaces per day, also for a sample of 100 days. +Assuming a common standard deviation of 4 purchaces per day. +Assuming that the groups are independent and that they days are iid, perform a Z test of +equivalence.
+ +The standard error is +\[ +s \sqrt{\frac{1}{n_1} + \frac{1}{n_2}} +\]
+ +m1 <- 10; m2 <- 11
+n1 <- n2 <- 100
+s <- 4
+se <- s * sqrt(1 / n1 + 1 / n2)
+ts <- (m2 - m1) / se
+pv <- 2 * pnorm(-abs(ts))
+The answer to 1 is 0.077
+The answer to 2 is 0
A confidence interval for the mean contains:
+ +This is directly from the notes. Note that a confidence interval gives +values of \(\mu\) that are supported by the data whereas a test rejects +for values of \(\mu\) different from \(\mu_0\).
+ +The only complicated part of this is the 2. Note that a 95% interval corresponds +to a 5% level two sided test. So it's \(\alpha = 1 - Conf.Level\). The confusion is that +for both the two sided test and confidence interval, one needs to calculate +\(Z_{1 - \alpha / 2}\) (or the relevant T quantile).
+ +Consider two problems previous. Assuming that 10 purchases per day is a benchmark null value, +that days are iid and that the standard deviation is 4 purchases for day. Suppose that you +plan on sampling 100 days. What would be the power for a one sided 5% +Z mean test that purchases per day +have increased under the alternative of \(\mu = 11\) purchase per day?
+ +Under \(H_0\) \(\bar X \sim N(10, .4)\). Under \(H_a\) \(\bar X \sim N(11, .4)\). We reject +when \(\bar X \geq 10 + Z_{.95} .4\).
+ +The hint prettty much gives it away.
+ +power <- pnorm(10 + qnorm(.95) * .4, mean = 11, sd = .4, lower.tail = FALSE)
+The answer is 0.804
+ +Researchers would like to conduct a study of healthy adults to detect a four year mean brain volume loss of .01 mm3. Assume that the standard deviation of four year volume loss in this population is .04 mm3.
+ +Under \(H_0\) \(\bar X\) is \(N(0, .05 / \sqrt{n})\) and is \(N(.01, .05 / \sqrt{n})\) under \(H_a\). +We will reject if +\[ +\bar X \geq Z_{.95} s / sqrt{n} +\] +which has probability 0.05 under \(H_0\). Under \(H_a\) it has probability +\[ +P\left( \frac{\bar X - 0.01}{s / \sqrt{n}} \geq \frac{.01}{s / \sqrt{n}} + z_{.95} \right) += P\left( Z \geq \frac{.01}{s / \sqrt{n}} + z_{.95}\right) +\]
+ +Looking at the hint we set +\[ +\frac{.01}{s / \sqrt{n}} + z_{.95} = z_{.2} +\] +\[ +n = \frac{(z_{.95} - z_{.2})^2 s^2}{.01^2} = \frac{ (z_{.95} + z_{.8})^2 s^2}{.01^2} +\] +So we get
+ +n <- (qnorm(.95) + qnorm(.8)) ^ 2 * .04 ^ 2 / .01^2
+The answer is 99
+ +In a court of law, all things being equal, if via policy you require a lower +standard of evidence to convict people then
+ +Think about it.
+ +If you require less evidence to convict, then you will convict more people, guilty and +innocent. Relate this property back to hypothesis tests.
+ +Consider the mtcars data set.
mpg8 <- mtcars$mpg[mtcars$cyl == 8]
+mpg6 <- mtcars$mpg[mtcars$cyl == 6]
+m8 <- mean(mpg8); m6 <- mean(mpg6)
+s8 <- sd(mpg8); s6 <- sd(mpg6)
+n8 <- length(mpg8); n6 <- length(mpg6)
+p <- t.test(mpg8, mpg6, paired = FALSE, alternative="two.sided", var.equal=TRUE)$p.value
+mixprob <- (n8 - 1) / (n8 + n6 - 2)
+s <- sqrt(mixprob * s8 ^ 2 + (1 - mixprob) * s6 ^ 2)
+z <- (m8 - m6) / (s * sqrt(1 / n8 + 1 / n6))
+pz <- 2 * pnorm(-abs(z))
+## Hand calculating the T just to check
+#2 * pt(-abs(z), df = n8 + n6 - 2)
+The Bonferonni correction controls this
+ +This is pretty much straight out of the notes
+ +The Bonferonni correction is the classic correction for the familywise error rate.
+ +
Johns Hopkins Bloomberg School of Public Health
Including math in a sentence involves wrapping the symbols in $ symbols. For -example if you write this in an R markdown file:
- -$\hat{\alpha} = 4$" Then you get the following text after running knit2html or slidify on your document.
- -Sometimes you have several equations you would like to line up. The -way that you do that is with a double dollar sign $$ and the align command.
- -For example if you write
- -
Then you get the following text after running knit2html or slidify on your document.
- -\[ - \begin{aligned} - y &= \beta_0 + \beta_1 + x_1 + \epsilon\\ - x &= \gamma z \\ - z &\sim N(0,1) - \end{aligned} -\]
- -$a_{b}$$a^{b}$$\alpha, \beta, \ldots$$\sum_{n=1}^N$$\times$$\prod_{n=1}^N$$<, \leq, \geq$$\sim$$\widehat{\alpha}$$\bar{x}$$\frac{a}{b}$$\left(\frac{a}{b}\right)$
-
-Then you get the following text after running _knit2html_ or _slidify_ on your document.
-
-$$
- \begin{aligned}
- y &= \beta_0 + \beta_1 + x_1 + \epsilon\\
- x &= \gamma z \\
- z &\sim N(0,1)
- \end{aligned}
-$$
-
----
-
-## Common symbols
-
-* Subscripts to get $a_{b}$ write: `$a_{b}$`
-* Superscripts write $a^{b}$ write: `$a^{b}$`
-* Greek letters like $\alpha, \beta, \ldots$ write: `$\alpha, \beta, \ldots$`
-* Sums like $\sum_{n=1}^N$ write: `$\sum_{n=1}^N$`
-* Multiplication like $\times$ write: `$\times$`
-* Products like $\prod_{n=1}^N$ write: `$\prod_{n=1}^N$`
-* Inequalities like $<, \leq, \geq$ write: `$<, \leq, \geq$`
-* Distributed like $\sim$ write: `$\sim$`
-* Hats like $\widehat{\alpha}$ write: `$\widehat{\alpha}$`
-* Averages like $\bar{x}$ write: `$\bar{x}$`
-* Fractions like $\frac{a}{b}$ write: `$\frac{a}{b}$`
-* Big parentheses like $\left(\frac{a}{b}\right)$ write: `$\left(\frac{a}{b}\right)$`
-
-
----
-
-## For more information
-
-* Rstudio's equations page:
- * http://www.rstudio.com/ide/docs/authoring/using_markdown_equations
-* Lists of Latex symbols
- * http://www.rpi.edu/dept/arc/training/latex/LaTeX_symbols.pdf
- * http://www.giss.nasa.gov/tools/latex/ltx-2.html
- * http://omega.albany.edu:8008/Symbols.html
-
-
-
diff --git a/09_DevelopingDataProducts/rStudioPresent/index.md b/09_DevelopingDataProducts/rStudioPresent/index.md
index b998542ae..6b5c3334e 100644
--- a/09_DevelopingDataProducts/rStudioPresent/index.md
+++ b/09_DevelopingDataProducts/rStudioPresent/index.md
@@ -1,7 +1,7 @@
RStudio Presenter
===
author: Brian Caffo, Jeff Leek Roger Peng
-date: May 21 2014
+date: May 24 2014
transition: rotate