--- title: "STAT 801A Lab 4: The t Distribution and Power" author: "Insert Name" date: today format: html editor: visual --- ```{r} library(ggplot2) ``` # Inference When $\sigma$ Is Unknown In Lab 3 we assumed the population standard deviation $\sigma$ was known. That kept the math simple, but is almost never true. If we do not know $\mu$, we almost never know $\sigma$ either. Both are parameters of the population, and we only have a sample. The fix looks obvious, replace $\sigma$ with the sample standard deviation $s$. $$SE = \frac{s}{\sqrt{n}}$$ Previously, the only thing that varied from sample to sample was $\bar{y}$. Now $s$ varies too, so the test statistic has two moving parts instead of one. That extra uncertainty means the statistic is no longer normal. ## The sample standard deviation $$s = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(y_i - \bar{y})^2}$$ In R this is `sd()`. ```{r} weights <- c(338.1, 341.2, 339.5, 342.0, 337.8, 340.6, 336.9, 343.1, 339.9, 338.4, 341.7, 340.2, 337.5, 342.8, 339.1, 340.8, 338.6, 341.4, 336.2, 340.0) mean(weights) sd(weights) ``` ------------------------------------------------------------------------ # The t Distribution Replacing $\sigma$ with $s$ changes the distribution of the test statistic. $$\frac{\bar{y} - \mu}{s / \sqrt{n}} \sim t_{n-1}$$ The t distribution is centered at zero and symmetric, like the standard normal, but with heavier tails. Extreme values are more likely under t than under $N(0,1)$, which is the extra uncertainty from estimating $s$. ## Degrees of freedom The t distribution has one parameter, the degrees of freedom. For a single mean $$df = n - 1$$ We lose one degree of freedom because we must estimate $\bar{y}$. ## The t functions in R The prefix system from previous labs is the same. The root is `t`. | Function | What it gives you | |-------------|-------------------------------------------------------------| | `dt(x, df)` | height of the curve | | `pt(x, df)` | the probability of that value or anything smaller | | `qt(x, df)` | works backwards: you give a probability, it returns a value | | `rt(x, df)` | random values drawn from the distribution | qnorm(0.975) ```{r} qt(0.975, df = 15) qnorm(0.975) ``` The t critical value is larger. That is the result of not knowing $\sigma$. ## Comparing t and the standard normal ```{r} x <- seq(-4, 4, length.out = 1000) curves <- data.frame( x = rep(x, 2), density = c(dnorm(x), dt(x, df = 3)), curve = rep(c("Normal", paste("t, df =", 3)), each = length(x))) ggplot(curves, aes(x = x, y = density, color = curve)) + geom_line(linewidth = 1) + labs(x = "Value", y = "Density", color = "") ``` As $n$ grows, $s$ becomes a better estimate of $\sigma$, so the extra uncertainty shrinks and t converges to $N(0,1)$. ------------------------------------------------------------------------ # Confidence Intervals with t Same structure, with two substitutions: $s$ for $\sigma$, and a t critical value for a z critical value. $$\bar{y} \pm t_{\alpha/2, \, n-1} \frac{s}{\sqrt{n}}$$ Note: a 95% interval puts 2.5% in each tail, so the critical value is `qt(0.975, df)`, not `qt(0.95, df)`. ## Computing one in R Return to the coffee weights, this time without pretending to know $\sigma$. ```{r} n <- length(weights) ybar <- mean(weights) s <- sd(weights) n ybar s ``` ```{r} se <- s / sqrt(n) crit <- qt(0.975,df = n - 1) margin <- crit * se se crit margin ``` ```{r} ybar - margin ybar + margin ``` We are 95% confident the true mean fill weight is between 338.88 and 340.70 grams. The t interval will be slightly wider than one calculated using the standard normal ------------------------------------------------------------------------ # The t Test The test statistic is the same idea as the z statistic, with $s$ in the denominator. $$t = \frac{\bar{y} - \mu_0}{s / \sqrt{n}}$$ When $H_0$ is true, this follows a t distribution with $n - 1$ degrees of freedom. ```{r} mu0 <- 340 t_stat <- (ybar - mu0) / (s / sqrt(n)) t_stat ``` The critical values and the p-value now come from `qt` and `pt` instead of `qnorm` and `pnorm`. ```{r} qt(0.025, df = n - 1) qt(0.975, df = n - 1) 2 * pt(abs(t_stat), df = n - 1, lower.tail = FALSE) #Two tailed p value ``` Our t is -0.48, nowhere near the rejection region, and the p-value is 0.64. We do not reject. ## The t.test function R packages all of that into one function. ```{r} t.test(weights, mu = 340) ``` Everything matches what we computed by hand ## Other arguments Use `alternative` for a one-tailed test and `conf.level` to change the confidence level. ```{r} t.test(weights, mu = 340, alternative = "less") ``` ```{r} t.test(weights, mu = 340, conf.level = 0.99) ``` ------------------------------------------------------------------------ # Example A seed company advertises that its new corn hybrid averages 192 bushels per acre. A grower cooperative suspects the real mean is lower than advertised, and plants the hybrid in 18 test plots. ```{r} yields <- c(182.4, 191.7, 178.9, 195.2, 187.6, 183.1, 199.4, 176.8, 190.3, 185.5, 193.8, 181.2, 188.9, 197.1, 179.6, 186.4, 192.5, 184.7) ``` ## Step 1: State the hypotheses The cooperative only cares about falling short, so this is one-tailed in the lower direction. $$H_0: \mu = 192 \qquad H_a: \mu < 192$$ ## Look at the data Always look before testing. ```{r} length(yields) mean(yields) sd(yields) ``` ```{r} ggplot(data.frame(yield = yields), aes(x = yield)) + geom_histogram(bins = 10, fill = "lightblue", color = "black") + geom_vline(xintercept = 192, linetype = "dashed") + labs(x = "Yield (bushels per acre)", y = "Count") ``` The sample mean is 187.5, below the dashed line at 192. Most plots fall short of the advertised value, but several exceed it. The data is also not strongly skewed in one direction. We can use a t test to see if the real mean yield is lower than advertised. ## Compute the test statistic ```{r} n_corn <- length(yields) ybar_corn <- mean(yields) s_corn <- sd(yields) t_corn <- (ybar_corn - 192) / (s_corn / sqrt(n_corn)) t_corn ``` ## Find the rejection region One-tailed at $\alpha = 0.05$, with all 5% in the lower tail ```{r} qt(0.05, df = n_corn - 1) ``` We reject if $t < -1.74$. Our t of -2.90 is below that, so we reject $H_0$. ## Compute the p-value For a lower-tailed test the p-value is the area to the left of our statistic, so we use `pt` with its default `lower.tail = TRUE`. ```{r} pt(t_corn, df = n_corn - 1) ``` The p-value is 0.005. If the true mean really were 192, a sample this low would turn up about five times in a thousand. ## Report a confidence interval ```{r} crit_corn <- qt(0.975, df = n_corn - 1) margin_corn <- crit_corn * (s_corn / sqrt(n_corn)) ybar_corn - margin_corn ybar_corn + margin_corn ``` We are 95% confident the true mean yield is between 184.2 and 190.8 bushels per acre ## Confirm with t.test ```{r} t.test(yields, mu = 192,alternative = "less") ``` Note: a one-tailed t.test returns a one-sided interval, so one endpoint is Inf. The test only looks in one direction, so the interval only bounds the mean from one side. Drop `alternative` to get the two-sided interval. ## Write the conclusion There is evidence that the new hybrid averages less than the advertised 192 bushels per acre ($t = -2.90$, $df = 17$, $p = 0.005$). The estimated mean yield is 187.5 bushels per acre, with a 95% confidence interval of 184.2 to 190.8. ------------------------------------------------------------------------ # Power A hypothesis test can fail in two directions. | | $H_0$ is true | $H_0$ is false | |-------------------|-----------------------------|------------------------------| | We reject $H_0$ | Type I error, rate $\alpha$ | Correct, probability = power | | We fail to reject | Correct | Type II error, rate $\beta$ | Type I error: we set $\alpha$ in advance and accept that rate of false positives. Power is the other corner, the probability of detecting a real effect. $$\text{power} = 1 - \beta = P(\text{reject } H_0 \mid H_0 \text{ is false})$$ Power is not one number for a test. It depends on how false $H_0$ is: detecting a huge effect is easy, detecting a tiny one is hard. ## Seeing what power is So far every test used one curve, the null distribution. Power needs a second curve for what happens when a specific alternative is true instead. Back to the corn example: $\mu = 192$ against $\mu < 192$, 18 plots, and suppose the truth is 188 with a SD of about 6.5. The gap between the truth and the null, in standard errors, is the noncentrality parameter. $$ncp = \frac{\mu_{true} - \mu_0}{SD / \sqrt{n}} = \frac{188 - 192}{6.5 / \sqrt{18}} = -2.61$$ ```{r} se_power <- 6.5 / sqrt(18) ncp_power <- (188 - 192) / se_power crit_power <- qt(0.05, df = 17) se_power ncp_power crit_power ``` The null curve is the usual $t_{17}$. The alternative curve is a t shifted to the noncentrality parameter, which `dt` draws when given `ncp`. ```{r} curves_power <- data.frame(t = seq(-7, 4, length.out = 1000)) curves_power$null_density <- dt(curves_power$t, df = 17) curves_power$alt_density <- dt(curves_power$t, df = 17, ncp = ncp_power) ggplot(curves_power, aes(x = t)) + geom_area(data = subset(curves_power, t <= crit_power), aes(y = alt_density), fill = "lightblue") + geom_line(aes(y = null_density), linewidth = 1) + geom_line(aes(y = alt_density), linewidth = 1, linetype = "dashed") + geom_vline(xintercept = crit_power, linetype = "dotted") + labs(x = "t statistic", y = "Density") ``` Solid is the null, dashed is the truth, and the dotted line is the critical value of -1.74. We reject anything to its left. The shaded area is power: the share of the dashed curve falling in the rejection region. What is left unshaded is $\beta$. ```{r} pt(crit_power, df = 17, ncp = ncp_power) ``` About 0.80, so four studies in five would catch this effect. ## What changes power Four things move power, and it helps to know which ones you control. | Factor | Effect on power | Under your control? | |--------------------|-----------------|--------------------------------------| | Larger effect size | Higher | No | | Larger $n$ | Higher | Yes | | Larger $\alpha$ | Higher | Yes, but also raises false positives | | Larger $\sigma$ | Lower | Rarely | All four work through the picture. Effect size and $n$ push the dashed curve away from the solid one, $\alpha$ slides the dotted line, and $\sigma$ pulls the curves back together. ## power.t.test R computes power directly. ```{r} power.t.test(n = 18, delta = 4, sd = 6.5, sig.level = 0.05, type = "one.sample", alternative = "one.sided") ``` The same 0.80 represented by the blue shaded area `delta` is the size of the difference to detect, not the true mean, and it is always positive For planning, leave `n` out and supply the power you want instead. The function will tell you what sample size is required to achieve the given power (round up). ```{r} power.t.test(delta = 4, sd = 6.5, sig.level = 0.05, power = 0.90, type = "one.sample", alternative = "one.sided") ``` ------------------------------------------------------------------------ # Practice Problem A sporting goods company measures the radius, in centimeters, of 36 volleyballs coming off a production line. Official FIVB rules require a circumference of 65 to 67 cm. Taking the midpoint of 66 cm, the regulation radius is $$r = \frac{66}{2\pi} \approx 10.5 \text{ cm}$$ so the production line should be centered on $\mu = 10.5$ cm. The population standard deviation is unknown. \(a\) Import `volleyball.csv` and assign it a logical name. ```{r} volleyball <- read.csv("volleyball.csv") ``` \(b\) Compute the sample size, sample mean, sample standard deviation, and standard error. ```{r} n_vb <- nrow(volleyball) ybar_vb <- mean(volleyball$radius) s_vb <- sd(volleyball$radius) se_vb <- s_vb / sqrt(n_vb) n_vb ybar_vb s_vb se_vb ``` \(c\) Make a histogram of the radii with a dashed line at 10.5 cm. ```{r} ggplot(volleyball, aes(x = radius)) + geom_histogram(bins = 10, fill = "lightblue", color = "black") + geom_vline(xintercept = 10.5, linetype = "dashed") + labs(x = "Radius (cm)", y = "Count") ``` \(d\) Find the critical value for a 95% confidence interval. How many degrees of freedom are there, and why is this a t critical value rather than a z? ```{r} qt(0.975, df = n_vb - 1) ``` 2.03, with $n - 1 = 35$ degrees of freedom. We use t because $\sigma$ is unknown and estimated by $s$. \(e\) Construct the 95% confidence interval by hand. ```{r} margin_vb <- qt(0.975, df = n_vb - 1) * se_vb ybar_vb - margin_vb ybar_vb + margin_vb ``` \(f\) Interpret the interval in a sentence. We are 95% confident the true mean radius is between 11.16 and 11.84 cm. \(g\) State the hypotheses for testing whether the mean radius differs from the regulation value. $$H_0: \mu = 10.5 \qquad H_a: \mu \neq 10.5$$ \(h\) Compute the test statistic by hand. ```{r} t_vb <- (ybar_vb - 10.5) / se_vb t_vb ``` \(i\) Find the rejection region at $\alpha = 0.05$ and state your conclusion. ```{r} qt(0.025, df = n_vb - 1) qt(0.975, df = n_vb - 1) ``` Reject if $t < -2.03$ or $t > 2.03$. Our $t = 6.0$, so we reject $H_0$: the mean radius differs from 10.5 cm. \(j\) Compute the two-tailed p-value. ```{r} 2 * pt(abs(t_vb), df = n_vb - 1, lower.tail = FALSE) ``` \(k\) Confirm your work with `t.test`. ```{r} t.test(volleyball$radius, mu = 10.5) ``` \(l\) The plant halts production if there is evidence the mean radius has drifted above 11.25 cm. State the hypotheses for this one-tailed test and run it with `t.test`. $$H_0: \mu = 11.25 \qquad H_a: \mu > 11.25$$ ```{r} t.test(volleyball$radius, mu = 11.25, alternative = "greater") ``` \(m\) Based on the p-value from part (l), do we reject at $\alpha = 0.05$? $t = 1.5$ and $p = 0.071 > 0.05$, so we fail to reject. There is not enough evidence the mean has drifted above 11.25 cm. \(n\) Use `power.t.test` to find the power of the test in part (m), assuming the true mean really is 11.5 and using $s = 1$ as an estimate of $\sigma$. The difference you are trying to detect is 0.25 cm. ```{r} power.t.test(n = 36, delta = 0.25, sd = 1, sig.level = 0.05, type = "one.sample", alternative = "one.sided") ``` Power is about 0.43, so this test would miss a real 0.25 cm drift more often than it catches one. \(o\) Use `power.t.test` to find the sample size needed for 90% power to detect that same 0.25 cm difference. Remember to round up. How does it compare to the 36 balls actually measured? ```{r} power.t.test(delta = 0.25, sd = 1, sig.level = 0.05, power = 0.90, type = "one.sample", alternative = "one.sided") ``` About 138.4, so 139 balls. Nearly four times the 36 measured. # Final Step Render your document. Click the Render button at the top of the editor, or press Ctrl + Shift + K. If the render fails, read the error message. It names the chunk that caused the problem. The most common causes are a missing closing parenthesis, or an object used before it was created.